Correlator.cc
1 #include "FilteredAnitaEvent.h"
2 #include "TString.h"
3 #include "AntennaPositions.h"
4 #include "DeltaT.h"
5 #include "TTree.h"
6 #include "TFile.h"
7 #include "TrigCache.h"
8 #include "FFTtools.h"
9 #include <assert.h>
10 #include "AnitaGeomTool.h"
11 #include "Correlator.h"
12 #include "AnalysisWaveform.h"
13 #include "TStopwatch.h"
14 
15 #ifdef UCORRELATOR_OPENMP
16 #include <omp.h>
17 
18 #define SECTIONS _Pragma("omp parallel sections")
19 #define SECTION _Pragma("omp section")
20 
21 #else
22 
23 #define SECTIONS if(true)
24 #define SECTION if(true)
25 
26 
27 #endif
28 
29 namespace UCorrelator
30 {
32  {
33  public:
34 
35 #ifdef UCORRELATOR_OPENMP
36  omp_lock_t waveform_locks[NANTENNAS];
37  omp_lock_t correlation_locks[NANTENNAS][NANTENNAS];
38 
39 
40 
42  {
43  for (int i = 0; i < NANTENNAS; i++)
44  {
45  omp_init_lock(&waveform_locks[i]);
46 
47  for (int j = 0; j < NANTENNAS; j++)
48  {
49  omp_init_lock(&correlation_locks[i][j]);
50  }
51  }
52 
53 
54 
55  }
56  ~CorrelatorLocks()
57  {
58  for (int i = 0; i < NANTENNAS; i++)
59  {
60  omp_destroy_lock(&waveform_locks[i]);
61 
62  for (int j = 0; j < NANTENNAS; j++)
63  {
64  omp_destroy_lock(&correlation_locks[i][j]);
65  }
66  }
67 
68  }
69 #endif
70  };
71 }
72 
73 static int nthreads()
74 {
75 #ifdef UCORRELATOR_OPENMP
76  return omp_get_max_threads();
77 #else
78  return 1;
79 #endif
80 }
81 
82 
83 
84 
85 static int get_thread_id()
86 {
87 #ifdef UCORRELATOR_OPENMP
88  return omp_get_thread_num();
89 #else
90  return 0;
91 #endif
92 }
93 
94 
95 
96 
97 static int count_the_correlators = 1;
98 
99 static int count_the_zoomed_correlators = 1;
100 
101 
102 
103 /*Add all histograms to the first one ... quickly... assuming they're the same size and type */
104 template<class T,class A>
105 static void combineHists(int N, T ** hists )
106 {
107  if (N <2) return;
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++ )
111  {
112  __restrict A* v = hists[i]->GetArray();
113 #pragma omp simd
114  for(int j = 0; j < nbins; j++)
115  {
116  sum[j]+=v[j];
117  }
118  }
119 }
120 
121 
122 
123 #ifndef RAD2DEG
124 #define RAD2DEG (180/M_PI)
125 #endif
126 
127 #define PHI_SECTOR_ANGLE (360. / NUM_PHI)
128 
129 
130 
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)
133 {
134  TString histname = TString::Format("ucorr_corr_%d",count_the_correlators);
135  TString normname = TString::Format("ucorr_norm_%d",count_the_correlators++);
136 
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);
141 
142  hist->GetXaxis()->SetTitle("#phi");
143  hist->GetYaxis()->SetTitle("-#theta");
144 
145  hists.resize(nthreads());
146  norms.resize(nthreads());
147  hists[0] = hist;
148  norms[0] = norm;
149 
150 // printf(":: %p %p\n",hists[0],norms[0]);
151  for (int i= 1; i < nthreads();i++)
152  {
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);
157 // printf(":: %p %p\n",hists[i],norms[i]);
158  }
159 
160 
161 
162  use_bin_center = use_center;
163  memset(trigcache,0, sizeof(trigcache));
164 
165 
166  disallowed_antennas = 0;
167  pad_factor = 3;
168  max_phi = 75;
169  max_phi2 = max_phi*max_phi;
170 
171  memset(padded_waveforms, 0, NANTENNAS * sizeof(AnalysisWaveform*));
172  memset(correlations, 0, NANTENNAS * NANTENNAS * sizeof(AnalysisWaveform*));
173 
174 #ifdef UCORRELATOR_OPENMP
175  locks = new CorrelatorLocks;
176 #endif
177  ev = 0;
178  groupDelayFlag = 1;
179 }
180 
181 static int allowedPhisPairOfAntennas(double &lowerAngle, double &higherAngle, double &centerTheta1, double &centerTheta2, double &centerPhi1, double &centerPhi2, int ant1, int ant2, double max_phi, AnitaPol::AnitaPol_t pol)
182 {
183 
184  int phi1=AnitaGeomTool::Instance()->getPhiFromAnt(ant1);
185  int phi2=AnitaGeomTool::Instance()->getPhiFromAnt(ant2);
186  int allowedFlag=0;
187 
188  int upperlimit=phi2+2;//2 phi sectors on either side
189  int lowerlimit=phi2-2;
190 
191  if(upperlimit>NUM_PHI-1)upperlimit-=NUM_PHI;
192  if(lowerlimit<0)lowerlimit+=NUM_PHI;
193 
194  if (upperlimit>lowerlimit){
195  if (phi1<=upperlimit && phi1>=lowerlimit){//within 2 phi sectors of eachother
196  allowedFlag=1;
197  }
198  }
199  if (upperlimit<lowerlimit){
200  if (phi1<=upperlimit || phi1>=lowerlimit){
201  allowedFlag=1;
202 
203  }
204  }
205 
206  double centerAngle1, centerAngle2;
207 
208  const UCorrelator::AntennaPositions * ap = UCorrelator::AntennaPositions::instance();
209 
210  if (allowedFlag==1)
211  {
212  centerAngle1=ap->phiAnt[pol][ant1];
213  centerAngle2=ap->phiAnt[pol][ant2];
214 // assert(centerAngle1 == ap->phiAnt[0][ant1]);
215 
216  if (centerAngle2>centerAngle1)
217  {
218  lowerAngle=centerAngle2-max_phi;
219  higherAngle=centerAngle1+max_phi;
220  }
221  else
222  {
223  lowerAngle=centerAngle1-max_phi;
224  higherAngle=centerAngle2+max_phi;
225  }
226 
227  if (lowerAngle<0) lowerAngle+=360;
228  if (higherAngle>360) higherAngle-=360;
229 
230  }
231  else
232  {
233 
234  centerAngle1= 0;
235  centerAngle2= 0;
236  }
237  centerTheta1=10;//degrees down
238  centerTheta2=10;//degrees down
239  centerPhi1=centerAngle1;
240  centerPhi2=centerAngle2;
241 
242  return allowedFlag;
243 
244 }
245 
246 void UCorrelator::Correlator::reset()
247 {
248 for (int i = 0; i < NANTENNAS; i++)
249  {
250  if (padded_waveforms[i])
251  {
252  delete padded_waveforms[i];
253  padded_waveforms[i] = 0;
254  }
255 
256  for (int j = 0; j < NANTENNAS; j++)
257  {
258  if (correlations[i][j])
259  {
260  delete correlations[i][j];
261  correlations[i][j] = 0;
262  }
263  }
264  }
265 
266 
267 }
268 
269 
270 
271 
272 AnalysisWaveform * UCorrelator::Correlator::getCorrelation(int ant1, int ant2)
273 {
274 // printf("%d %d / %d \n",ant1,ant2, NANTENNAS);
275 
276 
277 #ifdef UCORRELATOR_OPENMP
278 #ifndef FFTTOOLS_COMPILED_WITH_OPENMP
279 #pragma omp critical (getCorrelation)
280  {
281 #endif
282 #endif
283 
284 
285 #ifdef UCORRELATOR_OPENMP
286  omp_set_lock(&locks->waveform_locks[ant1]);
287 #endif
288  if (!padded_waveforms[ant1])
289  {
290 // printf("Copying and padding %d\n",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);
294  }
295 #ifdef UCORRELATOR_OPENMP
296  omp_unset_lock(&locks->waveform_locks[ant1]);
297  omp_set_lock(&locks->waveform_locks[ant2]);
298 #endif
299 
300  if (!padded_waveforms[ant2])
301  {
302 // printf("Copying and padding %d\n",ant2);
303  padded_waveforms[ant2] = new AnalysisWaveform(*ev->getFilteredGraph(ant2, pol));
304 // printf("Computing rms!\n");
305  rms[ant2] = padded_waveforms[ant2]->even()->GetRMS(2);
306  padded_waveforms[ant2]->padEven(1);
307  }
308 
309 #ifdef UCORRELATOR_OPENMP
310  omp_unset_lock(&locks->waveform_locks[ant2]);
311  omp_set_lock(&locks->correlation_locks[ant1][ant2]);
312 #endif
313 
314  if (!correlations[ant1][ant2])
315  {
316 #ifdef UCORRELATOR_OPENMP
317  omp_set_lock(&locks->waveform_locks[ant1]);
318  omp_set_lock(&locks->waveform_locks[ant2]);
319 #endif
320 // printf("Computing correlation %d %d\n", ant1, ant2);
321  correlations[ant1][ant2] = AnalysisWaveform::correlation(padded_waveforms[ant1],padded_waveforms[ant2],pad_factor, rms[ant1] * rms[ant2]);
322 
323 #ifdef UCORRELATOR_OPENMP
324  omp_unset_lock(&locks->waveform_locks[ant2]);
325  omp_unset_lock(&locks->waveform_locks[ant1]);
326 #endif
327  }
328 
329 
330 #ifdef UCORRELATOR_OPENMP
331  omp_unset_lock(&locks->correlation_locks[ant1][ant2]);
332 #ifndef FFTTOOLS_COMPILED_WITH_OPENMP
333  }
334 #endif
335 #endif
336  return correlations[ant1][ant2];
337 }
338 
339 
340 
341 TH2D * UCorrelator::Correlator::computeZoomed(double phi, double theta, int nphi, double dphi, int ntheta, double dtheta, int nant, TH2D * answer)
342 {
343 
344  if (!ev)
345  {
346  fprintf(stderr, "Must call Correlator::compute() prior to Correlator::computeZoomed!!!!");
347  return 0;
348  }
349 
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",
355  nphi, phi0,phi1,
356  ntheta, theta0, theta1);
357  zoomed_norm->SetDirectory(0);
358 
359  if (answer)
360  {
361  answer->SetBins(nphi, phi0, phi1,
362  ntheta, theta0, theta1);
363  answer->Reset();
364  }
365  else
366  {
367  answer = new TH2D(TString::Format("zoomed_corr_%d", count_the_zoomed_correlators++), "Zoomed Correlation",
368  nphi, phi0, phi1,
369  ntheta, theta0, theta1);
370  answer->SetDirectory(0);
371  }
372 
373  answer->GetXaxis()->SetTitle("#phi");
374  answer->GetYaxis()->SetTitle("-#theta");
375 
376  const AntennaPositions * ap = AntennaPositions::instance();
377 
378  int closest[nant]; // is it a problem if nant is 0?
379  if (nant)
380  {
381  memset(closest,0,sizeof(closest));
382  nant = ap->getClosestAntennas(phi, nant, closest, disallowed_antennas);
383  }
384 
385  TrigCache cache(nphi, dphi, phi0, ntheta,dtheta,theta0, ap, true,nant, nant ? closest : 0);
386 
387  int n2loop = nant ? nant : NANTENNAS;
388 
389  std::vector<std::pair<int,int> > pairs;
390  pairs.reserve(n2loop);
391 
392  double center_point[2];
393  center_point[0] = phi;
394  center_point[1] = theta;
395 
396  for (int ant_i = 0; ant_i < n2loop; ant_i++)
397  {
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++)
401  {
402  int ant2 = nant ? closest[ant_j] : ant_j;
403  if (!nant && disallowed_antennas & (1ul << ant2)) continue;
404 
405  pairs.push_back(std::pair<int,int>(ant1,ant2));
406  }
407  }
408 
409  unsigned nit = pairs.size();
410 
411 
412 
413  /* lock contention for the hist / norm lock is killing parallelization */
414  std::vector<TH2D*> zoomed_hists(nthreads());
415  std::vector<TH2D*> zoomed_norms(nthreads());
416 
417  zoomed_hists[0] = answer;
418  zoomed_norms[0] = zoomed_norm;
419 
420 // printf(":: %p %p\n",hists[0],norms[0]);
421  for (int i= 1; i < nthreads();i++)
422  {
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);
427 // zoomed_ printf(":: %p %p\n",hists[i],norms[i]);
428  }
429 
430 
431 #ifdef UCORRELATOR_OPENMP
432 #pragma omp parallel for
433 #endif
434  for (unsigned it = 0; it < nit; it++)
435  {
436  doAntennas(pairs[it].first, pairs[it].second, &zoomed_hists[0], &zoomed_norms[0], &cache, center_point);
437  }
438 
439 
440  /* combine histograms */
441 
442 SECTIONS
443  {
444 SECTION
445  combineHists<TH2D,double>(nthreads(), &zoomed_hists[0]);
446 
447 SECTION
448  combineHists<TH2D,double>(nthreads(), &zoomed_norms[0]);
449  }
450 
451  for (int i =1; i < nthreads(); i++)
452  {
453  delete zoomed_hists[i];
454  delete zoomed_norms[i];
455  }
456 
457 
458  int nonzero = 0;
459  //only keep values with contributing antennas
460  for (int i = 0; i < (answer->GetNbinsX()+2) * (answer->GetNbinsY()+2); i++)
461  {
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;
466  nonzero++;
467  }
468 
469  answer->SetEntries(nonzero);
470 
471  delete zoomed_norm;
472  return answer;
473 }
474 
475 
476 static inline bool between(double phi, double low, double high)
477 {
478 
479  double diff_highlow = fmod(high - low + 360, 360);
480  double diff_philow = fmod(phi - low + 360, 360);
481 
482  return diff_philow < diff_highlow;
483 }
484 
485 
486 inline void UCorrelator::Correlator::doAntennas(int ant1, int ant2, TH2D ** these_hists,
487  TH2D ** these_norms, const UCorrelator::TrigCache * cache ,
488  const double * center_point )
489 {
490  int allowedFlag;
491  double lowerAngleThis, higherAngleThis, centerTheta1, centerTheta2, centerPhi1, centerPhi2;
492 
493  allowedFlag=allowedPhisPairOfAntennas(lowerAngleThis,higherAngleThis,
494  centerTheta1, centerTheta2, centerPhi1, centerPhi2,
495  ant1,ant2, max_phi, pol);
496 
497 
498  assert(ant1 < 48);
499  assert(ant2 < 48);
500  if(!allowedFlag) return;
501 
502  TH2D * the_hist = these_hists[get_thread_id()];
503  TH2D * the_norm = these_norms[get_thread_id()];
504 
505 // printf("lowerAngleThis: %g higherAngleThis: %g\n", lowerAngleThis, higherAngleThis);
506  // More stringent check if we have a center point
507  if (center_point && !between(center_point[0], lowerAngleThis, higherAngleThis)) return;
508 
509  AnalysisWaveform * correlation = getCorrelation(ant1,ant2);
510 
511  int nphibins = the_hist->GetNbinsX() + 2;
512 
513  //find phi bin corresponding to lowerAngleThis and higherAngleThis
514 
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);
517 
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) ;
521 
522 
523  //So the maximum number of bins is going to be the total number of bins in the histogram. We probably won't fill all of them,
524  //but memory is cheap and std::vector is slow
525 
526  int maxsize = the_hist->GetNbinsY() * the_hist->GetNbinsX();
527 
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;
531  baseline_phi = FFTtools::wrap(baseline_phi, 360);
532 
533  double baseline_theta = (centerTheta1 + centerTheta2) / 2;
534 
535 
536  //This is bikeshedding, but allocate it all contiguosly
537  int * alloc = new int[3*maxsize];
538  int * phibins = alloc;
539  int * thetabins = alloc + maxsize;
540  int * bins_to_fill = alloc + 2 *maxsize;
541 
542  int nbins_used =0;
543 
544  double * gain_weights = 0;
545  if (gainSigma && !center_point) gain_weights = new double[maxsize];
546 
547 
548  for (int phibin = first_phi_bin; (phibin <= last_phi_bin) || must_wrap; phibin++)
549  {
550  if (must_wrap && phibin == nphibins-1)
551  {
552  phibin = 1;
553  must_wrap = false;
554  }
555 
556  double phi = cache->phi[phibin-1] ;
557 // double phi4width = center_point ? center_point[0] : phi;
558  double dphi1 = center_point ? 0 : FFTtools::wrap(phi - centerPhi1,360,0);
559  double dphi2 = center_point ? 0 : FFTtools::wrap(phi - centerPhi2,360,0);
560 
561 
562  //Check if in beam width in phi
563  if (!center_point && fabs(dphi1) > max_phi) continue;
564  if (!center_point && fabs(dphi2) > max_phi) continue;
565 
566  int ny = the_hist->GetNbinsY();
567 
568 
569  for (int thetabin = 1; thetabin <= ny; thetabin++)
570  {
571  double theta = cache->theta[thetabin-1];
572 // double theta4width = center_point ? center_point[1] : theta;
573  double dtheta1 = center_point ? 0 : FFTtools::wrap(theta - centerTheta1, 360, 0);
574  double dtheta2 = center_point ? 0 : FFTtools::wrap(theta - centerTheta2, 360, 0);
575 
576  // check if in beam width
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;
579 
580  phibins[nbins_used] = phibin;
581  thetabins[nbins_used] = thetabin;
582  bins_to_fill[nbins_used] = phibin + thetabin * nphibins;
583  if (gainSigma && !center_point)
584  {
585  //Matt Mottram weighted by the baseline angle difference somehow
586 
587  double Dphi = FFTtools::wrap(phi - baseline_phi, 360, 0);
588  double Dtheta = theta - baseline_theta;
589 
590  gain_weights[nbins_used] = TMath::Gaus(TMath::Sqrt(Dphi * Dphi + Dtheta * Dtheta), 0, gainSigma);
591  }
592  nbins_used++;
593  }
594  }
595 
596 
597  double * dalloc = new double[2 *nbins_used];
598  double * __restrict__ vals_to_fill = dalloc;
599  double * __restrict__ times_to_fill = dalloc + nbins_used;
600 
601  //TODO vectorize this
602  for (int i = 0; i < nbins_used; i++)
603  {
604 
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);
608  }
609 
610  correlation->evalEven(nbins_used, times_to_fill, vals_to_fill);
611 
612 
613  if (scale_cos_theta)
614  {
615  for(int i = 0; i < nbins_used; i++)
616  {
617  vals_to_fill[i] *= cache->cos_theta[thetabins[i]-1];
618  }
619  }
620 
621 
622  if (baselineWeight)
623  {
624  double wgt = TMath::Power(cache->ap->distance(ant1, ant2, pol), baselineWeight);
625  for (int i = 0; i < nbins_used; i++)
626  {
627  vals_to_fill[i] *= wgt;
628  }
629  }
630 
631 
632 
633  double * __restrict__ the_arr = the_hist->GetArray();
634  double * __restrict__ the_norm_arr = the_norm->GetArray();
635 
636  for (int bi = 0; bi < nbins_used; bi++)
637  {
638  double val = vals_to_fill[bi];
639  int bin = bins_to_fill[bi];
640 
641  if (gainSigma && !center_point)
642  {
643  the_arr[bin] += val * gain_weights[bi];
644  the_norm_arr[bin] += gain_weights[bi];
645  }
646  else
647  {
648  the_arr[bin]+= val;
649  the_norm_arr[bin]+= 1;
650  }
651  }
652 
653 
654  delete [] alloc;
655  delete [] dalloc;
656  if (gain_weights) delete [] gain_weights;
657 }
658 
659 
661 {
662 
663 // TStopwatch sw;
664 
665  pol = whichpol;
666  ev = event;
667  hist->Reset();
668  norm->Reset();
669  reset();
670 
671  //alright, we have to be able to dispatch between different versions
672  //because who knows what crazy things we might be asked to do.
673  // I suppose we could just require the user to make a different Correlator for each ANITA version... but whatever
674 
675  int v = event->getAnitaVersion();
676  AnitaVersion::set(v);
677 
678  if (! trigcache[v])
679  {
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);
687  }
688 
689 
690  //precompute antenna combinations
691 
692 
693  std::vector<std::pair<int,int> > pairs;
694  pairs.reserve(NANTENNAS *NANTENNAS/2);
695 
696  for (int ant1 = 0; ant1 < NANTENNAS; ant1++)
697  {
698  if (disallowed_antennas & (1ul << ant1)) continue;
699 
700  for (int ant2 = ant1+1; ant2 < NANTENNAS; ant2++)
701  {
702  if (disallowed_antennas & (1ul << ant2)) continue;
703 
704  pairs.push_back(std::pair<int,int>(ant1,ant2));;
705  }
706 
707  }
708 
709  unsigned nit = pairs.size();
710 
711  for (int i= 1; i < nthreads();i++)
712  {
713  hists[i]->Reset();
714  norms[i]->Reset();
715  }
716 
717 #ifdef UCORRELATOR_OPENMP
718  #pragma omp parallel for
719 #endif
720  for (unsigned it = 0; it < nit; it++)
721  {
722  doAntennas(pairs[it].first, pairs[it].second, &hists[0], &norms[0], trigcache[v]);
723  }
724 
725  /* combine histograms */
726 SECTIONS
727 {
728 
729  SECTION
730  combineHists<TH2D,double>(nthreads(),&hists[0]);
731  SECTION
732  combineHists<TH2D,double>(nthreads(),&norms[0]);
733 }
734 
735 
736  int nonzero = 0;
737  //only keep values with at least 3 contributing antennas
738  for (int i = 0; i < (hist->GetNbinsX()+2) * (hist->GetNbinsY()+2); i++)
739  {
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;
744 
745  // printf("%d %g %d\n", i, val, this_norm);
746  nonzero++;
747  }
748 
749  hist->SetEntries(nonzero);
750  norm->SetEntries(nonzero);
751 // sw.Print("u");
752 }
753 
754 
755 
756 UCorrelator::Correlator::~Correlator()
757 {
758 
759  for (int i = 0; i < NUM_ANITAS+1; i++)
760  {
761  if (trigcache[i])
762  {
763  delete trigcache[i];
764  }
765  }
766 
767  reset();
768  delete hist;
769  delete norm;
770 
771 #ifdef UCORRELATOR_OPENMP
772  delete locks;
773  for (int i = 1; i < nthreads(); i++)
774  {
775  delete hists[i];
776  delete norms[i];
777 
778  }
779 #endif
780 
781 }
782 
783 
784 void UCorrelator::Correlator::dumpDeltaTs(const char * fname) const
785 {
786 
787  TFile f(fname,"RECREATE");
788 
789  TTree * positions = new TTree("positions","Positions");
790 
791 
792  TTree * tree = new TTree("delays","Delays");
793 
794  int ant1, ant2, ipol;
795  double phi, theta, delta_t, group_delay;
796 
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);
804 
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);
811 
812  const AntennaPositions * ap = AntennaPositions::instance();
813 
814  for (ipol = 0; ipol < 2; ipol++)
815  {
816  for (ant1= 0; ant1 < NANTENNAS; ant1++)
817  {
818  ant_phi = ap->phiAnt[ipol][ant1];
819  ant_r = ap->rAnt[ipol][ant1];
820  ant_z = ap->zAnt[ipol][ant1];
821  positions->Fill();
822 
823  for (ant2 = ant1+1; ant2 < NANTENNAS; ant2++)
824  {
825  for (phi = 0; phi <=360; phi += 2)
826  {
827  for (theta = -90; theta <=90; theta +=2)
828  {
829  delta_t = getDeltaT(ant1,ant2,phi,theta,(AnitaPol::AnitaPol_t) ipol, groupDelayFlag);
830 
831  if (groupDelayFlag)
832  {
833  double dphi1 = FFTtools::wrap(phi - ap->phiAnt[ipol][ant1],360,0);
834  double dphi2 = FFTtools::wrap(phi - ap->phiAnt[ipol][ant2],360,0);
835  double delay1=getAntennaGroupDelay(dphi1,theta);
836  double delay2=getAntennaGroupDelay(dphi2,theta);
837  group_delay = delay1-delay2;
838  }
839  f.cd();
840  tree->Fill();
841  }
842  }
843  }
844  }
845  }
846 
847  f.cd();
848  positions->Write();
849  tree->Write();
850  f.Close();
851 
852 
853 }
854 
855 
double evalEven(double t, EvenEvaluationType=EVAL_LINEAR) const
void dumpDeltaTs(const char *file) const
Definition: Correlator.cc:784
TH2D * computeZoomed(double phi, double theta, int nphi, double dphi, int ntheta, double dtheta, int nant=0, TH2D *useme=0)
Definition: Correlator.cc:341
double getAntennaGroupDelay(double phidiff, double theta)
Definition: DeltaT.h:17
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)
Definition: Correlator.cc:131
double getDeltaT(int ant1, int ant2, double phi, double theta, AnitaPol::AnitaPol_t pol, bool includeGroupDelay=false)
Definition: DeltaT.h:68
double getDeltaTFast(int ant1, int ant2, int phibin, int thetabin, AnitaPol::AnitaPol_t pol, const TrigCache *cache, bool includeGroupDelay=false)
Definition: DeltaT.h:96
double zAnt[2][NUM_SEAVEYS]
int getClosestAntennas(double phi, int N, int *closest, ULong64_t disallowed=0, AnitaPol::AnitaPol_t pol=AnitaPol::kHorizontal) const
static Int_t getPhiFromAnt(Int_t ant)
get phi from ant
double rAnt[2][NUM_SEAVEYS]
void wrap(size_t N, double *vals, double period=360)
Definition: FFTtools.cxx:2851
double phiAnt[2][NUM_SEAVEYS]
static AnalysisWaveform * correlation(const AnalysisWaveform *A, const AnalysisWaveform *B, int npadfreq=0, double scale=1, int window_normalize=0)
const TGraphAligned * even() const
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)
Definition: Correlator.cc:660
enum AnitaPol::EAnitaPol AnitaPol_t
Polarisation enumeration.
static AnitaGeomTool * Instance(int anita_version=0)
Instance generator. If version_number == 0, uses AnitaVersion::get();.
This class is intended to be the main storage vessel for ANITA waveforms. It is similar in principle ...