1 #include "FilteredAnitaEvent.h" 3 #include "AntennaPositions.h" 10 #include "AnitaGeomTool.h" 11 #include "Correlator.h" 12 #include "AnalysisWaveform.h" 13 #include "TStopwatch.h" 15 #ifdef UCORRELATOR_OPENMP 18 #define SECTIONS _Pragma("omp parallel sections") 19 #define SECTION _Pragma("omp section") 23 #define SECTIONS if(true) 24 #define SECTION if(true) 35 #ifdef UCORRELATOR_OPENMP 36 omp_lock_t waveform_locks[NANTENNAS];
37 omp_lock_t correlation_locks[NANTENNAS][NANTENNAS];
43 for (
int i = 0; i < NANTENNAS; i++)
45 omp_init_lock(&waveform_locks[i]);
47 for (
int j = 0; j < NANTENNAS; j++)
49 omp_init_lock(&correlation_locks[i][j]);
58 for (
int i = 0; i < NANTENNAS; i++)
60 omp_destroy_lock(&waveform_locks[i]);
62 for (
int j = 0; j < NANTENNAS; j++)
64 omp_destroy_lock(&correlation_locks[i][j]);
75 #ifdef UCORRELATOR_OPENMP 76 return omp_get_max_threads();
85 static int get_thread_id()
87 #ifdef UCORRELATOR_OPENMP 88 return omp_get_thread_num();
97 static int count_the_correlators = 1;
99 static int count_the_zoomed_correlators = 1;
104 template<
class T,
class A>
105 static void combineHists(
int N, T ** hists )
108 int nbins = (hists[0]->GetNbinsX() +2) * (hists[0]->GetNbinsY()+2);
109 __restrict A* sum = hists[0]->GetArray();
110 for (
int i = 1; i <N; i++ )
112 __restrict A* v = hists[i]->GetArray();
114 for(
int j = 0; j < nbins; j++)
124 #define RAD2DEG (180/M_PI) 127 #define PHI_SECTOR_ANGLE (360. / NUM_PHI) 131 UCorrelator::Correlator::Correlator(
int nphi,
double phi_min,
double phi_max,
int ntheta,
double theta_min,
double theta_max,
bool use_center,
bool scale_by_cos_theta,
double baseline_weight,
double gain_sigma)
132 : scale_cos_theta(scale_by_cos_theta) , baselineWeight(baseline_weight), gainSigma(gain_sigma)
134 TString histname = TString::Format(
"ucorr_corr_%d",count_the_correlators);
135 TString normname = TString::Format(
"ucorr_norm_%d",count_the_correlators++);
137 hist =
new TH2D(histname.Data(),
"Correlator", nphi, phi_min, phi_max, ntheta, theta_min, theta_max);
138 hist->SetDirectory(0);
139 norm =
new TH2D(normname.Data(),
"Normalization", nphi, phi_min, phi_max, ntheta, theta_min, theta_max);
140 norm->SetDirectory(0);
142 hist->GetXaxis()->SetTitle(
"#phi");
143 hist->GetYaxis()->SetTitle(
"-#theta");
145 hists.resize(nthreads());
146 norms.resize(nthreads());
151 for (
int i= 1; i < nthreads();i++)
153 hists[i] = (TH2D*) hist->Clone(TString(hist->GetName()) + TString::Format(
"_%d",i));
154 hists[i]->SetDirectory(0);
155 norms[i] = (TH2D*) norm->Clone(TString(norm->GetName()) + TString::Format(
"_%d",i));
156 norms[i]->SetDirectory(0);
162 use_bin_center = use_center;
163 memset(trigcache,0,
sizeof(trigcache));
166 disallowed_antennas = 0;
169 max_phi2 = max_phi*max_phi;
174 #ifdef UCORRELATOR_OPENMP 181 static int allowedPhisPairOfAntennas(
double &lowerAngle,
double &higherAngle,
double ¢erTheta1,
double ¢erTheta2,
double ¢erPhi1,
double ¢erPhi2,
int ant1,
int ant2,
double max_phi,
AnitaPol::AnitaPol_t pol)
188 int upperlimit=phi2+2;
189 int lowerlimit=phi2-2;
191 if(upperlimit>NUM_PHI-1)upperlimit-=NUM_PHI;
192 if(lowerlimit<0)lowerlimit+=NUM_PHI;
194 if (upperlimit>lowerlimit){
195 if (phi1<=upperlimit && phi1>=lowerlimit){
199 if (upperlimit<lowerlimit){
200 if (phi1<=upperlimit || phi1>=lowerlimit){
206 double centerAngle1, centerAngle2;
212 centerAngle1=ap->
phiAnt[pol][ant1];
213 centerAngle2=ap->
phiAnt[pol][ant2];
216 if (centerAngle2>centerAngle1)
218 lowerAngle=centerAngle2-max_phi;
219 higherAngle=centerAngle1+max_phi;
223 lowerAngle=centerAngle1-max_phi;
224 higherAngle=centerAngle2+max_phi;
227 if (lowerAngle<0) lowerAngle+=360;
228 if (higherAngle>360) higherAngle-=360;
239 centerPhi1=centerAngle1;
240 centerPhi2=centerAngle2;
246 void UCorrelator::Correlator::reset()
248 for (
int i = 0; i < NANTENNAS; i++)
250 if (padded_waveforms[i])
252 delete padded_waveforms[i];
253 padded_waveforms[i] = 0;
256 for (
int j = 0; j < NANTENNAS; j++)
258 if (correlations[i][j])
260 delete correlations[i][j];
261 correlations[i][j] = 0;
272 AnalysisWaveform * UCorrelator::Correlator::getCorrelation(
int ant1,
int ant2)
277 #ifdef UCORRELATOR_OPENMP 278 #ifndef FFTTOOLS_COMPILED_WITH_OPENMP 279 #pragma omp critical (getCorrelation) 285 #ifdef UCORRELATOR_OPENMP 286 omp_set_lock(&locks->waveform_locks[ant1]);
288 if (!padded_waveforms[ant1])
291 padded_waveforms[ant1] =
new AnalysisWaveform(*ev->getFilteredGraph(ant1, pol));
292 rms[ant1] = padded_waveforms[ant1]->
even()->GetRMS(2);
293 padded_waveforms[ant1]->padEven(1);
295 #ifdef UCORRELATOR_OPENMP 296 omp_unset_lock(&locks->waveform_locks[ant1]);
297 omp_set_lock(&locks->waveform_locks[ant2]);
300 if (!padded_waveforms[ant2])
303 padded_waveforms[ant2] =
new AnalysisWaveform(*ev->getFilteredGraph(ant2, pol));
305 rms[ant2] = padded_waveforms[ant2]->
even()->GetRMS(2);
306 padded_waveforms[ant2]->padEven(1);
309 #ifdef UCORRELATOR_OPENMP 310 omp_unset_lock(&locks->waveform_locks[ant2]);
311 omp_set_lock(&locks->correlation_locks[ant1][ant2]);
314 if (!correlations[ant1][ant2])
316 #ifdef UCORRELATOR_OPENMP 317 omp_set_lock(&locks->waveform_locks[ant1]);
318 omp_set_lock(&locks->waveform_locks[ant2]);
323 #ifdef UCORRELATOR_OPENMP 324 omp_unset_lock(&locks->waveform_locks[ant2]);
325 omp_unset_lock(&locks->waveform_locks[ant1]);
330 #ifdef UCORRELATOR_OPENMP 331 omp_unset_lock(&locks->correlation_locks[ant1][ant2]);
332 #ifndef FFTTOOLS_COMPILED_WITH_OPENMP 336 return correlations[ant1][ant2];
346 fprintf(stderr,
"Must call Correlator::compute() prior to Correlator::computeZoomed!!!!");
350 double phi0 = phi - dphi * nphi/2;
351 double phi1 = phi + dphi * nphi/2;
352 double theta0 = theta - dtheta * ntheta/2;
353 double theta1 = theta + dtheta * ntheta/2;
354 TH2D* zoomed_norm =
new TH2D(TString::Format(
"zoomed_norm_%d",count_the_zoomed_correlators),
"Zoomed Correlation Normalization",
356 ntheta, theta0, theta1);
357 zoomed_norm->SetDirectory(0);
361 answer->SetBins(nphi, phi0, phi1,
362 ntheta, theta0, theta1);
367 answer =
new TH2D(TString::Format(
"zoomed_corr_%d", count_the_zoomed_correlators++),
"Zoomed Correlation",
369 ntheta, theta0, theta1);
370 answer->SetDirectory(0);
373 answer->GetXaxis()->SetTitle(
"#phi");
374 answer->GetYaxis()->SetTitle(
"-#theta");
381 memset(closest,0,
sizeof(closest));
385 TrigCache cache(nphi, dphi, phi0, ntheta,dtheta,theta0, ap,
true,nant, nant ? closest : 0);
387 int n2loop = nant ? nant : NANTENNAS;
389 std::vector<std::pair<int,int> > pairs;
390 pairs.reserve(n2loop);
392 double center_point[2];
393 center_point[0] = phi;
394 center_point[1] = theta;
396 for (
int ant_i = 0; ant_i < n2loop; ant_i++)
398 int ant1 = nant ? closest[ant_i] : ant_i;
399 if (!nant && disallowed_antennas & (1ul << ant1))
continue;
400 for (
int ant_j = ant_i +1; ant_j < n2loop; ant_j++)
402 int ant2 = nant ? closest[ant_j] : ant_j;
403 if (!nant && disallowed_antennas & (1ul << ant2))
continue;
405 pairs.push_back(std::pair<int,int>(ant1,ant2));
409 unsigned nit = pairs.size();
414 std::vector<TH2D*> zoomed_hists(nthreads());
415 std::vector<TH2D*> zoomed_norms(nthreads());
417 zoomed_hists[0] = answer;
418 zoomed_norms[0] = zoomed_norm;
421 for (
int i= 1; i < nthreads();i++)
423 zoomed_hists[i] = (TH2D*) answer->Clone(TString(answer->GetName()) + TString::Format(
"_%d",i));
424 zoomed_hists[i]->SetDirectory(0);
425 zoomed_norms[i] = (TH2D*) zoomed_norm->Clone(TString(norm->GetName()) + TString::Format(
"_%d",i));
426 norms[i]->SetDirectory(0);
431 #ifdef UCORRELATOR_OPENMP 432 #pragma omp parallel for 434 for (
unsigned it = 0; it < nit; it++)
436 doAntennas(pairs[it].first, pairs[it].second, &zoomed_hists[0], &zoomed_norms[0], &cache, center_point);
445 combineHists<TH2D,double>(nthreads(), &zoomed_hists[0]);
448 combineHists<TH2D,double>(nthreads(), &zoomed_norms[0]);
451 for (
int i =1; i < nthreads(); i++)
453 delete zoomed_hists[i];
454 delete zoomed_norms[i];
460 for (
int i = 0; i < (answer->GetNbinsX()+2) * (answer->GetNbinsY()+2); i++)
462 double val = answer->GetArray()[i];
463 if (val == 0)
continue;
464 double this_norm = zoomed_norm->GetArray()[i];
465 answer->GetArray()[i] = this_norm > 0 ? val/this_norm : 0;
469 answer->SetEntries(nonzero);
476 static inline bool between(
double phi,
double low,
double high)
479 double diff_highlow = fmod(high - low + 360, 360);
480 double diff_philow = fmod(phi - low + 360, 360);
482 return diff_philow < diff_highlow;
486 inline void UCorrelator::Correlator::doAntennas(
int ant1,
int ant2, TH2D ** these_hists,
488 const double * center_point )
491 double lowerAngleThis, higherAngleThis, centerTheta1, centerTheta2, centerPhi1, centerPhi2;
493 allowedFlag=allowedPhisPairOfAntennas(lowerAngleThis,higherAngleThis,
494 centerTheta1, centerTheta2, centerPhi1, centerPhi2,
495 ant1,ant2, max_phi, pol);
500 if(!allowedFlag)
return;
502 TH2D * the_hist = these_hists[get_thread_id()];
503 TH2D * the_norm = these_norms[get_thread_id()];
507 if (center_point && !between(center_point[0], lowerAngleThis, higherAngleThis))
return;
511 int nphibins = the_hist->GetNbinsX() + 2;
515 int first_phi_bin = center_point ? 1 : the_hist->GetXaxis()->FindFixBin(lowerAngleThis);
516 int last_phi_bin = center_point ? the_hist->GetNbinsX() : the_hist->GetXaxis()->FindFixBin(higherAngleThis);
518 if (first_phi_bin == 0) first_phi_bin = 1;
519 if (last_phi_bin == the_hist->GetNbinsX()+1) last_phi_bin = the_hist->GetNbinsX();
520 bool must_wrap = (last_phi_bin < first_phi_bin) ;
526 int maxsize = the_hist->GetNbinsY() * the_hist->GetNbinsX();
528 double phi_diff1 =
FFTtools::wrap(centerPhi1 - centerPhi2, 360, 0);
529 double phi_diff2 =
FFTtools::wrap(centerPhi2 - centerPhi1, 360, 0);
530 double baseline_phi = (fabs(phi_diff1) < fabs(phi_diff2)) ? centerPhi2 + phi_diff1 / 2 : centerPhi1 + phi_diff2 / 2;
533 double baseline_theta = (centerTheta1 + centerTheta2) / 2;
537 int * alloc =
new int[3*maxsize];
538 int * phibins = alloc;
539 int * thetabins = alloc + maxsize;
540 int * bins_to_fill = alloc + 2 *maxsize;
544 double * gain_weights = 0;
545 if (gainSigma && !center_point) gain_weights =
new double[maxsize];
548 for (
int phibin = first_phi_bin; (phibin <= last_phi_bin) || must_wrap; phibin++)
550 if (must_wrap && phibin == nphibins-1)
556 double phi = cache->phi[phibin-1] ;
558 double dphi1 = center_point ? 0 :
FFTtools::wrap(phi - centerPhi1,360,0);
559 double dphi2 = center_point ? 0 :
FFTtools::wrap(phi - centerPhi2,360,0);
563 if (!center_point && fabs(dphi1) > max_phi)
continue;
564 if (!center_point && fabs(dphi2) > max_phi)
continue;
566 int ny = the_hist->GetNbinsY();
569 for (
int thetabin = 1; thetabin <= ny; thetabin++)
571 double theta = cache->theta[thetabin-1];
573 double dtheta1 = center_point ? 0 :
FFTtools::wrap(theta - centerTheta1, 360, 0);
574 double dtheta2 = center_point ? 0 :
FFTtools::wrap(theta - centerTheta2, 360, 0);
577 if (!center_point && dphi1*dphi1 + dtheta1*dtheta1 > max_phi * max_phi)
continue;
578 if (!center_point && dphi2*dphi2 + dtheta2*dtheta2 > max_phi * max_phi)
continue;
580 phibins[nbins_used] = phibin;
581 thetabins[nbins_used] = thetabin;
582 bins_to_fill[nbins_used] = phibin + thetabin * nphibins;
583 if (gainSigma && !center_point)
588 double Dtheta = theta - baseline_theta;
590 gain_weights[nbins_used] = TMath::Gaus(TMath::Sqrt(Dphi * Dphi + Dtheta * Dtheta), 0, gainSigma);
597 double * dalloc =
new double[2 *nbins_used];
598 double * __restrict__ vals_to_fill = dalloc;
599 double * __restrict__ times_to_fill = dalloc + nbins_used;
602 for (
int i = 0; i < nbins_used; i++)
605 int phibin = phibins[i];;
606 int thetabin = thetabins[i];
607 times_to_fill[i] =
getDeltaTFast(ant1, ant2, phibin-1, thetabin-1,pol,cache, groupDelayFlag);
610 correlation->
evalEven(nbins_used, times_to_fill, vals_to_fill);
615 for(
int i = 0; i < nbins_used; i++)
617 vals_to_fill[i] *= cache->cos_theta[thetabins[i]-1];
624 double wgt = TMath::Power(cache->ap->distance(ant1, ant2, pol), baselineWeight);
625 for (
int i = 0; i < nbins_used; i++)
627 vals_to_fill[i] *= wgt;
633 double * __restrict__ the_arr = the_hist->GetArray();
634 double * __restrict__ the_norm_arr = the_norm->GetArray();
636 for (
int bi = 0; bi < nbins_used; bi++)
638 double val = vals_to_fill[bi];
639 int bin = bins_to_fill[bi];
641 if (gainSigma && !center_point)
643 the_arr[bin] += val * gain_weights[bi];
644 the_norm_arr[bin] += gain_weights[bi];
649 the_norm_arr[bin]+= 1;
656 if (gain_weights)
delete [] gain_weights;
675 int v =
event->getAnitaVersion();
676 AnitaVersion::set(v);
680 int nphi = hist->GetNbinsX();
681 double phi_max = hist->GetXaxis()->GetXmax();
682 double phi_min = hist->GetXaxis()->GetXmin();
683 int ntheta = hist->GetNbinsY();
684 double theta_max = hist->GetYaxis()->GetXmax();
685 double theta_min = hist->GetYaxis()->GetXmin();
686 trigcache[v] =
new TrigCache(nphi, (phi_max-phi_min)/nphi, phi_min, ntheta, (theta_max - theta_min)/ntheta,theta_min, UCorrelator::AntennaPositions::instance(v), use_bin_center);
693 std::vector<std::pair<int,int> > pairs;
694 pairs.reserve(NANTENNAS *NANTENNAS/2);
696 for (
int ant1 = 0; ant1 < NANTENNAS; ant1++)
698 if (disallowed_antennas & (1ul << ant1))
continue;
700 for (
int ant2 = ant1+1; ant2 < NANTENNAS; ant2++)
702 if (disallowed_antennas & (1ul << ant2))
continue;
704 pairs.push_back(std::pair<int,int>(ant1,ant2));;
709 unsigned nit = pairs.size();
711 for (
int i= 1; i < nthreads();i++)
717 #ifdef UCORRELATOR_OPENMP 718 #pragma omp parallel for 720 for (
unsigned it = 0; it < nit; it++)
722 doAntennas(pairs[it].first, pairs[it].second, &hists[0], &norms[0], trigcache[v]);
730 combineHists<TH2D,double>(nthreads(),&hists[0]);
732 combineHists<TH2D,double>(nthreads(),&norms[0]);
738 for (
int i = 0; i < (hist->GetNbinsX()+2) * (hist->GetNbinsY()+2); i++)
740 double val = hist->GetArray()[i];
741 if (val == 0)
continue;
742 int this_norm = norm->GetArray()[i];
743 hist->GetArray()[i] = this_norm > 2 ? val/this_norm : 0;
749 hist->SetEntries(nonzero);
750 norm->SetEntries(nonzero);
756 UCorrelator::Correlator::~Correlator()
759 for (
int i = 0; i < NUM_ANITAS+1; i++)
771 #ifdef UCORRELATOR_OPENMP 773 for (
int i = 1; i < nthreads(); i++)
787 TFile f(fname,
"RECREATE");
789 TTree * positions =
new TTree(
"positions",
"Positions");
792 TTree * tree =
new TTree(
"delays",
"Delays");
794 int ant1, ant2, ipol;
795 double phi, theta, delta_t, group_delay;
797 tree->Branch(
"pol",&ipol);
798 tree->Branch(
"ant1",&ant1);
799 tree->Branch(
"ant2",&ant2);
800 tree->Branch(
"phi",&phi);
801 tree->Branch(
"theta",&theta);
802 tree->Branch(
"delta_t",&delta_t);
803 tree->Branch(
"group_delay",&group_delay);
805 double ant_phi, ant_r, ant_z;
806 positions->Branch(
"ant",&ant1);
807 positions->Branch(
"phi",&ant_phi);
808 positions->Branch(
"r",&ant_r);
809 positions->Branch(
"z",&ant_z);
810 positions->Branch(
"pol",&ipol);
814 for (ipol = 0; ipol < 2; ipol++)
816 for (ant1= 0; ant1 < NANTENNAS; ant1++)
818 ant_phi = ap->
phiAnt[ipol][ant1];
819 ant_r = ap->
rAnt[ipol][ant1];
820 ant_z = ap->
zAnt[ipol][ant1];
823 for (ant2 = ant1+1; ant2 < NANTENNAS; ant2++)
825 for (phi = 0; phi <=360; phi += 2)
827 for (theta = -90; theta <=90; theta +=2)
837 group_delay = delay1-delay2;
void dumpDeltaTs(const char *file) const
TH2D * computeZoomed(double phi, double theta, int nphi, double dphi, int ntheta, double dtheta, int nant=0, TH2D *useme=0)
double getAntennaGroupDelay(double phidiff, double theta)
Correlator(int nphi, double phimin, double phimax, int ntheta, double theta_lowest, double theta_highest, bool use_bin_center=false, bool scale_by_cos_theta=false, double baseline_weight=0, double gain_sigma=0)
double getDeltaT(int ant1, int ant2, double phi, double theta, AnitaPol::AnitaPol_t pol, bool includeGroupDelay=false)
double getDeltaTFast(int ant1, int ant2, int phibin, int thetabin, AnitaPol::AnitaPol_t pol, const TrigCache *cache, bool includeGroupDelay=false)
double zAnt[2][NUM_SEAVEYS]
int getClosestAntennas(double phi, int N, int *closest, ULong64_t disallowed=0, AnitaPol::AnitaPol_t pol=AnitaPol::kHorizontal) const
double rAnt[2][NUM_SEAVEYS]
double phiAnt[2][NUM_SEAVEYS]
This class is intended to store all the necessary data about an ANITA event for filtering and analysi...
void compute(const FilteredAnitaEvent *event, AnitaPol::AnitaPol_t pol)
enum AnitaPol::EAnitaPol AnitaPol_t
Polarisation enumeration.