AnitaGeomTool.cxx
1 
9 #include "AnitaGeomTool.h"
10 #include "AnitaVersion.h"
11 #include <assert.h>
12 #include "TMutex.h"
13 
14 #define INCHTOMETER 0.0254
15 
16 Double_t deltaRL = 0.0;
17 Double_t deltaUD = 0.0;
18 
19 
20 static AnitaGeomTool * instances[NUM_ANITAS+1] = {0,0,0,0,0};
21 
22 
23 
25 
29 namespace AnitaGeom {
30 
32  Int_t antToSurfMap[NUM_SEAVEYS]={11,5,10,4,11,4,10,5,11,5,10,4,11,4,10,5,
33  9,3,8,2,8,3,9,2,9,3,8,2,8,3,9,2,
34  6,0,7,1,6,1,7,0,6,0,7,1,6,1,7,0};
35 
37  Int_t vAntToChanAnita3[NUM_SEAVEYS]={3,1,3,5,1,3,1,3,2,0,2,0,0,2,0,2,
38  1,3,1,3,3,1,3,1,0,2,0,2,2,0,2,0,
39  3,1,3,1,1,3,1,3,2,0,2,0,0,2,0,2};
40  Int_t vAntToChan[NUM_SEAVEYS]={3,1,3,5,1,3,1,3,2,0,2,0,0,2,0,2,
41  1,3,1,3,3,1,3,1,0,2,0,2,2,0,2,0,
42  3,1,3,1,1,3,1,3,2,0,2,0,0,2,0,2};
43 
45  Int_t hAntToChanAnita3[NUM_SEAVEYS]={7,5,7,1,5,7,5,7,6,4,6,4,4,6,4,6,
46  5,7,5,7,7,5,7,5,4,6,4,6,6,4,6,4,
47  7,5,7,5,5,7,5,7,6,4,6,4,4,6,4,6};
48  Int_t hAntToChan[NUM_SEAVEYS]={7,5,7,1,5,7,5,7,6,4,6,4,4,6,4,6,
49  5,7,5,7,7,5,7,5,4,6,4,6,6,4,6,4,
50  7,5,7,5,5,7,5,7,6,4,6,4,4,6,4,6};
51 
52 
54  // Apparently -2 is a 90 degree flip.
55  Int_t antOrientationMapAnita3[NUM_SEAVEYS] =
56  {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
57  1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
58  1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
59 
60  Int_t antOrientationMap[NUM_SEAVEYS] =
61  {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
62  1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
63  1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
64 
65 
69  Int_t surfToAntMapAnita4[ACTIVE_SURFS][RFCHAN_PER_SURF]= {{-42,-34,-48,-40,42,34,48,40},
70  {-44,-36,-46,-38,44,36,46,38},
71  {-32,-24,-28,-20,32,24,28,20},
72  {-30,-22,-26,-18,30,22,26,18},
73  {-12, 4,-14,-6,12,-4,14,6}, //ABL changed -/+4
74  {-10,-2,-16,-8,10,2,16,8},
75  {-45,-37,-41,-33,45,37,41,33},
76  {-47,-39,-43,-35,47,39,43,35},
77  {-27,-19,-29,-21,27,19,29,21},
78  {-25,-17,-31,-23,25,17,31,23},
79  {-15,-7,-11,-3,15,7,11,3},
80  {-13,-5,-9,-1,13,5,9,1}};
81 
82  Int_t surfToAntMapAnita3[ACTIVE_SURFS][RFCHAN_PER_SURF]= {{-42,-34,-48,-40,42,34,48,40},
83  {-44,-36,-46,-38,44,36,46,38},
84  {-32,-24,-28,-20,32,24,28,20},
85  {-30,-22,-26,-18,30,22,26,18},
86  {-12,4,-14,-6,12,-4,14,6},
87  {-10,-2,-16,-8,10,2,16,8},
88  {-45,-37,-41,-33,45,37,41,33},
89  {-47,-39,-43,-35,47,39,43,35},
90  {-27,-19,-29,-21,27,19,29,21},
91  {-25,-17,-31,-23,25,17,31,23},
92  {-15,-7,-11,-3,15,7,11,3},
93  {-13,-5,-9,-1,13,5,9,1}};
94 
96  Int_t topAntNums[NUM_PHI] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
97  Int_t middleAntNums[NUM_PHI] = {16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31};
98  Int_t bottomAntNums[NUM_PHI] = {32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47};
99 
101  Int_t topPhiNums[NUM_PHI] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
102  Int_t middlePhiNums[NUM_PHI] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
103  Int_t bottomPhiNums[NUM_PHI] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15};
104 
105 
107  Int_t surfToPhiTriggerMapAnita3[ACTIVE_SURFS][2]={{-1,-1}, {-1,-1}, { 0, 4}, { 2, 6},
108  { 1, 5}, { 3, 7}, {15,11}, {13, 9},
109  {14,10}, {12, 8}, {-1,-1}, {-1,-1}};
110  Int_t surfToPhiTriggerMapAnita4[ACTIVE_SURFS][2]={{-1,-1}, {-1,-1}, { 0, 1}, { 8, 9},
111  { 2, 3}, { 10,11}, {4,5}, {12, 13},
112  {6,7}, {14, 15}, {-1,-1}, {-1,-1}};
113 
114 
115  using AnitaRing::kTopRing;
118 
120  AnitaRing::AnitaRing_t surfTriggerChanToRing[SCALERS_PER_SURF]={kTopRing, kMiddleRing, kBottomRing,
121  kTopRing, kMiddleRing, kBottomRing,
122  kTopRing, kMiddleRing, kBottomRing,
123  kTopRing, kMiddleRing, kBottomRing};
124 
125 
129  AnitaTrigPol::AnitaTrigPol_t surfTriggerChanToPolAnita3[SCALERS_PER_SURF]={kVertical, kVertical, kVertical,
130  kVertical, kVertical, kVertical,
131  kHorizontal, kHorizontal, kHorizontal,
132  kHorizontal, kHorizontal, kHorizontal};
133 
134  using AnitaTrigPol::kRCP;
135  using AnitaTrigPol::kLCP;
136  AnitaTrigPol::AnitaTrigPol_t surfTriggerChanToPolAnita4[SCALERS_PER_SURF]={kRCP,kRCP,kRCP,
137  kRCP,kRCP,kRCP,
138  kLCP,kLCP,kLCP,
139  kLCP,kLCP,kLCP};
140 
141  Int_t phiToSurfTriggerMapAnita4[NUM_PHI] = {2,2,4,4,6,6,8,8,3,3,5,5,7,7,9,9};
142  Int_t phiToSurfHalfAnita4[NUM_PHI] = {0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1};
143 
144 
145 
146  //These are from ANITA3
147  Int_t phiToSurfTriggerMapAnita3[NUM_PHI] = {2,4,3,5,2,4,3,5,9,7,8,6,9,7,8,6};
148  Int_t phiToSurfHalfAnita3[NUM_PHI] = {0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0};
149 
150 }
151 
152 
153 
157 AnitaGeomTool::AnitaGeomTool()
158 {
159  int v = AnitaVersion::get();
160 
161  //Default constructor
162  ringPhaseCentreOffset[0]=0.2-0.042685;
163  ringPhaseCentreOffset[1]=0.2+0.00653;
164  ringPhaseCentreOffset[2]=0.2+0.1927;
165 
166  // readSimonsNumbers();
167  // readAnita2Photogrammetry();
168 
169 
170  fUsePhotogrammetryNumbers=0;
171  //TODO need to add phase-centers for A4
172  readPhaseCenterNumbers(v); // This one has to come before reading photogrammetry
173  readPhotogrammetry(v);
174 
175  // Moved from fillAntPositionsFromPrioritizerdConfig
176  // Who knows where else these get set in this crazy place...
177  fHeadingRotationAxis.SetXYZ(0.,0.,1.);
178  fPitchRotationAxis.SetXYZ(0.,1.,0.);
179  fRollRotationAxis=fPitchRotationAxis.Cross(fHeadingRotationAxis);
180  aftForeOffsetAngleVertical=TMath::DegToRad()*45;
181  // for(int ant=0; ant<NUM_SEAVEYS; ant++){
182  // for(int pol=0; pol<2; pol++){
183  // std::cout << ant << "\t" << pol<< "\t"
184  // << rPhaseCentreFromVerticalHorn[ant][pol] << "\t"
185  // << rPhaseCentreFromVerticalHornPhotogrammetry[ant][pol] << "\t"
186  // << zPhaseCentreFromVerticalHorn[ant][pol] << "\t"
187  // << zPhaseCentreFromVerticalHornPhotogrammetry[ant][pol] << "\t"
188  // << azPhaseCentreFromVerticalHorn[ant][pol] << "\t"
189  // << azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol] << "\t"
190  // << std::endl;
191  // }
192  // }
193 
194  // fillAntPositionsFromPrioritizerdConfig();
195 
196 
197 }
198 
199 
200 
205 {
206  //Default destructor
207 }
208 
209 
210 
215 
216  Double_t phiArrayDeg[NUM_SEAVEYS]={0,22.5,45,67.5,90,112.5,135,157.5,180,
217  202.5,225,247.5,270,292.5,315,337.5,
218  0,22.5,45,67.5,90,112.5,135,157.5,180,
219  202.5,225,247.5,270,292.5,315,337.5,
220  0,22.5,45,67.5,90,112.5,135,157.5,180,
221  202.5,225,247.5,270,292.5,315,337.5};
222 
223  Double_t rArray[NUM_SEAVEYS]={0.9675,0.7402,0.9675,0.7402,0.9675,0.7402,0.9675,0.7402,
224  0.9675,0.7402,0.9675,0.7402,0.9675,0.7402,0.9675,0.7402,
225  2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,
226  2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,
227  2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,
228  2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447,2.0447};
229  Double_t zArray[NUM_SEAVEYS]={-1.4407,-2.4135,-1.4407,-2.4135,-1.4407,-2.4135,-1.4407,-2.4135,
230  -1.4407,-2.4135,-1.4407,-2.4135,-1.4407,-2.4135,-1.4407,-2.4135,
231  -5.1090,-5.1090,-5.1090,-5.1090,-5.1090,-5.1090,-5.1090,-5.1090,
232  -5.1090,-5.1090,-5.1090,-5.1090,-5.1090,-5.1090,-5.1090,-5.1090,
233  -6.1951,-6.1951,-6.1951,-6.1951,-6.1951,-6.1951,-6.1951,-6.1951,
234  -6.1951,-6.1951,-6.1951,-6.1951,-6.1951,-6.1951,-6.1951,-6.1951};
235 
236 
237  for(Int_t ant=0;ant<NUM_SEAVEYS;ant++) {
238  for(Int_t pol=0;pol<2;pol++) {
239  zPhaseCentreFromVerticalHorn[ant][pol]=zArray[ant];
240  rPhaseCentreFromVerticalHorn[ant][pol]=rArray[ant];
241  azPhaseCentreFromVerticalHorn[ant][pol]=phiArrayDeg[ant]*TMath::DegToRad();
242  xPhaseCentreFromVerticalHorn[ant][pol]=rArray[ant]*TMath::Cos(azPhaseCentreFromVerticalHorn[ant][pol]);
243  yPhaseCentreFromVerticalHorn[ant][pol]=rArray[ant]*TMath::Sin(azPhaseCentreFromVerticalHorn[ant][pol]);
244  }
245  }
246 
247  aftForeOffsetAngleVertical=TMath::DegToRad()*45;
248 
249  fHeadingRotationAxis.SetXYZ(0.,0.,1.);
250  fPitchRotationAxis.SetXYZ(0.,1.,0.);
251  fRollRotationAxis=fPitchRotationAxis.Cross(fHeadingRotationAxis);
252 
253 }
254 
255 
256 
257 
261 Int_t AnitaGeomTool::getPhiRingPolFromSurfChanTrigger(Int_t surf,Int_t chan, Int_t &phi,
263 {
264  if(surf<0 || surf>=ACTIVE_SURFS) return -1;
265  if(chan<0 || chan>=SCALERS_PER_SURF) return -1;
266 
267  Int_t surfHalf=0;
268  if((chan%6)>=3) surfHalf=1;
269 
270  Int_t v = AnitaVersion::get();
271 
272  phi = v == 4 ? AnitaGeom::surfToPhiTriggerMapAnita4[surf][surfHalf] :
273  v == 3 ? AnitaGeom::surfToPhiTriggerMapAnita3[surf][surfHalf] :
274  -1;
275 
277 
278  pol = v == 4 ? AnitaGeom::surfTriggerChanToPolAnita4[chan] :
281 
282 
283  return phi;
284 }
285 
289 Int_t AnitaGeomTool::getPhiRingFromSurfL1Chan(Int_t surf, Int_t l1Chan, Int_t &phi, AnitaRing::AnitaRing_t &ring)
290 {
291  if(phi<0 || phi>=NUM_PHI) return -1;
292  Int_t surfHalf=l1Chan>2;
293 
294  Int_t v = AnitaVersion::get();
295  // this method only makes sends for A4
296  assert(v==4);
297 
298  phi = AnitaGeom::surfToPhiTriggerMapAnita4[surf][surfHalf] ;
299 
300  if((l1Chan%3)==0) ring=AnitaRing::kTopRing;
301  if((l1Chan%3)==1) ring=AnitaRing::kMiddleRing;
302  if((l1Chan%3)==2) ring=AnitaRing::kBottomRing;
303  return phi;
304 }
305 
306 
310 Int_t AnitaGeomTool::getPhiPolFromSurfL1Chan(Int_t surf, Int_t l1Chan, Int_t &phi, AnitaPol::AnitaPol_t & pol)
311 {
312  if(phi<0 || phi>=NUM_PHI) return -1;
313 
314  Int_t v = AnitaVersion::get();
315 
316  //this method only makes sense for A3
317  assert(v == 3);
318 
319  Int_t surfHalf=l1Chan%2;
320  phi=AnitaGeom::surfToPhiTriggerMapAnita3[surf][surfHalf];
322  if(l1Chan>=2) pol=AnitaPol::kHorizontal;
323  return phi;
324 }
325 
330  if(phi<0 || phi>=NUM_PHI) return -1;
331 
332  Int_t v = AnitaVersion::get();
333 
334  surf = v == 4 ? AnitaGeom::phiToSurfTriggerMapAnita4[phi] :
335  v == 3 ? AnitaGeom::phiToSurfTriggerMapAnita3[phi] :
336  -1 ;
337 
338  Int_t surfHalf = v == 4 ? AnitaGeom::phiToSurfHalfAnita4[phi] :
339  v == 3 ? AnitaGeom::phiToSurfHalfAnita3[phi] :
340  -1;
341 
342  chan = (6-6*(pol % 2))+ ring + 3*surfHalf;
343 
344 
345  if(v == 3 && phi==3 && ring==AnitaRing::kTopRing) {
346  if(pol==AnitaTrigPol::kVertical) chan=6;
347  if(pol==AnitaTrigPol::kHorizontal) chan=0;
348  }
349  return surf;
350 }
351 
352 
356 Int_t AnitaGeomTool::getSurfL1TriggerChanFromPhiRing(Int_t phi, AnitaRing::AnitaRing_t ring, Int_t &surf, Int_t &l1Chan)
357 {
358  if(phi<0 || phi>=NUM_PHI) return -1;
359 
360  Int_t v = AnitaVersion::get();
361 
362 
363  surf = v == 4 ? AnitaGeom::phiToSurfTriggerMapAnita4[phi] :
364  v == 3 ? AnitaGeom::phiToSurfTriggerMapAnita3[phi] :
365  -1 ;
366 
367  Int_t surfHalf = v == 4 ? AnitaGeom::phiToSurfHalfAnita4[phi] :
368  v == 3 ? AnitaGeom::phiToSurfHalfAnita3[phi] :
369  -1;
370 
371  switch(ring) {
372  case AnitaRing::kTopRing:
373  l1Chan=3*surfHalf;
374  break;
376  l1Chan=1+(3*surfHalf);
377  break;
379  l1Chan=2+(3*surfHalf);
380  break;
381  default:
382  std::cerr << "Unknown ring " << ring << "\n";
383  }
384  return surf;
385 }
386 
387 Int_t AnitaGeomTool::getSurfL1TriggerChanFromPhiPol(Int_t phi, AnitaPol::AnitaPol_t pol, Int_t &surf, Int_t &l1Chan)
388 {
389  if(phi<0 || phi>=NUM_PHI) return -1;
390 
391  Int_t v = AnitaVersion::get();
392  assert(v==3);
393 
394  surf=AnitaGeom::phiToSurfTriggerMapAnita3[phi];
395  Int_t surfHalf=AnitaGeom::phiToSurfHalfAnita3[phi];
396  if(pol==AnitaPol::kVertical)
397  {
398  l1Chan=surfHalf;
399  }
400  else
401  {
402  l1Chan=2+surfHalf;
403  }
404 
405  return surf;
406 }
407 
411 Int_t AnitaGeomTool::getSurfL2TriggerChanFromPhi(Int_t phi, Int_t &surf, Int_t &l2Chan)
412 {
413  if(phi<0 || phi>=NUM_PHI) return -1;
414 
415  Int_t v = AnitaVersion::get();
416  assert(v==4);
417  surf = AnitaGeom::phiToSurfTriggerMapAnita4[phi] ;
418  Int_t surfHalf= AnitaGeom::phiToSurfHalfAnita4[phi] ;
419 
420  l2Chan=surfHalf;
421  return surf;
422 }
423 
424 
425 
426 
427 static TMutex instance_lock;
432 {
433 
434  /**************************************************************************************************
435  * *
436  * Yeah, the classic singleton pattern isn't thread-safe. *
437  * We got away it before because we usually statically initialized AnitaGeomTool instances, but *
438  * with multiple versions, we can't necessarily do that. *
439  * *
440  * If we required c++11, static initialization is guaranteed thread-safe, although *
441  * we'd have to use templates or something like that for this multiton pattern. *
442  * *
443  **************************************************************************************************/
444 
445  if (!v) v = AnitaVersion::get();
446 
447  AnitaGeomTool * tmp = instances[v];
448 
449  __asm__ __volatile__ ("" ::: "memory"); //memory fence!
450  if (!tmp)
451  {
452  instance_lock.Lock();
453  tmp = instances[v];
454  if (!tmp)
455  {
456  tmp = new AnitaGeomTool();
457  __asm__ __volatile__ ("" ::: "memory");
458  instances[v] = tmp;
459  }
460  instance_lock.UnLock();
461  }
462 
463  return instances[v];
464 }
465 
466 
467 
468 Int_t AnitaGeomTool::getChanIndex(Int_t surf, Int_t chan){
469  return chan+(CHANNELS_PER_SURF*surf);
470 }
471 
472 
473 
474 void AnitaGeomTool::getCartesianCoords(Double_t lat, Double_t lon, Double_t alt, Double_t p[3])
475 {
476  if(lat<0) lat*=-1;
477  //Note that x and y are switched to conform with previous standards
478  lat*=TMath::DegToRad();
479  lon*=TMath::DegToRad();
480  //calculate x,y,z coordinates
481  Double_t C2=pow(cos(lat)*cos(lat)+(1-FLATTENING_FACTOR)*(1-FLATTENING_FACTOR)*sin(lat)*sin(lat),-0.5);
482  Double_t Q2=(1-FLATTENING_FACTOR)*(1-FLATTENING_FACTOR)*C2;
483  p[1]=(R_EARTH*C2+alt)*TMath::Cos(lat)*TMath::Cos(lon);
484  p[0]=(R_EARTH*C2+alt)*TMath::Cos(lat)*TMath::Sin(lon);
485  p[2]=(R_EARTH*Q2+alt)*TMath::Sin(lat);
486 }
487 
488 void AnitaGeomTool::getLatLonAltFromCartesian(Double_t p[3], Double_t &lat, Double_t &lon, Double_t &alt)
489 {
490  //Here again x and y are flipped for confusions sake
491  Double_t xt=p[1];
492  Double_t yt=p[0];
493  Double_t zt=p[2]; //And flipped z for a test
494 
495  static Double_t cosaeSq=(1-FLATTENING_FACTOR)*(1-FLATTENING_FACTOR);
496  Double_t lonVal=TMath::ATan2(yt,xt);
497  Double_t xySq=TMath::Sqrt(xt*xt+yt*yt);
498  Double_t tanPsit=zt/xySq;
499  Double_t latGuess=TMath::ATan(tanPsit/cosaeSq);
500  Double_t nextLat=latGuess;
501  Double_t geomBot=R_EARTH*R_EARTH*xySq;
502  do {
503  latGuess=nextLat;
504  Double_t N=R_EARTH/TMath::Sqrt(cos(latGuess)*cos(latGuess)+(1-FLATTENING_FACTOR)*(1-FLATTENING_FACTOR)*sin(latGuess)*sin(latGuess));
505  Double_t top=(R_EARTH*R_EARTH*zt + (1-cosaeSq)*cosaeSq*TMath::Power(N*TMath::Sin(latGuess),3));
506  Double_t bottom=geomBot-(1-cosaeSq)*TMath::Power(N*TMath::Cos(latGuess),3);
507  nextLat=TMath::ATan(top/bottom);
508  // std::cout << latGuess << "\t" << nextLat << "\n";
509 
510  } while(TMath::Abs(nextLat-latGuess)>0.0001);
511  latGuess=nextLat;
512  Double_t N=R_EARTH/TMath::Sqrt(cos(latGuess)*cos(latGuess)+(1-FLATTENING_FACTOR)*(1-FLATTENING_FACTOR)*sin(latGuess)*sin(latGuess));
513  Double_t height=(xySq/TMath::Cos(nextLat))-N;
514 
515  lat=latGuess*TMath::RadToDeg();
516  lon=lonVal*TMath::RadToDeg();
517  alt=height;
518  if(lat>0) lat*=-1;
519 
520 }
521 
522 Double_t AnitaGeomTool::getDistanceToCentreOfEarth(Double_t lat)
523 {
524  Double_t pVec[3];
525  this->getCartesianCoords(lat,0,0,pVec);
526 // Double_t cosLat=TMath::Cos(lat);
527 // Double_t sinLat=TMath::Sin(lat);
528 // Double_t a=R_EARTH;
529 // Double_t b=a-FLATTENING_FACTOR*a;
530 // Double_t radSq=(a*a*cosLat)*(a*a*cosLat)+(b*b*sinLat)*(b*b*sinLat);
531 // radSq/=(a*cosLat)*(a*cosLat)+(b*sinLat)*(b*sinLat);
532  // Double_t cosSqAe=(1-FLATTENING_FACTOR)*(1-FLATTENING_FACTOR);
533 // Double_t N=R_EARTH/TMath::Sqrt(cosLat*cosLat+cosSqAe*sinLat*sinLat);
534 // Double_t radSq=N*N*(cosLat*cosLat+cosSqAe*cosSqAe*sinLat*sinLat);
535  return TMath::Sqrt(pVec[0]*pVec[0]+pVec[1]*pVec[1]+pVec[2]*pVec[2]);
536 }
537 
538 
539 // void AnitaGeomTool::getPhiWave(Double_t balloonLon, Double_t balloonLat, Double_t balloonAlt, Double_t balloonHeading, Double_t sourceLon, Double_t sourceLat, Double_t sourceAlt, Double_t &thetaWave, Double_t &phiWave)
540 // {
541 // Double_t thetaBalloon=getThetaFromLat(TMath::Abs(balloonLat));
542 // Double_t phiBalloon=getPhiFromLon(balloonLon);
543 // Double_t balloonHeight=getGeoid(thetaBalloon)+balloonAlt;
544 
545 // Double_t thetaSource=getThetaFromLat(TMath::Abs(sourceLat));
546 // Double_t phiSource=getPhiFromLon(sourceLon);
547 // Double_t radiusSource=getGeoid(thetaSource)+sourceAlt;
548 
549 // //Get vector from Earth's centre to source
550 // TVector3 fSourcePos;
551 // fSourcePos.SetX(radiusSource*TMath::Sin(thetaSource)*TMath::Cos(phiSource));
552 // fSourcePos.SetY(radiusSource*TMath::Sin(thetaSource)*TMath::Sin(phiSource));
553 // fSourcePos.SetZ(radiusSource*TMath::Cos(thetaSource));
554 
555 // //Rotate such that balloon is at 0,0,balloonHeight
556 // fSourcePos.RotateZ(-1*phiBalloon);
557 // fSourcePos.RotateY(-1*thetaBalloon);
558 
559 // //Now find thetaWave and phiWave
560 // thetaWave=TMath::ATan((balloonHeight-fSourcePos.Z())/TMath::Sqrt(fSourcePos.X()*fSourcePos.X() + fSourcePos.Y()*fSourcePos.Y()));
561 
562 // //phiWave is just atan(yp/xp) only looks confusing to make sure I get the sign and 0-360 convention
563 // phiWave=0;
564 // if(fSourcePos.X()==0) {
565 // phiWave=TMath::PiOver2();
566 // if(fSourcePos.Y()<0)
567 // phiWave+=TMath::Pi();
568 // }
569 // else if(fSourcePos.X()<0) {
570 // phiWave=TMath::Pi()+TMath::ATan(fSourcePos.Y()/fSourcePos.X());
571 // }
572 // else {
573 // phiWave=TMath::ATan(fSourcePos.Y()/fSourcePos.X());
574 // if(fSourcePos.Y()<0) {
575 // phiWave+=TMath::TwoPi();
576 // }
577 // }
578 
579 // //Now need to take account of balloon heading
580 // //Will have to check heading at some point
581 // if(balloonHeading>=0 && balloonHeading<=360) {
582 // phiWave+=balloonHeading*TMath::DegToRad();
583 // if(phiWave>TMath::TwoPi())
584 // phiWave-=TMath::TwoPi();
585 // }
586 
587 // }
588 
589 Double_t AnitaGeomTool::getPhiDiff(Double_t firstPhi, Double_t secondPhi)
590 {
591  Double_t phiDiff=firstPhi-secondPhi;
592  if(TMath::Abs(phiDiff)>TMath::Abs(phiDiff+TMath::TwoPi()))
593  phiDiff+=TMath::TwoPi();
594  if(TMath::Abs(phiDiff)>TMath::Abs(phiDiff-TMath::TwoPi()))
595  phiDiff-=TMath::TwoPi();
596  return phiDiff;
597 }
598 
599 
601  Int_t phi,
603 {
604  Int_t ant;
605  if(phi<0 || phi>15) return -1;
606  switch(ring) {
607  case AnitaRing::kTopRing:
608  ant=AnitaGeom::topAntNums[phi];
609  break;
611  ant=AnitaGeom::middleAntNums[phi];
612  break;
614  ant=AnitaGeom::bottomAntNums[phi];
615  break;
616  default:
617  return -1;
618  }
619  return getChanIndexFromAntPol(ant,pol);
620 }
621 
622 
625 {
626  if(ant<0 || ant>(NUM_SEAVEYS-1)) return -1;
627  Int_t surf,chan;
628  surf=AnitaGeom::antToSurfMap[ant];
629  if(pol==AnitaPol::kHorizontal)
630  chan=AnitaGeom::hAntToChan[ant];
631  else if(pol==AnitaPol::kVertical)
632  chan=AnitaGeom::vAntToChan[ant];
633  else
634  return -1;
635  return getChanIndex(surf,chan);
636 }
637 
638 
639 // Int_t AnitaGeomTool::getRFCMFromAnt(Int_t ant)
640 // {
641 // if(ant<0 || ant>(NUM_SEAVEYS-1)) return -1;
642 // return AnitaGeom::antToRFCM[ant];
643 // }
644 
645 // Int_t AnitaGeomTool::getRFCMChannelFromAntPol(Int_t ant,AnitaPol::AnitaPol_t pol)
646 // {
647 // if(ant<0 || ant>(NUM_SEAVEYS-1)) return -1;
648 // if(pol==AnitaPol::kHorizontal)
649 // return AnitaGeom::hAntToRFCMChan[ant];
650 // if(pol==AnitaPol::kVertical)
651 // return AnitaGeom::vAntToRFCMChan[ant];
652 
653 // return -1;
654 // }
655 
656 
657 Int_t AnitaGeomTool::getSurfFromAnt(Int_t ant)
658 {
659  if(ant<0 || ant>(NUM_SEAVEYS-1)) return -1;
660  Int_t surf;
661  surf=AnitaGeom::antToSurfMap[ant];
662  return surf;
663 }
664 
665 
666 Int_t AnitaGeomTool::getChanFromAntPol(Int_t ant,AnitaPol::AnitaPol_t pol)
667 {
668  Int_t chan;
669  if(pol==AnitaPol::kHorizontal)
670  chan=AnitaGeom::hAntToChan[ant];
671  else if(pol==AnitaPol::kVertical)
672  chan=AnitaGeom::vAntToChan[ant];
673  else
674  return -1;
675  return chan;
676 }
677 
678 
680 {
681  //Need to update for bottom ring
682 
683  if (rx<16)
684  return AnitaGeom::middleAntNums[AnitaGeom::topPhiNums[rx]];
685  else if(rx<32)
686  return AnitaGeom::topAntNums[AnitaGeom::middlePhiNums[rx-16]];
687  else
688  return AnitaGeom::middleAntNums[AnitaGeom::bottomPhiNums[rx-32]];
689 
690  return -1;
691 }
692 
693 void AnitaGeomTool::getNeighbors(Int_t rx,int& rxleft,int& rxright)
694 {
695 
696  /* for A3 and A4, for the top ring, the neighbors are offset by 2 and in the
697  * bottom / mid ring they are offset by 1
698  */
699 
700  int ring = rx / 16;
701 
702  if (ring == 0)
703  {
704  rxleft = (rx + 14)%16;
705  rxright = (rx + 2)%16;
706  }
707 
708  else
709  {
710  int ant = rx % 16;
711  rxleft = ring * 16 + (ant +15) % 16;
712  rxright = ring * 16 + (ant +1) % 16;
713  }
714 
715 
716 }
717 
718 void AnitaGeomTool::getThetaPartners(Int_t rx,int& rxleft,int& rxright)
719 {
720  if (rx<0 || rx>(NUM_SEAVEYS-1))
721  std::cerr << "Antenna number out of range!\n";
722  if (rx<16) {
723  Int_t phi=AnitaGeom::topPhiNums[rx];
724  Int_t phiLeft=phi-1;
725  if(phiLeft<0) phiLeft=15;
726  Int_t phiRight=phi+1;
727  if(phiRight>15) phiRight=0;
728  rxleft=AnitaGeom::topAntNums[phiLeft];
729  rxright=AnitaGeom::topAntNums[phiRight];
730  }
731  else if (rx<32){
732  Int_t phi=AnitaGeom::middlePhiNums[rx-16];
733  Int_t phiLeft=phi-1;
734  if(phiLeft<0) phiLeft=15;
735  Int_t phiRight=phi+1;
736  if(phiRight>15) phiRight=0;
737  rxleft=AnitaGeom::middleAntNums[phiLeft];
738  rxright=AnitaGeom::middleAntNums[phiRight];
739  }
740  else{
741  Int_t phi=AnitaGeom::bottomPhiNums[rx-32];
742  Int_t phiLeft=phi-1;
743  if(phiLeft<0) phiLeft=15;
744  Int_t phiRight=phi+1;
745  if(phiRight>15) phiRight=0;
746  rxleft=AnitaGeom::bottomAntNums[phiLeft];
747  rxright=AnitaGeom::bottomAntNums[phiRight];
748  }
749 
750 }
751 
753 {
754  if (rx<16)
755  return AnitaGeom::topPhiNums[rx];
756  else if(rx<32)
757  return AnitaGeom::middlePhiNums[rx-16];
758  else if(rx<48) {
759  return AnitaGeom::bottomPhiNums[rx-32];
760  }
761 
762  return 0;
763 }
764 
765 Double_t AnitaGeomTool::getDirectionWrtNorth(Int_t phiSector, Double_t heading) {
766 
767  if(phiSector<0 || phiSector>=NUM_PHI) return -1.;
768 
769  //Copied this change from Amy don't (yet) know if it is correct
770  // Float_t direction=(1.*360./NUM_PHI)+heading-(phiSector*360./NUM_PHI);
771 
773  Double_t direction = heading + aftForeOffsetAngleVertical*TMath::RadToDeg() - (phiSector*360./NUM_PHI);
774 
775  if(direction>=360){
776  direction -= 360;
777  }
778  else if(direction < 0){
779  direction += 360;
780  }
781 
782  return direction;
783 
784 }
785 
786 
787 
788 Int_t AnitaGeomTool::getSurfChanFromChanIndex(Int_t chanIndex, // input channel index
789  Int_t &surf,Int_t &chan) // output surf and channel
790 {
791  if(chanIndex<0 || chanIndex>=NUM_DIGITZED_CHANNELS)
792  return 0;
793 
794  // std::cout << "chanIndex is " << chanIndex << "\n";
795  chan=chanIndex%CHANNELS_PER_SURF; // 0 to 8
796  surf=(chanIndex-chan)/CHANNELS_PER_SURF; // 0 to 8
797 
798  return 1;
799 
800 }
801 Int_t AnitaGeomTool::getAntPolFromSurfChan(Int_t surf,Int_t chan,Int_t &ant, AnitaPol::AnitaPol_t &pol)
802 {
803 
804  int v = AnitaVersion::get();
805 
806  if (v >4 && v < 3) return 0;
807 
808  ant= v == 4 ? AnitaGeom::surfToAntMapAnita4[surf][chan] :
809  v == 3 ? AnitaGeom::surfToAntMapAnita3[surf][chan] :
810  0 ;
811 
812 
813  // std::cout << "surf, chan, ant are " << surf << " " << chan << "\n";
814 
815  if (ant>0)
817  else
819 
820  ant=abs(ant)-1; // so that it's from 0 to 31
821 
822  return 1;
823 }
824 
826  if(ant<0 || ant>=NUM_SEAVEYS)
827  return 0;
828 
829  return AnitaVersion::get() == 3 ? AnitaGeom::antOrientationMapAnita3[ant] : AnitaGeom::antOrientationMap[ant];
830 }
831 
832 Int_t AnitaGeomTool::getLayer(Int_t irx)
833 {
834  if (irx<8)
835  return 0;
836  else if (irx>=8 && irx<16)
837  return 1;
838  else if (irx>=16)
839  return 2;
840 
841  return 0;
842 
843 }
844 
846  if (ant<8)
847  return AnitaRing::kTopRing;
848  else if (ant>=8 && ant<16)
849  return AnitaRing::kTopRing;
850  else if (ant<32)
851  return AnitaRing::kMiddleRing;
852  else if (ant<48)
853  return AnitaRing::kBottomRing;
854  return AnitaRing::kNotARing;
855 
856 }
857 
860  Int_t &ant,
862  Int_t &phi)
863 {
864  AnitaGeomTool::getAntPolFromSurfChan(surf,chan,ant,pol);
867 }
868 
870  Int_t phi,
872  Int_t &surf, Int_t &chan, Int_t &ant)
873 {
874  Int_t chanIndex=AnitaGeomTool::getChanIndexFromRingPhiPol(ring,phi,pol);
875  AnitaGeomTool::getSurfChanFromChanIndex(chanIndex,surf,chan);
876  ant=AnitaGeomTool::getAntFromPhiRing(phi,ring);
877 }
878 
879 
880 void AnitaGeomTool::readAnita2Photogrammetry()
881 {
882  // return; //Hack for now
883  char calibDir[FILENAME_MAX];
884  char fileName[FILENAME_MAX];
885  char *calibEnv=getenv("ANITA_CALIB_DIR");
886  if(!calibEnv) {
887  char *utilEnv=getenv("ANITA_UTIL_INSTALL_DIR");
888  if(!utilEnv)
889  sprintf(calibDir,"calib");
890  else
891  sprintf(calibDir,"%s/share/anitaCalib",utilEnv);
892  }
893  else {
894  strncpy(calibDir,calibEnv,FILENAME_MAX);
895  }
896 
897  sprintf(fileName,"%s/antennaPositionTables.csv",calibDir);
898  std::ifstream AntennaFile(fileName);
899  if(!AntennaFile) {
900  std::cerr << "Couldn't open:\t" << fileName << "\n";
901  return;
902  }
903 
904 
905  //First up are the antenna positions
906  TString line;
907  for(Int_t i=0;i<4;i++) {
908  line.ReadLine(AntennaFile);
909  //std::cout << line.Data() << "\n";
910  }
911 
912  for(Int_t ant=0;ant<32;ant++) {
913  line.ReadLine(AntennaFile);
914  //std::cout << "Seavey:\t" << line.Data() << "\n";
915  TObjArray *tokens = line.Tokenize(",");
916  for(Int_t j=0;j<8;j++) {
917  const TString subString = ((TObjString*)tokens->At(j))->GetString();
918  // TString *subString = (TString*) tokens->At(j);
919  // std::cout << j << "\t" << subString.Data() << "\n";
920  switch(j) {
921  case 0:
922  if (ant+1 != subString.Atoi()) {
923  std::cerr << "Antenna number mismatch\n";
924  }
925  break;
926  case 1:
927  xAntFromDeckHorn[ant]=subString.Atof()*INCHTOMETER; //m
928  break;
929  case 2:
930  yAntFromDeckHorn[ant]=subString.Atof()*INCHTOMETER; //m
931  break;
932  case 3:
933  zAntFromDeckHorn[ant]=subString.Atof()*INCHTOMETER; //m
934  break;
935  case 4:
936  rAntFromDeckHorn[ant]=subString.Atof()*INCHTOMETER; //m
937  // std::cout << ant << "\t" << rAntFromDeckHorn[ant] << "\n";
938  break;
939  case 5:
940  azCentreFromDeckHorn[ant]=subString.Atof()*TMath::DegToRad(); //radians
941  break;
942  case 6:
943  apertureAzFromDeckHorn[ant]=subString.Atof()*TMath::DegToRad(); //radians
944  break;
945  case 7:
946  apertureElFromDeckHorn[ant]=subString.Atof()*TMath::DegToRad(); //radians
947  break;
948  default:
949  break;
950  }
951 
952  }
953  tokens->Delete();
954 
955  }
956 
957 
958  //Now we'll add a hack for the drop-down antennas
959  for(Int_t ant=32;ant<48;ant++) {
960  rAntFromDeckHorn[ant]=(101*INCHTOMETER)-0.38;
961  zAntFromDeckHorn[ant]=(-21*INCHTOMETER)-1.8;
962  azCentreFromDeckHorn[ant]=(-67.5 + 45*(ant-32))*TMath::DegToRad();
963  apertureAzFromDeckHorn[ant]=azCentreFromDeckHorn[ant];
964  apertureElFromDeckHorn[ant]=-10*TMath::DegToRad();
965  xAntFromDeckHorn[ant]=rAntFromDeckHorn[ant]*TMath::Cos(azCentreFromDeckHorn[ant]);
966  yAntFromDeckHorn[ant]=rAntFromDeckHorn[ant]*TMath::Sin(azCentreFromDeckHorn[ant]);
967  // std::cout << ant << "\t" << xAntFromDeckHorn[ant] << "\t" << yAntFromDeckHorn[ant]
968  // << "\t" << zAntFromDeckHorn[ant] << "\t" << apertureAzFromDeckHorn[ant]*TMath::RadToDeg() << "\n";
969  }
970  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
971  for(Int_t ant=0;ant<NUM_SEAVEYS;ant++) {
972  Double_t phaseCentreToAntFront=ringPhaseCentreOffset[Int_t(this->getRingFromAnt(ant))];
973 
974 
975  Double_t deltaXRL = -deltaRL*TMath::Sin(apertureAzFromDeckHorn[ant]);
976  Double_t deltaYRL = deltaRL*TMath::Cos(apertureAzFromDeckHorn[ant]);
977 
978  Double_t deltaZUD = deltaUD*TMath::Cos(apertureElFromDeckHorn[ant]);
979  Double_t deltaXUD = deltaUD*TMath::Sin(apertureElFromDeckHorn[ant])*TMath::Cos(apertureAzFromDeckHorn[ant]);
980  Double_t deltaYUD = deltaUD*TMath::Sin(apertureElFromDeckHorn[ant])*TMath::Sin(apertureAzFromDeckHorn[ant]);
981 
982  deltaXRL = 0;
983  deltaYRL = 0;
984 
985  deltaZUD = 0;
986  deltaXUD = 0;
987  deltaYUD = 0;
988 
989  // std::cout << deltaXRL << " " << deltaXUD << " "<< deltaYUD << " "<< deltaYRL << " " << deltaZUD << std::endl;
990 
991  // Now try to find the location of the antenna phaseCentre poInt_t
992  xPhaseCentreFromDeckHorn[ant][pol]=(xAntFromDeckHorn[ant] + deltaXRL + deltaXUD) -
993  phaseCentreToAntFront*TMath::Cos( apertureAzFromDeckHorn[ant])*TMath::Cos(apertureElFromDeckHorn[ant]); //m
994  yPhaseCentreFromDeckHorn[ant][pol]=(yAntFromDeckHorn[ant] + deltaYRL + deltaYUD) -
995  phaseCentreToAntFront*TMath::Sin( apertureAzFromDeckHorn[ant])*TMath::Cos(apertureElFromDeckHorn[ant]); //m
996  zPhaseCentreFromDeckHorn[ant][pol]=(zAntFromDeckHorn[ant] + deltaZUD) -phaseCentreToAntFront*TMath::Sin(apertureElFromDeckHorn[ant]); //m
997 
998  // xPhaseCentreFromDeckHorn[ant]=xAntFromDeckHorn[ant] -
999  // phaseCentreToAntFront*TMath::Cos( apertureAzFromDeckHorn[ant])*TMath::Cos(apertureElFromDeckHorn[ant]); //m
1000  // yPhaseCentreFromDeckHorn[ant]=yAntFromDeckHorn[ant] -
1001  // phaseCentreToAntFront*TMath::Sin( apertureAzFromDeckHorn[ant])*TMath::Cos(apertureElFromDeckHorn[ant]); //m
1002  // zPhaseCentreFromDeckHorn[ant]=zAntFromDeckHorn[ant] - phaseCentreToAntFront*TMath::Sin(apertureElFromDeckHorn[ant]); //m
1003 
1004  rPhaseCentreFromDeckHorn[ant][pol]=TMath::Sqrt(xPhaseCentreFromDeckHorn[ant][pol]*xPhaseCentreFromDeckHorn[ant][pol]+
1005  yPhaseCentreFromDeckHorn[ant][pol]*yPhaseCentreFromDeckHorn[ant][pol]);//m
1006  azPhaseCentreFromDeckHorn[ant][pol]=TMath::ATan(yPhaseCentreFromDeckHorn[ant][pol]/xPhaseCentreFromDeckHorn[ant][pol]); //radians
1007  if(xPhaseCentreFromDeckHorn[ant][pol]<0)
1008  azPhaseCentreFromDeckHorn[ant][pol]+=TMath::Pi();
1009  if(azPhaseCentreFromDeckHorn[ant][pol]<0)
1010  azPhaseCentreFromDeckHorn[ant][pol]+=TMath::TwoPi();
1011 
1012  // std::cout << "Face v PhaseCentre -- Ant " << ant << std::endl;
1013  // std::cout << apertureAzFromDeckHorn[ant] << "\t" << apertureElFromDeckHorn[ant] << "\n";
1014  // std::cout << "x\t" << xAntFromDeckHorn[ant] << "\t" << xPhaseCentreFromDeckHorn[ant][pol] << std::endl;
1015  // std::cout << "y\t" << yAntFromDeckHorn[ant] << "\t" << yPhaseCentreFromDeckHorn[ant][pol] << std::endl;
1016  // std::cout << "z\t" << zAntFromDeckHorn[ant] << "\t" << zPhaseCentreFromDeckHorn[ant][pol] << std::endl;
1017  // std::cout << "r\t" << rAntFromDeckHorn[ant] << "\t" << rPhaseCentreFromDeckHorn[ant][pol] << std::endl;
1018  // std::cout << "phi\t" << azCentreFromDeckHorn[ant] << "\t" << azPhaseCentreFromDeckHorn[ant][pol] << std::endl;
1019 
1020  }
1021  }
1022 
1023  //Now bicones and discones
1024  for(Int_t i=0;i<2;i++) {
1025  line.ReadLine(AntennaFile);
1026  //std::cout << line.Data() << "\n";
1027  }
1028 
1029 
1030  for(Int_t ant=0;ant<4;ant++) {
1031  line.ReadLine(AntennaFile);
1032  //std::cout << "Bicones:\t" << line.Data() << "\n";
1033  TObjArray *tokens = line.Tokenize(",");
1034  for(Int_t j=0;j<4;j++) {
1035  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1036  // TString *subString = (TString*) tokens->At(j);
1037  // std::cout << j << "\t" << subString.Data() << "\n";
1038  switch(j) {
1039  case 0:
1040  break;
1041  case 1:
1042  xAntFromDeckBicone[ant]=subString.Atof()*INCHTOMETER; //m
1043  break;
1044  case 2:
1045  yAntFromDeckBicone[ant]=subString.Atof()*INCHTOMETER; //m
1046  break;
1047  case 3:
1048  zAntFromDeckBicone[ant]=subString.Atof()*INCHTOMETER; //m
1049  break;
1050  default:
1051  break;
1052  }
1053 
1054  }
1055  tokens->Delete();
1056  }
1057 
1058  for(Int_t ant=0;ant<4;ant++) {
1059  line.ReadLine(AntennaFile);
1060  //std::cout << "Discones:\t" << line.Data() << "\n";
1061  TObjArray *tokens = line.Tokenize(",");
1062  for(Int_t j=0;j<4;j++) {
1063  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1064  // TString *subString = (TString*) tokens->At(j);
1065  // std::cout << j << "\t" << subString.Data() << "\n";
1066  switch(j) {
1067  case 0:
1068  break;
1069  case 1:
1070  xAntFromDeckDiscone[ant]=subString.Atof()*INCHTOMETER; //m
1071  break;
1072  case 2:
1073  yAntFromDeckDiscone[ant]=subString.Atof()*INCHTOMETER; //m
1074  break;
1075  case 3:
1076  zAntFromDeckDiscone[ant]=subString.Atof()*INCHTOMETER; //m
1077  break;
1078  default:
1079  break;
1080  }
1081  }
1082  tokens->Delete();
1083  }
1084 
1085  //Now box enclosures and the like
1086  for(Int_t i=0;i<2;i++) {
1087  line.ReadLine(AntennaFile);
1088  //std::cout << line.Data() << "\n";
1089  }
1090  for(Int_t corner=0;corner<4;corner++) {
1091  line.ReadLine(AntennaFile);
1092  //std::cout << "Anita Box:\t" << line.Data() << "\n";
1093  TObjArray *tokens = line.Tokenize(",");
1094  for(Int_t j=0;j<4;j++) {
1095  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1096  // TString *subString = (TString*) tokens->At(j);
1097  // std::cout << j << "\t" << subString.Data() << "\n";
1098  switch(j) {
1099  case 0:
1100  break;
1101  case 1:
1102  xAnitaBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1103  break;
1104  case 2:
1105  yAnitaBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1106  break;
1107  case 3:
1108  zAnitaBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1109  break;
1110  default:
1111  break;
1112  }
1113 
1114  }
1115  tokens->Delete();
1116  }
1117  for(Int_t i=0;i<2;i++) {
1118  line.ReadLine(AntennaFile);
1119  //std::cout << line.Data() << "\n";
1120  }
1121  for(Int_t corner=0;corner<4;corner++) {
1122  line.ReadLine(AntennaFile);
1123  //std::cout << "Battery Box:\t" << line.Data() << "\n";
1124  TObjArray *tokens = line.Tokenize(",");
1125  for(Int_t j=0;j<4;j++) {
1126  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1127  // TString *subString = (TString*) tokens->At(j);
1128  // std::cout << j << "\t" << subString.Data() << "\n";
1129  switch(j) {
1130  case 0:
1131  break;
1132  case 1:
1133  xBatteryBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1134  break;
1135  case 2:
1136  yBatteryBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1137  break;
1138  case 3:
1139  zBatteryBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1140  break;
1141  default:
1142  break;
1143  }
1144 
1145  }
1146  tokens->Delete();
1147  }
1148  for(Int_t i=0;i<2;i++) {
1149  line.ReadLine(AntennaFile);
1150  //std::cout << line.Data() << "\n";
1151  }
1152  for(Int_t corner=0;corner<4;corner++) {
1153  line.ReadLine(AntennaFile);
1154  //std::cout << "Sip Box:\t" << line.Data() << "\n";
1155  TObjArray *tokens = line.Tokenize(",");
1156  for(Int_t j=0;j<4;j++) {
1157  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1158  // TString *subString = (TString*) tokens->At(j);
1159  // std::cout << j << "\t" << subString.Data() << "\n";
1160  switch(j) {
1161  case 0:
1162  break;
1163  case 1:
1164  xSipBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1165  break;
1166  case 2:
1167  ySipBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1168  break;
1169  case 3:
1170  zSipBoxFromDeckCorner[corner]=subString.Atof()*INCHTOMETER; //m
1171  break;
1172  default:
1173  break;
1174  }
1175 
1176  }
1177  tokens->Delete();
1178  }
1179  //Now gps stuff
1180  for(Int_t i=0;i<2;i++) {
1181  line.ReadLine(AntennaFile);
1182  //std::cout << line.Data() << "\n";
1183  }
1184  {
1185  line.ReadLine(AntennaFile);
1186  //std::cout << "GPS Plane:\t" << line.Data() << "\n";
1187  TObjArray *tokens = line.Tokenize(",");
1188  for(Int_t j=0;j<4;j++) {
1189  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1190  // TString *subString = (TString*) tokens->At(j);
1191  // std::cout << j << "\t" << subString.Data() << "\n";
1192  switch(j) {
1193  case 0:
1194  break;
1195  case 1:
1196  case 2:
1197  case 3:
1198  gpsPlaneFromDeck[j-1]=subString.Atof(); //unit vector
1199  break;
1200  default:
1201  break;
1202  }
1203 
1204  }
1205  tokens->Delete();
1206  }
1207  {
1208  line.ReadLine(AntennaFile);
1209  //std::cout << "GPS Heading:\t" << line.Data() << "\n";
1210  TObjArray *tokens = line.Tokenize(",");
1211  for(Int_t j=0;j<4;j++) {
1212  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1213  // TString *subString = (TString*) tokens->At(j);
1214  // std::cout << j << "\t" << subString.Data() << "\n";
1215  switch(j) {
1216  case 0:
1217  break;
1218  case 1:
1219  case 2:
1220  case 3:
1221  gpsHeadingFromDeck[j-1]=subString.Atof(); //unit vector
1222  break;
1223  default:
1224  break;
1225  }
1226 
1227  }
1228  tokens->Delete();
1229  }
1230 
1231  aftForeOffsetAngleDeck=TMath::ATan(gpsHeadingFromDeck[1]/gpsHeadingFromDeck[0]);
1232 
1233  //Now switch to the (more useful) vertical measurements
1234  //First up are the antenna positions
1235  for(Int_t i=0;i<5;i++) {
1236  line.ReadLine(AntennaFile);
1237  //std::cout << line.Data() << "\n";
1238  }
1239 
1240  for(Int_t ant=0;ant<32;ant++) {
1241  line.ReadLine(AntennaFile);
1242  //std::cout << "Seavey:\t" << line.Data() << "\n";
1243  TObjArray *tokens = line.Tokenize(",");
1244  for(Int_t j=0;j<8;j++) {
1245  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1246  // TString *subString = (TString*) tokens->At(j);
1247  // //std::cout << j << "\t" << subString.Data() << "\n";
1248  switch(j) {
1249  case 0:
1250  if (ant+1 != subString.Atoi()) {
1251  std::cerr << "Antenna number mismatch\n";
1252  }
1253  break;
1254  case 1:
1255  xAntFromVerticalHorn[ant]=subString.Atof()*INCHTOMETER; //m
1256  break;
1257  case 2:
1258  yAntFromVerticalHorn[ant]=subString.Atof()*INCHTOMETER; //m
1259  break;
1260  case 3:
1261  zAntFromVerticalHorn[ant]=subString.Atof()*INCHTOMETER; //m
1262  break;
1263  case 4:
1264  rAntFromVerticalHorn[ant]=subString.Atof()*INCHTOMETER; //m
1265  break;
1266  case 5:
1267  azCentreFromVerticalHorn[ant]=subString.Atof()*TMath::DegToRad(); //radians
1268  if(azCentreFromVerticalHorn[ant]<0)
1269  azCentreFromVerticalHorn[ant]+=TMath::TwoPi();
1270  break;
1271  case 6:
1272  apertureAzFromVerticalHorn[ant]=subString.Atof()*TMath::DegToRad(); //radians
1273  break;
1274  case 7:
1275  apertureElFromVerticalHorn[ant]=subString.Atof()*TMath::DegToRad(); //radians
1276  break;
1277  default:
1278  break;
1279  }
1280 
1281  }
1282  tokens->Delete();
1283  }
1284  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
1285  for(Int_t ant=0;ant<32;ant++) {
1286 
1287  Double_t phaseCentreToAntFront=ringPhaseCentreOffset[Int_t(this->getRingFromAnt(ant))];
1288 
1289  Double_t deltaXRL = -deltaRL*TMath::Sin(apertureAzFromVerticalHorn[ant]);
1290  Double_t deltaYRL = deltaRL*TMath::Cos(apertureAzFromVerticalHorn[ant]);
1291 
1292  Double_t deltaZUD = deltaUD*TMath::Cos(apertureElFromVerticalHorn[ant]);
1293  Double_t deltaXUD = -deltaUD*TMath::Sin(apertureElFromVerticalHorn[ant])*TMath::Cos(apertureAzFromVerticalHorn[ant]);
1294  Double_t deltaYUD = -deltaUD*TMath::Sin(apertureElFromVerticalHorn[ant])*TMath::Sin(apertureAzFromVerticalHorn[ant]);
1295 
1296  //std::cout << deltaXRL << " " << deltaXUD << " "<< deltaYUD << " "<< deltaYRL << " " << deltaZUD << std::endl;
1297 
1298  //Now try to find the location of the antenna phaseCentre poInt_t
1299  xPhaseCentreFromVerticalHorn[ant][pol]=(xAntFromVerticalHorn[ant] + deltaXRL + deltaXUD) -
1300  phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1301  yPhaseCentreFromVerticalHorn[ant][pol]=(yAntFromVerticalHorn[ant] + deltaYRL + deltaYUD) -
1302  phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1303  zPhaseCentreFromVerticalHorn[ant][pol]=(zAntFromVerticalHorn[ant] + deltaZUD) -phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1304 
1305  // //Now try to find the location of the antenna phaseCentre poInt_t
1306  // xPhaseCentreFromVerticalHorn[ant]=xAntFromVerticalHorn[ant]-
1307  // phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1308  // yPhaseCentreFromVerticalHorn[ant]=yAntFromVerticalHorn[ant]-
1309  // phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1310  // zPhaseCentreFromVerticalHorn[ant]=zAntFromVerticalHorn[ant]-phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1311  rPhaseCentreFromVerticalHorn[ant][pol]=TMath::Sqrt(xPhaseCentreFromVerticalHorn[ant][pol]*xPhaseCentreFromVerticalHorn[ant][pol]+
1312  yPhaseCentreFromVerticalHorn[ant][pol]*yPhaseCentreFromVerticalHorn[ant][pol]);//m
1313  azPhaseCentreFromVerticalHorn[ant][pol]=TMath::ATan(yPhaseCentreFromVerticalHorn[ant][pol]/xPhaseCentreFromVerticalHorn[ant][pol]); //radians
1314  if(xPhaseCentreFromVerticalHorn[ant][pol]<0)
1315  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::Pi();
1316  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
1317  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
1318 
1319 
1320  // std::cout << "Face v PhaseCentre -- Ant " << ant << std::endl;
1321  // std::cout << "x\t" << xAntFromVerticalHorn[ant] << "\t" << xPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
1322  // std::cout << "y\t" << yAntFromVerticalHorn[ant] << "\t" << yPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
1323  // std::cout << "z\t" << zAntFromVerticalHorn[ant] << "\t" << zPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
1324  // std::cout << "r\t" << rAntFromVerticalHorn[ant] << "\t" << rPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
1325  // std::cout << "phi\t" << azCentreFromVerticalHorn[ant] << "\t" << azPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
1326 
1327  }
1328  }
1329 
1330 
1331 
1332  //Now we'll add a hack for the drop-down antennas
1333  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
1334  for(Int_t ant=32;ant<NUM_SEAVEYS;ant++) {
1335  Double_t phaseCentreToAntFront=ringPhaseCentreOffset[Int_t(this->getRingFromAnt(ant))];
1336  // std::cout << ant << "\t" << phaseCentreToAntFront << "\n";
1337  Double_t deltaXRL = 0;
1338  Double_t deltaYRL = 0;
1339 
1340  Double_t deltaZUD = 0;
1341  Double_t deltaXUD = 0;
1342  Double_t deltaYUD = 0;
1343 
1344  // RJN changed the 0.38 number to 0.3 to match simons update function
1345  //Here's the little bastard that was causing all the problems
1346  rAntFromVerticalHorn[ant]=(101*INCHTOMETER)-0.3;
1347  // rAntFromVerticalHorn[ant]=(101*INCHTOMETER)-0.38;
1348  zAntFromVerticalHorn[ant]=(-21*INCHTOMETER)-1.8;
1349  azCentreFromVerticalHorn[ant]=(-67.5 + 45*(ant-32))*TMath::DegToRad();
1350 
1351  apertureAzFromVerticalHorn[ant]=azCentreFromVerticalHorn[ant];
1352  apertureElFromVerticalHorn[ant]=-10*TMath::DegToRad();
1353  xAntFromVerticalHorn[ant]=rAntFromVerticalHorn[ant]*TMath::Cos(azCentreFromVerticalHorn[ant]);
1354  yAntFromVerticalHorn[ant]=rAntFromVerticalHorn[ant]*TMath::Sin(azCentreFromVerticalHorn[ant]);
1355 
1356 
1357  xPhaseCentreFromVerticalHorn[ant][pol]=(xAntFromVerticalHorn[ant] + deltaXRL + deltaXUD) -
1358  phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1359  yPhaseCentreFromVerticalHorn[ant][pol]=(yAntFromVerticalHorn[ant] + deltaYRL + deltaYUD) -
1360  phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1361  zPhaseCentreFromVerticalHorn[ant][pol]=(zAntFromVerticalHorn[ant] + deltaZUD) -phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1362 
1363  // //Now try to find the location of the antenna phaseCentre poInt_t
1364  // xPhaseCentreFromVerticalHorn[ant][pol]=xAntFromVerticalHorn[ant]-
1365  // phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1366  // yPhaseCentreFromVerticalHorn[ant][pol]=yAntFromVerticalHorn[ant]-
1367  // phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1368  // zPhaseCentreFromVerticalHorn[ant][pol]=zAntFromVerticalHorn[ant]-phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1369  rPhaseCentreFromVerticalHorn[ant][pol]=TMath::Sqrt(xPhaseCentreFromVerticalHorn[ant][pol]*xPhaseCentreFromVerticalHorn[ant][pol]+
1370  yPhaseCentreFromVerticalHorn[ant][pol]*yPhaseCentreFromVerticalHorn[ant][pol]);//m
1371  azPhaseCentreFromVerticalHorn[ant][pol]=TMath::ATan(yPhaseCentreFromVerticalHorn[ant][pol]/xPhaseCentreFromVerticalHorn[ant][pol]); //radians
1372  if(xPhaseCentreFromVerticalHorn[ant][pol]<0)
1373  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::Pi();
1374  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
1375  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
1376 
1377 
1378  // std::cout << "Ryan's:\t" << ant << " "
1379 // << rPhaseCentreFromVerticalHorn[ant][pol] << " "
1380 // << azPhaseCentreFromVerticalHorn[ant][pol] << " "
1381 // << zPhaseCentreFromVerticalHorn[ant][pol] << "\n";
1382 
1383  }
1384  }
1385 
1386 
1387  //Now bicones and discones
1388  for(Int_t i=0;i<2;i++) {
1389  line.ReadLine(AntennaFile);
1390  //std::cout << line.Data() << "\n";
1391  }
1392 
1393 
1394  for(Int_t ant=0;ant<4;ant++) {
1395  line.ReadLine(AntennaFile);
1396  //std::cout << "Bicones:\t" << line.Data() << "\n";
1397  TObjArray *tokens = line.Tokenize(",");
1398  for(Int_t j=0;j<4;j++) {
1399  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1400  // TString *subString = (TString*) tokens->At(j);
1401  // //std::cout << j << "\t" << subString.Data() << "\n";
1402  switch(j) {
1403  case 0:
1404  break;
1405  case 1:
1406  xAntFromVerticalBicone[ant]=subString.Atof()*INCHTOMETER; //m
1407  break;
1408  case 2:
1409  yAntFromVerticalBicone[ant]=subString.Atof()*INCHTOMETER; //m
1410  break;
1411  case 3:
1412  zAntFromVerticalBicone[ant]=subString.Atof()*INCHTOMETER; //m
1413  break;
1414  default:
1415  break;
1416  }
1417 
1418  }
1419  tokens->Delete();
1420  }
1421 
1422  for(Int_t ant=0;ant<4;ant++) {
1423  line.ReadLine(AntennaFile);
1424  //std::cout << "Discones:\t" << line.Data() << "\n";
1425  TObjArray *tokens = line.Tokenize(",");
1426  for(Int_t j=0;j<4;j++) {
1427  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1428  // TString *subString = (TString*) tokens->At(j);
1429  // //std::cout << j << "\t" << subString.Data() << "\n";
1430  switch(j) {
1431  case 0:
1432  break;
1433  case 1:
1434  xAntFromVerticalDiscone[ant]=subString.Atof()*INCHTOMETER; //m
1435  break;
1436  case 2:
1437  yAntFromVerticalDiscone[ant]=subString.Atof()*INCHTOMETER; //m
1438  break;
1439  case 3:
1440  zAntFromVerticalDiscone[ant]=subString.Atof()*INCHTOMETER; //m
1441  break;
1442  default:
1443  break;
1444  }
1445 
1446  }
1447  tokens->Delete();
1448  }
1449 
1450  //Now box enclosures and the like
1451  for(Int_t i=0;i<2;i++) {
1452  line.ReadLine(AntennaFile);
1453  //std::cout << line.Data() << "\n";
1454  }
1455  for(Int_t corner=0;corner<4;corner++) {
1456  line.ReadLine(AntennaFile);
1457  //std::cout << "Anita Box:\t" << line.Data() << "\n";
1458  TObjArray *tokens = line.Tokenize(",");
1459  for(Int_t j=0;j<4;j++) {
1460  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1461  // TString *subString = (TString*) tokens->At(j);
1462  // //std::cout << j << "\t" << subString.Data() << "\n";
1463  switch(j) {
1464  case 0:
1465  break;
1466  case 1:
1467  xAnitaBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1468  break;
1469  case 2:
1470  yAnitaBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1471  break;
1472  case 3:
1473  zAnitaBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1474  break;
1475  default:
1476  break;
1477  }
1478 
1479  }
1480  tokens->Delete();
1481  }
1482  for(Int_t i=0;i<2;i++) {
1483  line.ReadLine(AntennaFile);
1484  //std::cout << line.Data() << "\n";
1485  }
1486  for(Int_t corner=0;corner<4;corner++) {
1487  line.ReadLine(AntennaFile);
1488  //std::cout << "Battery Box:\t" << line.Data() << "\n";
1489  TObjArray *tokens = line.Tokenize(",");
1490  for(Int_t j=0;j<4;j++) {
1491  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1492  // TString *subString = (TString*) tokens->At(j);
1493  // //std::cout << j << "\t" << subString.Data() << "\n";
1494  switch(j) {
1495  case 0:
1496  break;
1497  case 1:
1498  xBatteryBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1499  break;
1500  case 2:
1501  yBatteryBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1502  break;
1503  case 3:
1504  zBatteryBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1505  break;
1506  default:
1507  break;
1508  }
1509 
1510  }
1511  tokens->Delete();
1512  }
1513  for(Int_t i=0;i<2;i++) {
1514  line.ReadLine(AntennaFile);
1515  //std::cout << line.Data() << "\n";
1516  }
1517  for(Int_t corner=0;corner<4;corner++) {
1518  line.ReadLine(AntennaFile);
1519  //std::cout << "Sip Box:\t" << line.Data() << "\n";
1520  TObjArray *tokens = line.Tokenize(",");
1521  for(Int_t j=0;j<4;j++) {
1522  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1523  // TString *subString = (TString*) tokens->At(j);
1524  // //std::cout << j << "\t" << subString.Data() << "\n";
1525  switch(j) {
1526  case 0:
1527  break;
1528  case 1:
1529  xSipBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1530  break;
1531  case 2:
1532  ySipBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1533  break;
1534  case 3:
1535  zSipBoxFromVerticalCorner[corner]=subString.Atof()*INCHTOMETER; //m
1536  break;
1537  default:
1538  break;
1539  }
1540 
1541  }
1542  tokens->Delete();
1543  }
1544  //Now gps stuff
1545  for(Int_t i=0;i<2;i++) {
1546  line.ReadLine(AntennaFile);
1547  //std::cout << line.Data() << "\n";
1548  }
1549  {
1550  line.ReadLine(AntennaFile);
1551  //std::cout << "GPS Plane:\t" << line.Data() << "\n";
1552  TObjArray *tokens = line.Tokenize(",");
1553  for(Int_t j=0;j<4;j++) {
1554  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1555  // TString *subString = (TString*) tokens->At(j);
1556  // //std::cout << j << "\t" << subString.Data() << "\n";
1557  switch(j) {
1558  case 0:
1559  break;
1560  case 1:
1561  case 2:
1562  case 3:
1563  gpsPlaneFromVertical[j-1]=subString.Atof(); //unit vector
1564  break;
1565  default:
1566  break;
1567  }
1568 
1569  }
1570  tokens->Delete();
1571  }
1572  {
1573  line.ReadLine(AntennaFile);
1574  //std::cout << "GPS Heading:\t" << line.Data() << "\n";
1575  TObjArray *tokens = line.Tokenize(",");
1576  for(Int_t j=0;j<4;j++) {
1577  const TString subString = ((TObjString*)tokens->At(j))->GetString();
1578  // TString *subString = (TString*) tokens->At(j);
1579  // //std::cout << j << "\t" << subString.Data() << "\n";
1580  switch(j) {
1581  case 0:
1582  break;
1583  case 1:
1584  case 2:
1585  case 3:
1586  gpsHeadingFromVertical[j-1]=subString.Atof(); // unit vector
1587  break;
1588  default:
1589  break;
1590  }
1591 
1592  }
1593  tokens->Delete();
1594  }
1595  aftForeOffsetAngleVertical=TMath::ATan(gpsHeadingFromVertical[1]/gpsHeadingFromVertical[0]);
1596  // fHeadingRotationAxis.SetXYZ(gpsPlaneFromVertical[0],gpsPlaneFromVertical[1],gpsPlaneFromVertical[2]);
1597  // fRollRotationAxis.SetXYZ(gpsHeadingFromVertical[0],gpsHeadingFromVertical[1],gpsHeadingFromVertical[2]);
1598 
1599  //Fix the heading and other axes to be 'ideal'
1600 
1601  //Now just for the sake of confusion we are going to redefine our definition
1602  //of phi equals zero to lie in the ADU5 fore direction.
1603  fHeadingRotationAxis.SetXYZ(0.,0.,1.);
1604  fPitchRotationAxis.SetXYZ(0.,1.,0.);
1605  fRollRotationAxis=fPitchRotationAxis.Cross(fHeadingRotationAxis);
1606 
1607  //ryans original pitch axis
1608  //fPitchRotationAxis=fRollRotationAxis.Cross(fHeadingRotationAxis);
1609  //reversed pitch axis
1610  //fPitchRotationAxis=fHeadingRotationAxis.Cross(fRollRotationAxis);
1611 
1612  //fPitchRotationAxis.SetXYZ(0.,1.,0.);
1613  //fRollRotationAxis=fPitchRotationAxis.Cross(fHeadingRotationAxis);
1614 
1615  //fitted gps axes
1616 // fPitchRotationAxis.SetXYZ(0,cos(-0.176*TMath::DegToRad()),sin(-0.176*TMath::DegToRad()));
1617 // fRollRotationAxis.SetXYZ(cos(0.460*TMath::DegToRad()),0,sin(0.460*TMath::DegToRad()));
1618 // fHeadingRotationAxis=fRollRotationAxis.Cross(fPitchRotationAxis);
1619 
1620  //2nd attempt fitted gps axes
1621 // fPitchRotationAxis.SetXYZ(0,cos(1.3*TMath::DegToRad()),sin(1.3*TMath::DegToRad()));
1622 // fRollRotationAxis.SetXYZ(cos(0.36*TMath::DegToRad()),0,sin(0.36*TMath::DegToRad()));
1623 // fHeadingRotationAxis=fRollRotationAxis.Cross(fPitchRotationAxis)
1624  ;
1625 
1626 
1627  //Now add in Simon's corrections
1628  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
1629  for(Int_t ant=0;ant<NUM_SEAVEYS;ant++) {
1630  rPhaseCentreFromVerticalHorn[ant][pol]+=deltaRPhaseCentre[ant][pol];
1631  // std::cout << rPhaseCentreFromVerticalHorn[ant][pol] << "\t" << deltaRPhaseCentre[ant][pol] << "\n";
1632  azPhaseCentreFromVerticalHorn[ant][pol]+=deltaPhiPhaseCentre[ant][pol];
1633  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
1634  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
1635  if(azPhaseCentreFromVerticalHorn[ant][pol]>TMath::TwoPi())
1636  azPhaseCentreFromVerticalHorn[ant][pol]-=TMath::TwoPi();
1637  zPhaseCentreFromVerticalHorn[ant][pol]+=deltaZPhaseCentre[ant][pol];
1638  xPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHorn[ant][pol]*TMath::Cos(azPhaseCentreFromVerticalHorn[ant][pol]);
1639  yPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHorn[ant][pol]*TMath::Sin(azPhaseCentreFromVerticalHorn[ant][pol]);
1640  }
1641  }
1642 
1643 // std::cout << " heading axis x " << fHeadingRotationAxis.x() << " y " << fHeadingRotationAxis.y() << " z " << fHeadingRotationAxis.z() << std::endl;
1644 // std::cout << " roll axis x " << fRollRotationAxis.x() << " y " << fRollRotationAxis.y() << " z " << fRollRotationAxis.z() << std::endl;
1645 // std::cout << " pitch axis x " << fPitchRotationAxis.x() << " y " << fPitchRotationAxis.y() << " z " << fPitchRotationAxis.z() << std::endl;
1646 
1647 
1648 }
1649 
1650 
1651 //Non static thingies
1652 void AnitaGeomTool::getAntXYZ(Int_t ant, Double_t &x, Double_t &y, Double_t &z,AnitaPol::AnitaPol_t pol)
1653 {
1654  if(ant>=0 && ant<NUM_SEAVEYS) {
1655  x=xPhaseCentreFromVerticalHorn[ant][pol];
1656  y=yPhaseCentreFromVerticalHorn[ant][pol];
1657  z=zPhaseCentreFromVerticalHorn[ant][pol];
1658  //std::cout << "ant " << x << " y " << y << " z " << std::endl;
1659  }
1660 }
1661 
1663  if(fUsePhotogrammetryNumbers==0) {
1664  if(ant>=0 && ant<NUM_SEAVEYS) {
1665  // std::cout << "Using simon z\n";
1666  return zPhaseCentreFromVerticalHorn[ant][pol];
1667  }
1668  }
1669  else {
1670  if(ant>=0 && ant<NUM_SEAVEYS) {
1671  // std::cout << "Using kurt z\n";
1672  return zPhaseCentreFromVerticalHornPhotogrammetry[ant][pol];
1673  }
1674  }
1675  return 0;
1676 }
1677 
1679  if(fUsePhotogrammetryNumbers==0) {
1680  // std::cout << "Using simon R\n";
1681  if(ant>=0 && ant<NUM_SEAVEYS) {
1682  return rPhaseCentreFromVerticalHorn[ant][pol];
1683  }
1684  }
1685  else {
1686  // std::cout << "Using kurt R\n";
1687  if(ant>=0 && ant<NUM_SEAVEYS) {
1688  return rPhaseCentreFromVerticalHornPhotogrammetry[ant][pol];
1689  }
1690  }
1691  return 0;
1692 }
1693 
1695 
1696  if(fUsePhotogrammetryNumbers==0) {
1697  // std::cout << "Using simon phi\n";
1698  if(ant>=0 && ant<NUM_SEAVEYS) {
1699  return azPhaseCentreFromVerticalHorn[ant][pol];
1700  }
1701  }
1702  else {
1703  if(ant>=0 && ant<NUM_SEAVEYS) {
1704  // std::cout << "Using kurt phi\n";
1705  return azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol];
1706  }
1707  }
1708  return 0;
1709 }
1710 
1712 
1713  if(fUsePhotogrammetryNumbers==0) {
1714  // std::cout << "Using simon phi-rel\n";
1715  if(ant>=0 && ant<NUM_SEAVEYS) {
1716  Double_t phi=azPhaseCentreFromVerticalHorn[ant][pol]-aftForeOffsetAngleVertical;
1717  if(phi<0)
1718  phi+=TMath::TwoPi();
1719  if(phi>TMath::TwoPi())
1720  phi-=TMath::TwoPi();
1721  return phi;
1722  }
1723  }
1724  else {
1725  if(ant>=0 && ant<NUM_SEAVEYS) {
1726  // std::cout << "Using kurt phi-rel\n";
1727  Double_t phi=azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]-aftForeOffsetAngleVerticalPhotogrammetry;
1728  // std::cout << ant << " " << phi*TMath::RadToDeg() << std::endl;
1729 
1730  if(phi<0)
1731  phi+=TMath::TwoPi();
1732  if(phi>TMath::TwoPi())
1733  phi-=TMath::TwoPi();
1734  return phi;
1735  }
1736  }
1737  return 0;
1738 }
1739 
1740 Double_t AnitaGeomTool::getMeanAntPairPhiRelToAftFore(Int_t firstAnt, Int_t secondAnt, AnitaPol::AnitaPol_t pol) {
1741 
1742  if(firstAnt>=0 && firstAnt<NUM_SEAVEYS && secondAnt>=0 && secondAnt<NUM_SEAVEYS) {
1743  //Calculating the phi of each antenna pair
1744  Double_t meanPhi=this->getAntPhiPositionRelToAftFore(firstAnt,pol);
1745  if(TMath::Abs(meanPhi-this->getAntPhiPositionRelToAftFore(secondAnt,pol))<TMath::PiOver2()) {
1746  meanPhi+=this->getAntPhiPositionRelToAftFore(secondAnt,pol);
1747  meanPhi*=0.5;
1748  }
1749  else if(TMath::Abs(meanPhi-(this->getAntPhiPositionRelToAftFore(secondAnt,pol)+TMath::TwoPi()))<TMath::PiOver2()) {
1750  meanPhi+=(this->getAntPhiPositionRelToAftFore(secondAnt,pol)+TMath::TwoPi());
1751  meanPhi*=0.5;
1752  }
1753  else if(TMath::Abs(meanPhi-(this->getAntPhiPositionRelToAftFore(secondAnt,pol)-TMath::TwoPi()))<TMath::PiOver2()) {
1754  meanPhi+=(this->getAntPhiPositionRelToAftFore(secondAnt,pol)-TMath::TwoPi());
1755  meanPhi*=0.5;
1756  }
1757  return meanPhi;
1758  }
1759  return -999;
1760 }
1761 
1762 
1763 
1765  if(phiWave<0) phiWave+=TMath::TwoPi();
1766  if(phiWave>TMath::TwoPi()) phiWave-=TMath::TwoPi();
1767  Double_t minDiff=TMath::TwoPi();
1768  Int_t minAnt=0;;
1769  for(Int_t ant=0;ant<16;ant++) {
1770  // std::cout << ant << "\t" << azPhaseCentreFromVerticalHorn[ant][pol] << "\t" << phiPrime << "\t" << TMath::Abs(azPhaseCentreFromVerticalHorn[ant][pol]-phiPrime) << "\n";
1771  if(TMath::Abs(getAntPhiPositionRelToAftFore(ant,pol)-phiWave)<minDiff) {
1772  minDiff=TMath::Abs(getAntPhiPositionRelToAftFore(ant,pol)-phiWave);
1773  minAnt=ant;
1774  }
1775 
1776  if(TMath::Abs((getAntPhiPositionRelToAftFore(ant,pol)-TMath::TwoPi())-phiWave)<minDiff) {
1777  minDiff=TMath::Abs((getAntPhiPositionRelToAftFore(ant,pol)-TMath::TwoPi())-phiWave);
1778  minAnt=ant;
1779  }
1780 
1781  if(TMath::Abs((getAntPhiPositionRelToAftFore(ant,pol)+TMath::TwoPi())-phiWave)<minDiff) {
1782  minDiff=TMath::Abs((getAntPhiPositionRelToAftFore(ant,pol)+TMath::TwoPi())-phiWave);
1783  minAnt=ant;
1784  }
1785 
1786  }
1787  return minAnt;
1788 }
1789 
1790 Int_t AnitaGeomTool::getUpperAntNearestPhiWave(Double_t phiWave, AnitaPol::AnitaPol_t pol) {
1791  return getTopAntNearestPhiWave(phiWave,pol);
1792 }
1793 
1794 
1795 //Non static thingies
1796 void AnitaGeomTool::getAntFaceXYZ(Int_t ant, Double_t &x, Double_t &y, Double_t &z)
1797 {
1798  if(ant>=0 && ant<NUM_SEAVEYS) {
1799  x=xAntFromVerticalHorn[ant];
1800  y=yAntFromVerticalHorn[ant];
1801  z=zAntFromVerticalHorn[ant];
1802  }
1803 }
1804 
1805 Double_t AnitaGeomTool::getAntFaceZ(Int_t ant) {
1806  if(ant>=0 && ant<NUM_SEAVEYS) {
1807  return zAntFromVerticalHorn[ant];
1808  // return zAntFaceFromDeckHorn[ant];
1809  }
1810  return 0;
1811 }
1812 
1813 Double_t AnitaGeomTool::getAntFaceR(Int_t ant) {
1814  if(ant>=0 && ant<NUM_SEAVEYS) {
1815  return rAntFromVerticalHorn[ant];
1816  //return rAntFromDeckHorn[ant];
1817  }
1818  return 0;
1819 }
1820 
1822  if(ant>=0 && ant<NUM_SEAVEYS) {
1823  return azCentreFromVerticalHorn[ant];
1824  }
1825  return 0;
1826 }
1827 
1829  if(ant>=0 && ant<NUM_SEAVEYS) {
1830  Double_t phi=azCentreFromVerticalHorn[ant]-aftForeOffsetAngleVertical;
1831  if(phi<0)
1832  phi+=TMath::TwoPi();
1833  return phi;
1834  }
1835  return 0;
1836 }
1837 
1839  Double_t phiPrime=phiWave+aftForeOffsetAngleVertical;
1840  if(phiPrime>TMath::TwoPi())
1841  phiPrime-=TMath::TwoPi();
1842  Double_t minDiff=TMath::TwoPi();
1843  Int_t minAnt=0;;
1844  for(Int_t ant=0;ant<16;ant++) {
1845  // std::cout << ant << "\t" << azPhaseCentreFromVerticalHorn[ant][pol] << "\t" << phiPrime << "\t" << TMath::Abs(azPhaseCentreFromVerticalHorn[ant][pol]-phiPrime) << "\n";
1846  if(TMath::Abs(azCentreFromVerticalHorn[ant]-phiPrime)<minDiff) {
1847  minDiff=TMath::Abs(azCentreFromVerticalHorn[ant]-phiPrime);
1848  minAnt=ant;
1849  }
1850  }
1851  return minAnt;
1852 }
1853 
1854 void AnitaGeomTool::updateAnt(Double_t deltaR,Double_t deltaRL,Double_t deltaUD){
1855 
1856  Double_t phaseCentreToAntFront=0.2+deltaR;
1857 
1858  Double_t deltaXRL = 0;
1859  Double_t deltaYRL = 0;
1860 
1861  Double_t deltaZUD = 0;
1862  Double_t deltaXUD = 0;
1863  Double_t deltaYUD = 0;
1864 
1865  std::cout << deltaXRL << " " << deltaXUD << " "<< deltaYUD << " "<< deltaYRL << " " << deltaZUD << " " << deltaRL << std::endl;
1866 
1867  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
1868  for(Int_t ant=0;ant<32;ant++) {
1869 
1870  //deltaXRL = -deltaRL*TMath::Sin(apertureAzFromVerticalHorn[ant]);
1871  // deltaYRL = deltaRL*TMath::Cos(apertureAzFromVerticalHorn[ant]);
1872 
1873  deltaZUD = deltaUD*TMath::Cos(apertureElFromVerticalHorn[ant]);
1874  deltaXUD = -deltaUD*TMath::Sin(apertureElFromVerticalHorn[ant])*TMath::Cos(apertureAzFromVerticalHorn[ant]);
1875  deltaYUD = -deltaUD*TMath::Sin(apertureElFromVerticalHorn[ant])*TMath::Sin(apertureAzFromVerticalHorn[ant]);
1876 
1877  //std::cout << deltaXRL << " " << deltaXUD << " "<< deltaYUD << " "<< deltaYRL << " " << deltaZUD << std::endl;
1878 
1879  //Now try to find the location of the antenna phaseCentre poInt_t
1880  xPhaseCentreFromVerticalHorn[ant][pol]=(xAntFromVerticalHorn[ant] + deltaXRL + deltaXUD) -
1881  phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1882  yPhaseCentreFromVerticalHorn[ant][pol]=(yAntFromVerticalHorn[ant] + deltaYRL + deltaYUD) -
1883  phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1884  zPhaseCentreFromVerticalHorn[ant][pol]=(zAntFromVerticalHorn[ant] + deltaZUD) -phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1885 
1886 
1887  rPhaseCentreFromVerticalHorn[ant][pol]=TMath::Sqrt(xPhaseCentreFromVerticalHorn[ant][pol]*xPhaseCentreFromVerticalHorn[ant][pol]+
1888  yPhaseCentreFromVerticalHorn[ant][pol]*yPhaseCentreFromVerticalHorn[ant][pol]);//m
1889  azPhaseCentreFromVerticalHorn[ant][pol]=TMath::ATan(yPhaseCentreFromVerticalHorn[ant][pol]/xPhaseCentreFromVerticalHorn[ant][pol]); //radians
1890  if(xPhaseCentreFromVerticalHorn[ant][pol]<0)
1891  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::Pi();
1892  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
1893  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
1894 
1895  }
1896 }
1897 
1898 
1899  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
1900  for(Int_t ant=32;ant<48;ant++) {
1901 
1902  Double_t deltaXRL = 0;
1903  Double_t deltaYRL = 0;
1904 
1905  Double_t deltaZUD = 0;
1906  Double_t deltaXUD = 0;
1907  Double_t deltaYUD = 0;
1908 
1909  deltaZUD = deltaUD*TMath::Cos(apertureElFromVerticalHorn[ant]);
1910  deltaXUD = -deltaUD*TMath::Sin(apertureElFromVerticalHorn[ant])*TMath::Cos(apertureAzFromVerticalHorn[ant]);
1911  deltaYUD = -deltaUD*TMath::Sin(apertureElFromVerticalHorn[ant])*TMath::Sin(apertureAzFromVerticalHorn[ant]);
1912 
1913  // rAntFromVerticalHorn[ant]=(101*INCHTOMETER)-0.38;
1914  rAntFromVerticalHorn[ant]=(101*INCHTOMETER)-0.3;
1915  zAntFromVerticalHorn[ant]=(-21*INCHTOMETER)-1.8;
1916  azCentreFromVerticalHorn[ant]=(-67.5 + 45*(ant-32))*TMath::DegToRad();
1917  apertureAzFromVerticalHorn[ant]=azCentreFromVerticalHorn[ant];
1918  apertureElFromVerticalHorn[ant]=-10*TMath::DegToRad();
1919  xAntFromVerticalHorn[ant]=rAntFromVerticalHorn[ant]*TMath::Cos(azCentreFromVerticalHorn[ant]);
1920  yAntFromVerticalHorn[ant]=rAntFromVerticalHorn[ant]*TMath::Sin(azCentreFromVerticalHorn[ant]);
1921 
1922 
1923  xPhaseCentreFromVerticalHorn[ant][pol]=(xAntFromVerticalHorn[ant] + deltaXRL + deltaXUD) -
1924  phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1925  yPhaseCentreFromVerticalHorn[ant][pol]=(yAntFromVerticalHorn[ant] + deltaYRL + deltaYUD) -
1926  phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1927  zPhaseCentreFromVerticalHorn[ant][pol]=(zAntFromVerticalHorn[ant] + deltaZUD) -phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1928 
1929  // //Now try to find the location of the antenna phaseCentre poInt_t
1930  // xPhaseCentreFromVerticalHorn[ant]=xAntFromVerticalHorn[ant]-
1931  // phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1932  // yPhaseCentreFromVerticalHorn[ant]=yAntFromVerticalHorn[ant]-
1933  // phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHorn[ant])*TMath::Cos(apertureElFromVerticalHorn[ant]); //m
1934  // zPhaseCentreFromVerticalHorn[ant]=zAntFromVerticalHorn[ant]-phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHorn[ant]); //m
1935  rPhaseCentreFromVerticalHorn[ant][pol]=TMath::Sqrt(xPhaseCentreFromVerticalHorn[ant][pol]*xPhaseCentreFromVerticalHorn[ant][pol]+
1936  yPhaseCentreFromVerticalHorn[ant][pol]*yPhaseCentreFromVerticalHorn[ant][pol]);//m
1937  azPhaseCentreFromVerticalHorn[ant][pol]=TMath::ATan(yPhaseCentreFromVerticalHorn[ant][pol]/xPhaseCentreFromVerticalHorn[ant][pol]); //radians
1938  if(xPhaseCentreFromVerticalHorn[ant][pol]<0)
1939  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::Pi();
1940  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
1941  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
1942 
1943  // std::cout << "Simon's:\t" << ant << " "
1944 // << rPhaseCentreFromVerticalHorn[ant][pol] << " "
1945 // << azPhaseCentreFromVerticalHorn[ant][pol] << " "
1946 // << zPhaseCentreFromVerticalHorn[ant][pol] << "\n";
1947  }
1948 }
1949 
1950 
1951 
1952 };
1953 
1954 
1955 
1956 void AnitaGeomTool::printAntPos(){
1957 
1958  for(Int_t pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
1959  for(Int_t ant = 0; ant < 32; ant++){
1960  std::cout << rPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
1961  }
1962  }
1963 }
1964 
1965 
1966 
1968 {
1969  if(ant<16)
1970  return AnitaGeom::topPhiNums[ant];
1971  else if(ant<32)
1972  return AnitaGeom::middlePhiNums[ant-16];
1973  else if(ant<48)
1974  return AnitaGeom::bottomPhiNums[ant-32];
1975  std::cerr << "There isn't an antenna " << ant << " (0-47 only)\n";
1976  return -1;
1977 
1978 }
1979 
1980 
1982 {
1983  switch(ring) {
1984  case AnitaRing::kTopRing:
1985  return AnitaGeom::topAntNums[phi];
1987  return AnitaGeom::middleAntNums[phi];
1989  return AnitaGeom::bottomAntNums[phi];
1990  default:
1991  return -1;
1992  }
1993 
1994 }
1995 
1996 void AnitaGeomTool::readSimonsNumbers() {
1997 
1998  char calibDir[FILENAME_MAX];
1999  char fileName[FILENAME_MAX];
2000  char *calibEnv=getenv("ANITA_CALIB_DIR");
2001  if(!calibEnv) {
2002  char *utilEnv=getenv("ANITA_UTIL_INSTALL_DIR");
2003  if(!utilEnv)
2004  sprintf(calibDir,"calib");
2005  else
2006  sprintf(calibDir,"%s/share/anitaCalib",utilEnv);
2007  }
2008  else {
2009  strncpy(calibDir,calibEnv,FILENAME_MAX);
2010  }
2011 
2012  sprintf(fileName,"%s/simonsPositionAndTimingOffsets.dat",calibDir);
2013  std::ifstream SimonFile(fileName);
2014  if(!SimonFile) {
2015  std::cerr << "Couldn't open:\t" << fileName << "\n";
2016  return;
2017  }
2018 
2019  Int_t antNum;
2020  Double_t deltaT,deltaR,deltaPhi,deltaZ;
2021  char firstLine[180];
2022 
2023  SimonFile.getline(firstLine,179);
2024  while(SimonFile >> antNum >> deltaT >> deltaR >> deltaPhi >> deltaZ) {
2025  deltaRPhaseCentre[antNum][AnitaPol::kVertical]=deltaR;
2026  deltaPhiPhaseCentre[antNum][AnitaPol::kVertical]=deltaPhi;
2027  deltaZPhaseCentre[antNum][AnitaPol::kVertical]=deltaZ;
2028  }
2029 
2030 
2031  //The H-pol file is (for now at least) just a delta file relative to the
2032  // V-pol
2033  sprintf(fileName,"%s/simonsHPolPositionAndTimingOffsets.dat",calibDir);
2034  std::ifstream SimonHPolFile(fileName);
2035  if(!SimonHPolFile) {
2036  std::cerr << "Couldn't open:\t" << fileName << "\n";
2037  return;
2038  }
2039 
2040  SimonHPolFile.getline(firstLine,179);
2041  while(SimonHPolFile >> antNum >> deltaT >> deltaR >> deltaPhi >> deltaZ) {
2042  deltaRPhaseCentre[antNum][AnitaPol::kHorizontal]=
2043  deltaRPhaseCentre[antNum][AnitaPol::kVertical]+deltaR;
2044  deltaPhiPhaseCentre[antNum][AnitaPol::kHorizontal]=
2045  deltaPhiPhaseCentre[antNum][AnitaPol::kVertical]+deltaPhi;
2046  deltaZPhaseCentre[antNum][AnitaPol::kHorizontal]=
2047  deltaZPhaseCentre[antNum][AnitaPol::kVertical]+deltaZ;
2048  }
2049 
2050 
2051 }
2052 
2053 
2054 
2055 void AnitaGeomTool::readPhaseCenterNumbers(int version) {
2056 
2057  char calibDir[FILENAME_MAX];
2058  char fileName[FILENAME_MAX];
2059  char *calibEnv=getenv("ANITA_CALIB_DIR");
2060  if(!calibEnv) {
2061  char *utilEnv=getenv("ANITA_UTIL_INSTALL_DIR");
2062  if(!utilEnv)
2063  sprintf(calibDir,"calib");
2064  else
2065  sprintf(calibDir,"%s/share/anitaCalib",utilEnv);
2066  }
2067  else {
2068  strncpy(calibDir,calibEnv,FILENAME_MAX);
2069  }
2070 
2071  switch(version) {
2072  case 3 : sprintf(fileName,"%s/phaseCenterPositionsRelativeToPhotogrammetryAnita3.dat",calibDir);
2073  break; // and exits the switch
2074  case 4 : sprintf(fileName,"%s/phaseCenterPositionsRelativeToPhotogrammetryAnita4.dat",calibDir);
2075  break;
2076  default :
2077  std::cerr << "Unknown AnitaVersion " << version << "\n";
2078  return;
2079  }
2080 
2081 
2082  std::ifstream PhaseCenterFile(fileName);
2083  if(!PhaseCenterFile) {
2084  std::cerr << "Couldn't open:\t" << fileName << "\n";
2085  return;
2086  }
2087 
2088  Int_t antNum, pol;
2089  Double_t deltaR,deltaPhi,deltaZ;
2090  char firstLine[180];
2091 
2092  PhaseCenterFile.getline(firstLine,179);
2093  while(PhaseCenterFile >> antNum >> pol >> deltaR >> deltaPhi >> deltaZ) {
2094  deltaRPhaseCentre[antNum][pol]=deltaR;
2095  deltaPhiPhaseCentre[antNum][pol]=deltaPhi*TMath::DegToRad();
2096  deltaZPhaseCentre[antNum][pol]=deltaZ;
2097  }
2098 
2099 }
2100 
2101 
2102 
2103 
2104 
2105 
2106 
2107 
2108 void AnitaGeomTool::readPhotogrammetry(int version)
2109 {
2110  // return; //Hack for now
2111  char calibDir[FILENAME_MAX];
2112  char fileName[FILENAME_MAX];
2113  char *calibEnv=getenv("ANITA_CALIB_DIR");
2114  if(!calibEnv) {
2115  char *utilEnv=getenv("ANITA_UTIL_INSTALL_DIR");
2116  if(!utilEnv)
2117  sprintf(calibDir,"calib");
2118  else
2119  sprintf(calibDir,"%s/share/anitaCalib",utilEnv);
2120  }
2121  else {
2122  strncpy(calibDir,calibEnv,FILENAME_MAX);
2123  }
2124 
2125 
2126  switch(version) {
2127  case 3 : sprintf(fileName,"%s/anitaIIIPhotogrammetry.csv",calibDir);
2128  break; // and exits the switch
2129  case 4 : sprintf(fileName,"%s/anitaIVPhotogrammetry.csv",calibDir);
2130  break;
2131  default :
2132  std::cerr << "Unknown AnitaVersion " << version << "\n";
2133  return;
2134  }
2135 
2136  std::ifstream AnitaPhotoFile(fileName);
2137  if(!AnitaPhotoFile) {
2138  std::cerr << "Couldn't open:\t" << fileName << "\n";
2139  return;
2140  }
2141 
2142 
2143  //First up are the antenna positions
2144  TString line;
2145  for(int i=0;i<2;i++) {
2146  line.ReadLine(AnitaPhotoFile);
2147  //std::cout << line.Data() << "\n";
2148  }
2149 
2150  for(int ant=0;ant<48;ant++) {
2151  line.ReadLine(AnitaPhotoFile);
2152  //std::cout << "Seavey:\t" << line.Data() << "\n";
2153  TObjArray *tokens = line.Tokenize(",");
2154  for(int j=0;j<8;j++) {
2155  const TString subString = ((TObjString*)tokens->At(j))->GetString();
2156  // TString *subString = (TString*) tokens->At(j);
2157  // std::cout << j << "\t" << subString.Data() << "\n";
2158  switch(j) {
2159  case 0:
2160  if (ant+1 != subString.Atoi()) {
2161  std::cerr << "Antenna number mismatch\n";
2162  }
2163  break;
2164  case 1:
2165  xAntFromVerticalHornPhotogrammetry[ant]=subString.Atof()*INCHTOMETER; //m
2166  break;
2167  case 2:
2168  yAntFromVerticalHornPhotogrammetry[ant]=subString.Atof()*INCHTOMETER; //m
2169  break;
2170  case 3:
2171  zAntFromVerticalHornPhotogrammetry[ant]=subString.Atof()*INCHTOMETER; //m
2172  break;
2173  case 4:
2174  rAntFromVerticalHornPhotogrammetry[ant]=subString.Atof()*INCHTOMETER; //m
2175  // std::cout << ant << "\t" << rAntFromVerticalHornKurtAnitaII[ant] << "\n";
2176  break;
2177  case 5:
2178  azCentreFromVerticalHornPhotogrammetry[ant]=subString.Atof()*TMath::DegToRad(); //radians
2179  break;
2180  case 6:
2181  apertureAzFromVerticalHornPhotogrammetry[ant]=subString.Atof()*TMath::DegToRad(); //radians
2182  break;
2183  case 7:
2184  apertureElFromVerticalHornPhotogrammetry[ant]=subString.Atof()*TMath::DegToRad(); //radians
2185  break;
2186  default:
2187  break;
2188  }
2189 
2190  }
2191  tokens->Delete();
2192 
2193  }
2194 
2195  for(int pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
2196  for(int ant=0;ant<NUM_SEAVEYS;ant++) {
2197  Double_t phaseCentreToAntFront=ringPhaseCentreOffset[Int_t(this->getRingFromAnt(ant))];
2198  //Need to think about this at some point, but for now will just hard code in 20cm
2199  phaseCentreToAntFront=0.2;
2200 
2201 
2202  // Now try to find the location of the antenna phaseCentre point
2203  xPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]=xAntFromVerticalHornPhotogrammetry[ant] -
2204  phaseCentreToAntFront*TMath::Cos( apertureAzFromVerticalHornPhotogrammetry[ant])*TMath::Cos(apertureElFromVerticalHornPhotogrammetry[ant]); //m
2205  yPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]=yAntFromVerticalHornPhotogrammetry[ant] -
2206  phaseCentreToAntFront*TMath::Sin( apertureAzFromVerticalHornPhotogrammetry[ant])*TMath::Cos(apertureElFromVerticalHornPhotogrammetry[ant]); //m
2207  zPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]=zAntFromVerticalHornPhotogrammetry[ant] - phaseCentreToAntFront*TMath::Sin(apertureElFromVerticalHornPhotogrammetry[ant]); //m
2208 
2209  rPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]=TMath::Sqrt(xPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]*xPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]+yPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]*yPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]);//m
2210  azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]=TMath::ATan(yPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]/xPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]); //radians
2211  if(xPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]<0)
2212  azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]+=TMath::Pi();
2213  if(azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]<0)
2214  azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]+=TMath::TwoPi();
2215 
2216 
2217  }
2218  }
2219 
2220  aftForeOffsetAngleVerticalPhotogrammetry= +45*TMath::DegToRad();
2221  aftForeOffsetAngleVertical=aftForeOffsetAngleVerticalPhotogrammetry;
2222  //Will now make the default numbers the ones that start with Simon's and modify the positions
2223 
2224  //Removed corrections for the moment
2225  for(int pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
2226  for(int ant=0;ant<NUM_SEAVEYS;ant++) {
2227  rPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHornPhotogrammetry[ant][pol];
2228  azPhaseCentreFromVerticalHorn[ant][pol]=azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol];
2229  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
2230  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
2231  if(azPhaseCentreFromVerticalHorn[ant][pol]>TMath::TwoPi())
2232  azPhaseCentreFromVerticalHorn[ant][pol]-=TMath::TwoPi();
2233  zPhaseCentreFromVerticalHorn[ant][pol]=zPhaseCentreFromVerticalHornPhotogrammetry[ant][pol];
2234  xPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHorn[ant][pol]*TMath::Cos(azPhaseCentreFromVerticalHorn[ant][pol]);
2235  yPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHorn[ant][pol]*TMath::Sin(azPhaseCentreFromVerticalHorn[ant][pol]);
2236  // std::cout << xPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
2237 
2238  }
2239  }
2240 
2241  if (fUsePhotogrammetryNumbers==0){
2242  addPhaseCenters();
2243  }
2244 
2245 
2246 }
2247 
2248 
2249 void AnitaGeomTool::addPhaseCenters(){
2250 
2251 
2252  // Now add in ANITA-3 phase centre corrections
2253  // std::cerr << "Inside AnitaGeomTool we have..." << std::endl;
2254  // std::cerr << "pol\tant\tr\tz\tphi" << std::endl;
2255  for(int pol=AnitaPol::kHorizontal;pol<=AnitaPol::kVertical;pol++) {
2256  for(int ant=0;ant<NUM_SEAVEYS;ant++) {
2257  rPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]+deltaRPhaseCentre[ant][pol];
2258  azPhaseCentreFromVerticalHorn[ant][pol]=azPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]+deltaPhiPhaseCentre[ant][pol];
2259  if(azPhaseCentreFromVerticalHorn[ant][pol]<0)
2260  azPhaseCentreFromVerticalHorn[ant][pol]+=TMath::TwoPi();
2261  if(azPhaseCentreFromVerticalHorn[ant][pol]>TMath::TwoPi())
2262  azPhaseCentreFromVerticalHorn[ant][pol]-=TMath::TwoPi();
2263  zPhaseCentreFromVerticalHorn[ant][pol]=zPhaseCentreFromVerticalHornPhotogrammetry[ant][pol]+deltaZPhaseCentre[ant][pol];
2264  xPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHorn[ant][pol]*TMath::Cos(azPhaseCentreFromVerticalHorn[ant][pol]);
2265  yPhaseCentreFromVerticalHorn[ant][pol]=rPhaseCentreFromVerticalHorn[ant][pol]*TMath::Sin(azPhaseCentreFromVerticalHorn[ant][pol]);
2266  // std::cout << xPhaseCentreFromVerticalHorn[ant][pol] << std::endl;
2267 
2268 
2269  // std::cerr << pol << "\t" << ant << "\t" << rPhaseCentreFromVerticalHorn[ant][pol] << "\t" << zPhaseCentreFromVerticalHorn[ant][pol] << "\t" << getAntPhiPositionRelToAftFore(ant, (AnitaPol::AnitaPol_t)pol)*TMath::RadToDeg() << std::endl;
2270  // azPhaseCentreFromVerticalHorn[ant][pol]*TMath::RadToDeg() << std::endl;
2271 
2272  }
2273  }
2274 
2275 }
static Int_t getSurfL2TriggerChanFromPhi(Int_t phi, Int_t &surf, Int_t &l2Chan)
Converts phi-sector to surf and l2 trigger channel.
static Int_t getPhiRingFromSurfL1Chan(Int_t surf, Int_t l1Chan, Int_t &phi, AnitaRing::AnitaRing_t &ring)
Converts the SURF and L1 trigger channel to phi-sector and polarization.
static Int_t getLayer(Int_t irx)
get layer given antenna number
static void getThetaPartners(Int_t rx, int &rxleft, int &rxright)
output the antennas that are in neighbouring phi sectors (not the neighbouring antennas) ...
Int_t getTopAntNearestPhiWave(Double_t phiWave, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
get antenna closest to given plane wave direction
Double_t getDirectionWrtNorth(Int_t phi, Double_t heading)
Get direction that a phi sector is pointing wrt north. Also takes heading as a input.
AnitaRing::AnitaRing_t surfTriggerChanToRing[SCALERS_PER_SURF]
< Map from SURF trigger channel to ring
Double_t getAntPhiPositionRelToAftFore(Int_t ant, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
get antenna phi position relative to ADU5 AFT-FORE line
Int_t surfToAntMapAnita3[ACTIVE_SURFS][RFCHAN_PER_SURF]
Map from phi-sector to antenna. Both start counting at zero.
static void getRingAntPolPhiFromSurfChan(Int_t surf, Int_t chan, AnitaRing::AnitaRing_t &ring, Int_t &ant, AnitaPol::AnitaPol_t &pol, Int_t &phi)
Convert surf-chan to ring-ant-pol-phi.
Double_t getAntZ(Int_t ant, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
get antenna z position
Int_t getTopAntFaceNearestPhiWave(Double_t phiWave)
get upper antenna closest to given plane wave direction
static void getNeighbors(Int_t rx, int &rxleft, int &rxright)
output the neighboring antennas
static Int_t getChanIndexFromAntPol(Int_t ant, AnitaPol::AnitaPol_t pol)
Convert ant-pol to logical index.
static Int_t getAntPolFromSurfChan(Int_t surf, Int_t chan, Int_t &ant, AnitaPol::AnitaPol_t &pol)
Convert surf-chan to ant-pol.
static Int_t getSurfChanFromChanIndex(Int_t chanIndex, Int_t &surf, Int_t &chan)
Convert logical index to surf and channel.
Left-circular polarisation (e.g. A4)
static Int_t getChanIndexFromRingPhiPol(AnitaRing::AnitaRing_t ring, Int_t phi, AnitaPol::AnitaPol_t pol)
Convert ring-phi-pol to logical index.
static Double_t getPhiDiff(Double_t firstPhi, Double_t secondPhi)
Returns the angular difference between two phi values (in radians)
Int_t vAntToChan[NUM_SEAVEYS]
Map for HPOL channel of antenna to channel on SURF. (HPOL channels are 4-7)
Horizontal Polarisation (e.g. A3)
USeful in for loops.
Double_t getAntR(Int_t ant, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
get antenna r position
static Int_t getPhiRingPolFromSurfChanTrigger(Int_t surf, Int_t chan, Int_t &phi, AnitaRing::AnitaRing_t &ring, AnitaTrigPol::AnitaTrigPol_t &pol)
Converts the SURF and channel numbers to phi-sector, ring and polarization.
static void getSurfChanAntFromRingPhiPol(AnitaRing::AnitaRing_t ring, Int_t phi, AnitaPol::AnitaPol_t pol, Int_t &surf, Int_t &chan, Int_t &ant)
Convert ring-phi-pol to surf-chan-ant.
static Int_t getAzimuthPartner(Int_t rx)
output the antenna that points to the same place in phi as the input antenna
void getAntFaceXYZ(Int_t ant, Double_t &x, Double_t &y, Double_t &z)
get location fo antenna face in balloon cartesian coordinates
Double_t getMeanAntPairPhiRelToAftFore(Int_t firstAnt, Int_t secondAnt, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
The mean of the two antenna phi positions.
Right-circular polarisation (e.g. A4)
Int_t bottomPhiNums[NUM_PHI]
Map from SURF and polarization to phi-sector.
Double_t getAntPhiPosition(Int_t ant, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
get antenna phi position
void getAntXYZ(Int_t ant, Double_t &x, Double_t &y, Double_t &z, AnitaPol::AnitaPol_t pol=AnitaPol::kVertical)
get antenna cartesian coordinates (from photogrammetry)
enum AnitaRing::EAnitaRing AnitaRing_t
Ring enumeration.
enum AnitaTrigPol::EAnitaTrigPol AnitaTrigPol_t
Polarisation enumeration.
AnitaTrigPol::AnitaTrigPol_t surfTriggerChanToPolAnita3[SCALERS_PER_SURF]
< Map from SURF trigger channel to polarization
AnitaGeom – Namespace used inside AnitaGeomTool.
Int_t hAntToChan[NUM_SEAVEYS]
1 is Normal orientation, -1 is 180 degree flip. (Top ring needs to be inverted in software when signa...
static Int_t getAntOrientation(Int_t ant)
Some of the antennas have their orientation reversed relative to nominal. The effect of this is to sw...
virtual ~AnitaGeomTool()
Default destructor.
static AnitaRing::AnitaRing_t getRingFromAnt(Int_t ant)
Get ring form antenna number (0-39)
static Int_t getSurfL1TriggerChanFromPhiRing(Int_t phi, AnitaRing::AnitaRing_t ring, Int_t &surf, Int_t &l1Chan)
Convert phi-sector and polarization from phi-sector and polarization.
Double_t getAntFaceR(Int_t ant)
get r position for antenna face
static Int_t getPhiPolFromSurfL1Chan(Int_t surf, Int_t l1Chan, Int_t &phi, AnitaPol::AnitaPol_t &pol)
Converts the SURF and L1 trigger channel to phi-sector and polarization.
Int_t antToSurfMap[NUM_SEAVEYS]
< Map from antenna to SURF. Both polarizations from an antenna go to the same SURF.
static Int_t getPhiFromAnt(Int_t ant)
get phi from ant
static Int_t getAntFromPhiRing(Int_t phi, AnitaRing::AnitaRing_t ring)
get antenna number from phi and ring
static Int_t getSurfChanTriggerFromPhiRingPol(Int_t phi, AnitaRing::AnitaRing_t ring, AnitaTrigPol::AnitaTrigPol_t pol, Int_t &surf, Int_t &chan)
Converts the SURF and L1 trigger channel to phi-sector and polarization.
Int_t bottomAntNums[NUM_PHI]
And the inverse, a map from antenna to phi-sector (using ANT-1 and ANT-17 and ANT-33 with the arrays)...
static Int_t getChanIndex(Int_t surf, Int_t chan)
Surf + channel to channel index.
Vertical Polarisation.
Useful in for loops.
Double_t getAntFacePhiPosition(Int_t ant)
get phi position for antenna face
Vertical Polarisation (e.g. A3)
static Int_t getPhiSector(Int_t rx)
phi sector of this antenna. rx runs from 0 to 31.
Int_t antOrientationMap[NUM_SEAVEYS]
Horizontal Polarisation.
Double_t getAntFaceZ(Int_t ant)
get z position for antenna face
enum AnitaPol::EAnitaPol AnitaPol_t
Polarisation enumeration.
static AnitaGeomTool * Instance(int anita_version=0)
Instance generator. If version_number == 0, uses AnitaVersion::get();.
AnitaGeomTool – The ANITA Geometry Tool.
Definition: AnitaGeomTool.h:48
Double_t getAntFacePhiPositionRelToAftFore(Int_t ant)
get phi position relative to ADU5 AFT-FORE direction
void fillAntPositionsFromPrioritizerdConfig()
Used for evaluating ANITA-3 Prioritizerd performance.