TH2DAntarctica.cxx
1 #include "TH2DAntarctica.h"
2 #include "RampdemReader.h"
3 #include "AntarcticaBackground.h"
4 #include "TGraphAntarctica.h"
5 #include "TRandom3.h"
6 #include "TVirtualPad.h"
7 #include "TPaletteAxis.h"
8 #include "TList.h"
9 #include "AntarcticaGeometry.h"
10 
11 ClassImp(TProfile2DAntarctica)
12 ClassImp(TH2DAntarctica)
13 
14 const double palX1 = AntarcticaBackgroundDefaults::zAxisRightMargin;
15 const double palX2 = AntarcticaBackgroundDefaults::zAxisRightMargin + AntarcticaBackgroundDefaults::zAxisWidth;
16 const double palY1 = 1 - AntarcticaBackgroundDefaults::zAxisTopBottomMargin - AntarcticaBackgroundDefaults::zAxisHeight;
17 const double palY2 = 1 - AntarcticaBackgroundDefaults::zAxisTopBottomMargin;
18 
19 TProfile2DAntarctica::TProfile2DAntarctica(Int_t nx, Int_t ny)
20  : TProfile2D(), fAntarcticaBackground(NULL){
21 
22  AntarcticaBackground* b = getBackground();
23 
24  nx = nx <= 0 ? b->GetNbinsX() : nx;
25  ny = ny <= 0 ? b->GetNbinsY() : ny;
26  SetBins(nx, b->GetXaxis()->GetBinLowEdge(1), b->GetXaxis()->GetBinUpEdge(b->GetNbinsX()),
27  ny, b->GetYaxis()->GetBinLowEdge(1), b->GetYaxis()->GetBinUpEdge(b->GetNbinsY()));
28  accept_stereographic = false;
29 }
30 
31 
32 TProfile2DAntarctica::TProfile2DAntarctica(const char* name, const char* title, Int_t nx, Int_t ny)
33  : TProfile2D(), fAntarcticaBackground(NULL){
34 
35  SetNameTitle(name, title);
36  AntarcticaBackground* b = getBackground();
37 
38  nx = nx <= 0 ? b->GetNbinsX() : nx;
39  ny = ny <= 0 ? b->GetNbinsY() : ny;
40  SetBins(nx, b->GetXaxis()->GetBinLowEdge(1), b->GetXaxis()->GetBinUpEdge(b->GetNbinsX()),
41  ny, b->GetYaxis()->GetBinLowEdge(1), b->GetYaxis()->GetBinUpEdge(b->GetNbinsY()));
42  accept_stereographic = false;
43 }
44 
45 
46 TProfile2DAntarctica::TProfile2DAntarctica(const char* name, const char* title, const std::vector<Double_t>& x, const std::vector<Double_t>& y)
47  : TProfile2D(), fAntarcticaBackground(NULL){
48 
49  SetNameTitle(name, title);
50  if(x.size() > 1 && y.size() > 1){
51  SetBins(x.size()-1, &x[0], y.size()-1, &y[0]);
52  }
53  else{
54  AntarcticaBackground* b = getBackground();
55 
56  Int_t nx = b->GetNbinsX();
57  Int_t ny = b->GetNbinsY();
58  SetBins(nx, b->GetXaxis()->GetBinLowEdge(1), b->GetXaxis()->GetBinUpEdge(b->GetNbinsX()),
59  ny, b->GetYaxis()->GetBinLowEdge(1), b->GetYaxis()->GetBinUpEdge(b->GetNbinsY()));
60 
61  }
62  accept_stereographic = false;
63 }
64 
65 
66 void TProfile2DAntarctica::Draw(Option_t* opt){
67 
68  TString opt2(opt);
69  if(!opt2.Contains("same")){
70  AntarcticaBackground* b = getBackground();
71  b->Draw();
72  b->SetBit(kCanDelete, false);
73  }
74  else{
75  TList* l = gPad->GetListOfPrimitives();
76  for(int i=0; i < l->GetEntries(); i++){
77  TString objName = l->At(i)->GetName();
78  if(objName.Contains("fAntarctica") && !objName.Contains("PalSetter")){
79  if(fAntarcticaBackground && fAntarcticaBackground != (AntarcticaBackground*) l->At(i)){
80  delete fAntarcticaBackground;
81  }
82  fAntarcticaBackground = (AntarcticaBackground*) l->At(i);
83  // break;
84  }
85  }
86  }
87  TString sameOpt = opt2.Contains("same") ? opt2 : TString::Format("%s same", opt);
88  TProfile2D::Draw(sameOpt);
89  ResetColorAxis();
90  RescaleBackground();
91  gPad->Modified();
92 }
93 
94 
95 Int_t TProfile2DAntarctica::Fill(Double_t lon, Double_t lat, Double_t val){
96 
97  Double_t easting, northing;
98  if (accept_stereographic)
99  {
100  easting = lon;
101  northing = lat;
102  }
103  else
104  {
105  RampdemReader::LonLatToEastingNorthing(lon, lat, easting, northing);
106  }
107  return TProfile2D::Fill(easting, northing, val);
108 }
109 
110 
116 void TProfile2DAntarctica::FillRandomly(Int_t nTimes){
117 
118  TRandom3 randy;
119 
120  for(int i=0; i < nTimes; i++){
121 
122  double lon = randy.Uniform(0, 360);
123  double lat = randy.Uniform(-90, -60);
124 
125  Fill(lon, lat);
126  }
127 }
128 
129 
130 
134 void TProfile2DAntarctica::ResetColorAxis(){
135  gPad->Modified();
136  gPad->Update();
137  TPaletteAxis *palette = (TPaletteAxis*) GetListOfFunctions()->FindObject("palette");
138  if(palette){
139  palette->SetX1NDC(palX1);
140  palette->SetX2NDC(palX2);
141  palette->SetY1NDC(palY1);
142  palette->SetY2NDC(palY2);
143 
144  TAxis* zAxis = GetZaxis();
145  // zAxis->SetTitle();
146  zAxis->SetTitleSize(AntarcticaBackgroundDefaults::zAxisTextSize);
147  zAxis->SetLabelSize(AntarcticaBackgroundDefaults::zAxisTextSize);
148  // std::cout << zAxis->GetTitleOffset() << std::endl;
149  // zAxis->SetTitleOffset(0.1);
150  gPad->Modified();
151  gPad->Update();
152  }
153 }
154 
155 
156 void TProfile2DAntarctica::ExecuteEvent(int event, int px, int py){
157  if(event==kButton1Double){
158  this->UnZoom();
159  }
160  AntarcticaBackground* b = getBackground();
161  if(b){
162  b->ExecuteEvent(event, px, py);
163  }
164  TProfile2D::ExecuteEvent(event, px, py);
165 }
166 
167 
168 TGraphAntarctica* TProfile2DAntarctica::findLocalMaxima() const {
169 
170  TGraphAntarctica* grLocalMaxima = new TGraphAntarctica();
171 
172  for(int by=2; by <= GetNbinsY()-1; by++){
173  for(int bx=2; bx <= GetNbinsX()-1; bx++){
174 
175  // check the 8 points around this point
176  double cVal = GetBinContent(bx, by);
177 
178  // Is this a local maximum?
179  if(cVal > GetBinContent(bx-1, by-1) &&
180  cVal > GetBinContent(bx-1, by ) &&
181  cVal > GetBinContent(bx-1, by+1) &&
182  cVal > GetBinContent(bx , by-1) &&
183  cVal > GetBinContent(bx , by+1) &&
184  cVal > GetBinContent(bx+1, by-1) &&
185  cVal > GetBinContent(bx+1, by ) &&
186  cVal > GetBinContent(bx+1, by+1)){
187 
188  Double_t easting = fXaxis.GetBinCenter(bx);
189  Double_t northing = fYaxis.GetBinCenter(by);
190  grLocalMaxima->TGraph::SetPoint(grLocalMaxima->GetN(), easting, northing);
191  }
192  }
193  }
194  return grLocalMaxima;
195 }
196 
197 TAxis* TProfile2DAntarctica::GetXaxis(){
198  AntarcticaBackground* b = getBackground();
199  if(b){
200  return b->GetXaxis();
201  }
202  else{
203  return &fXaxis;
204  }
205 }
206 
207 TAxis* TProfile2DAntarctica::GetYaxis(){
208  AntarcticaBackground* b = getBackground();
209  if(b){
210  return b->GetYaxis();
211  }
212  else{
213  return &fYaxis;
214  }
215 }
216 
217 
218 
219 
220 
221 
222 
223 
224 
225 
226 
227 
228 
229 
230 
231 
232 
233 
234 
235 
236 
237 TH2DAntarctica::TH2DAntarctica(Int_t nx, Int_t ny)
238  : TH2D(), fAntarcticaBackground(NULL){
239 
240  AntarcticaBackground* b = getBackground();
241 
242  nx = nx <= 0 ? b->GetNbinsX() : nx;
243  ny = ny <= 0 ? b->GetNbinsY() : ny;
244  SetBins(nx, b->GetXaxis()->GetBinLowEdge(1), b->GetXaxis()->GetBinUpEdge(b->GetNbinsX()),
245  ny, b->GetYaxis()->GetBinLowEdge(1), b->GetYaxis()->GetBinUpEdge(b->GetNbinsY()));
246 
247  accept_stereographic = false;
248 }
249 
250 
251 TH2DAntarctica::TH2DAntarctica(const char* name, const char* title, Int_t nx, Int_t ny)
252  : TH2D(), fAntarcticaBackground(NULL){
253 
254  SetNameTitle(name, title);
255  AntarcticaBackground* b = getBackground();
256 
257  nx = nx <= 0 ? b->GetNbinsX() : nx;
258  ny = ny <= 0 ? b->GetNbinsY() : ny;
259  SetBins(nx, b->GetXaxis()->GetBinLowEdge(1), b->GetXaxis()->GetBinUpEdge(b->GetNbinsX()),
260  ny, b->GetYaxis()->GetBinLowEdge(1), b->GetYaxis()->GetBinUpEdge(b->GetNbinsY()));
261  accept_stereographic = false;
262 }
263 
264 
265 
266 TH2DAntarctica::TH2DAntarctica(const char* name, const char* title, const std::vector<Double_t>& x, const std::vector<Double_t>& y)
267  : TH2D(), fAntarcticaBackground(NULL){
268 
269  SetNameTitle(name, title);
270  if(x.size() > 1 && y.size() > 1){
271  SetBins(x.size()-1, &x[0], y.size()-1, &y[0]);
272  }
273  else{
274  AntarcticaBackground* b = getBackground();
275 
276  Int_t nx = b->GetNbinsX();
277  Int_t ny = b->GetNbinsY();
278  SetBins(nx, b->GetXaxis()->GetBinLowEdge(1), b->GetXaxis()->GetBinUpEdge(b->GetNbinsX()),
279  ny, b->GetYaxis()->GetBinLowEdge(1), b->GetYaxis()->GetBinUpEdge(b->GetNbinsY()));
280 
281  }
282  accept_stereographic = false;
283 }
284 
285 
286 
287 
288 void TH2DAntarctica::Draw(Option_t* opt){
289 
290  TString opt2(opt);
291  if(!opt2.Contains("same")){
292  AntarcticaBackground* b = getBackground();
293  b->Draw();
294  b->SetBit(kCanDelete, false);
295  }
296  else{
297  TList* l = gPad->GetListOfPrimitives();
298  for(int i=0; i < l->GetEntries(); i++){
299  TString objName = l->At(i)->GetName();
300  if(objName.Contains("fAntarctica") && !objName.Contains("PalSetter")){
301  if(fAntarcticaBackground && fAntarcticaBackground != (AntarcticaBackground*) l->At(i)){
302  delete fAntarcticaBackground;
303  }
304  fAntarcticaBackground = (AntarcticaBackground*) l->At(i);
305  // break;
306  }
307  }
308  }
309  TString sameOpt = opt2.Contains("same") ? opt2 : TString::Format("%s same", opt);
310  TH2D::Draw(sameOpt);
311  ResetColorAxis();
312  RescaleBackground();
313  gPad->Modified();
314 }
315 
316 
317 Int_t TH2DAntarctica::Fill(Double_t lon, Double_t lat, Double_t val){
318 
319  Double_t easting, northing;
320 
321  if (accept_stereographic)
322  {
323  easting = lon;
324  northing = lat;
325  }
326  else
327  {
328  RampdemReader::LonLatToEastingNorthing(lon, lat, easting, northing);
329  }
330  return TH2D::Fill(easting, northing, val);
331 }
332 
333 void getAngularResolution(Double_t snr, Double_t& sigmaTheta, Double_t& sigmaPhi)
334 {
335  //this is just here for the error contour filling
336  sigmaPhi = (AnitaVersion::get() == 3) ? exp(-2.50414e-01 * snr + 3.02406e-01) + 2.43376e-01 : 5.09057/(pow(snr, 8.01369e-01) + 1.);
337  sigmaTheta = (AnitaVersion::get() == 3) ? exp(-3.83773e-01 * snr + -3.00964e-01) + 1.64537e-01 : 1.34307/(pow(snr, 7.09382e-01) + 1.);
338 }
339 
340 Int_t TH2DAntarctica::FillWithErrorContours(Double_t lon, Double_t lat, Double_t phi, Double_t theta,Double_t snr, Double_t ll_thresh, UsefulAdu5Pat upat, Double_t dist_thresh){
341  Double_t sigmaPhi, sigmaTheta;
342  getAngularResolution(snr, sigmaTheta, sigmaPhi);
343  return FillWithErrorContours(lon,lat,phi,theta,sigmaPhi,sigmaTheta,ll_thresh,upat,dist_thresh);
344 }
345 
346 
347 Int_t TH2DAntarctica::FillWithErrorContours(Double_t lon, Double_t lat, Double_t phi, Double_t theta,Double_t sigmaPhi, Double_t sigmaTheta, Double_t ll_thresh, UsefulAdu5Pat upat, Double_t dist_thresh){
348  //maxval is the value of the reconstructed point on the continent, other points are maxval - sqrt(ll)
349  Int_t maxval = ceil(sqrt(ll_thresh)) + 1;
350 
351  Double_t easting, northing;
352  RampdemReader::LonLatToEastingNorthing(lon, lat, easting, northing);
353 
354  Int_t binX = TH2D::GetXaxis()->FindBin(easting);
355  Int_t binY = TH2D::GetYaxis()->FindBin(northing);
356 
357  bool checking = true;
358  Int_t n_added = 1;
359 
360  //want to keep track of which ways it's still worth checking
361  int udlr[4];
362  bool updone = 0;
363  bool downdone = 0;
364  bool leftdone = 0;
365  bool rightdone = 0;
366  Int_t nrow;
367  Int_t ncol;
368  int level = 1;
369 
370  while(checking)
371  {
372  int miss_count = 0;
373  int curr = 0;
374  udlr[0] = 0;
375  udlr[1] = 0;
376  udlr[2] = 0;
377  udlr[3] = 0;
378  int dist_added = 0;
379  for(int i = 0; i < 8*level; i++)
380  {
381  Double_t thetaSource, phiSource;
382  Double_t sourceEasting, sourceNorthing;
383  Double_t sourceLon, sourceLat;
384  nrow = 3 + 2*(level-1);
385  ncol = nrow - 2;
386 
387  Int_t newBinX, newBinY;
388 
389  //search outwards from a point in squares
390  if(i < nrow){
391  if(updone){
392  miss_count++;
393  continue;
394  }
395  curr = 0;
396  newBinX = binX-level+i;
397  newBinY = binY+level;
398  }
399  else if(i < 2 * nrow){
400  if(downdone){
401  miss_count++;
402  continue;
403  }
404  curr = 1;
405  newBinX = binX-level+i-nrow;
406  newBinY = binY-level;
407  }
408  else{
409  Int_t alternator = i%2 ? 1 : -1;
410  if(alternator == 1 && rightdone){
411  miss_count++;
412  continue;
413  }
414  if(alternator == -1 && leftdone){
415  miss_count++;
416  continue;
417  }
418  curr = alternator == 1 ? 3 : 2;
419  Int_t ypos = level - 1 - (i - 2*nrow)/2;
420  newBinX = binX+(alternator*level);
421  newBinY = binY+ypos;
422  }
423  sourceEasting = TH2D::GetXaxis()->GetBinCenter(newBinX);
424  sourceNorthing = TH2D::GetYaxis()->GetBinCenter(newBinY);
425  RampdemReader::EastingNorthingToLonLat(sourceEasting, sourceNorthing, sourceLon, sourceLat);
426  Double_t sourceAlt = RampdemReader::SurfaceAboveGeoid(sourceLon, sourceLat);
427 
428  if(dist_thresh > 0)
429  {
430  Double_t surfaceDist = 1e-3*AntarcticCoord::surfaceDistance(sourceLat, lat, sourceLon, lon);
431  if(surfaceDist < dist_thresh)
432  {
433  Fill(sourceLon, sourceLat, 0.01);
434  dist_added = 1;
435  }
436  }
437 
438  upat.getThetaAndPhiWave2(sourceLon, sourceLat, sourceAlt, thetaSource, phiSource);
439 
440  thetaSource *= TMath::RadToDeg();
441  phiSource *= TMath::RadToDeg();
442 
443  Double_t dTheta = (thetaSource - theta)/sigmaTheta;
444  Double_t deltaPhi = phiSource - phi;
445  if(deltaPhi > 180) {
446  while(deltaPhi >= 180)
447  deltaPhi -= 360;
448  while(deltaPhi < -180)
449  deltaPhi += 360;
450  }
451  Double_t dPhi = deltaPhi/sigmaPhi;
452 
453  Double_t ll = dTheta * dTheta + dPhi * dPhi;
454  if(ll > ll_thresh)
455  {
456  if(!dist_added)
457  {
458  miss_count++;
459  udlr[curr]++;
460  }
461  else
462  {
463  n_added++;
464  }
465  continue;
466  }
467  Fill(sourceLon, sourceLat, maxval - ceil(sqrt(ll)));
468  n_added++;
469  }
470  if(miss_count == 8*level) checking = false;
471  if(udlr[0] == nrow) updone = true;
472  if(udlr[1] == nrow) downdone = true;
473  if(udlr[2] == ncol) leftdone = true;
474  if(udlr[3] == ncol) rightdone = true;
475  level++;
476  }
477  Fill(lon, lat, maxval);
478  return n_added;
479 }
480 
481 
487 void TH2DAntarctica::FillRandomly(Int_t nTimes){
488 
489  TRandom3 randy;
490 
491  for(int i=0; i < nTimes; i++){
492 
493  double lon = randy.Uniform(0, 360);
494  double lat = randy.Uniform(-90, -60);
495 
496  Fill(lon, lat);
497  }
498 }
499 
500 
504 void TH2DAntarctica::ResetColorAxis(){
505  if(gPad){
506  gPad->Modified();
507  gPad->Update();
508  TPaletteAxis *palette = (TPaletteAxis*) GetListOfFunctions()->FindObject("palette");
509  if(palette){
510  palette->SetX1NDC(palX1);
511  palette->SetX2NDC(palX2);
512  palette->SetY1NDC(palY1);
513  palette->SetY2NDC(palY2);
514  TAxis* zAxis = GetZaxis();
515  // zAxis->SetTitle();
516  zAxis->SetTitleSize(AntarcticaBackgroundDefaults::zAxisTextSize);
517  zAxis->SetLabelSize(AntarcticaBackgroundDefaults::zAxisTextSize);
518  // std::cout << zAxis->GetTitleOffset() << std::endl;
519  // zAxis->SetTitleOffset(0.1);
520  gPad->Modified();
521  gPad->Update();
522  }
523  }
524 }
525 
526 
527 void TH2DAntarctica::ExecuteEvent(int event, int px, int py){
528  if(event==kButton1Double){
529  this->UnZoom();
530  }
531  AntarcticaBackground* b = getBackground();
532  if(b){
533  b->ExecuteEvent(event, px, py);
534  }
535  TH2D::ExecuteEvent(event, px, py);
536 }
537 
538 
539 TGraphAntarctica* TH2DAntarctica::findLocalMaxima() const {
540 
541  TGraphAntarctica* grLocalMaxima = new TGraphAntarctica();
542 
543  for(int by=2; by <= GetNbinsY()-1; by++){
544  for(int bx=2; bx <= GetNbinsX()-1; bx++){
545 
546  // check the 8 points around this point
547  double cVal = GetBinContent(bx, by);
548 
549  // Is this a local maximum?
550  if(cVal > GetBinContent(bx-1, by-1) &&
551  cVal > GetBinContent(bx-1, by ) &&
552  cVal > GetBinContent(bx-1, by+1) &&
553  cVal > GetBinContent(bx , by-1) &&
554  cVal > GetBinContent(bx , by+1) &&
555  cVal > GetBinContent(bx+1, by-1) &&
556  cVal > GetBinContent(bx+1, by ) &&
557  cVal > GetBinContent(bx+1, by+1)){
558 
559  Double_t easting = fXaxis.GetBinCenter(bx);
560  Double_t northing = fYaxis.GetBinCenter(by);
561  grLocalMaxima->TGraph::SetPoint(grLocalMaxima->GetN(), easting, northing);
562  }
563  }
564  }
565  return grLocalMaxima;
566 }
567 
568 
569 TAxis* TH2DAntarctica::GetXaxis(){
570  AntarcticaBackground* b = getBackground();
571  if(b){
572  return b->GetXaxis();
573  }
574  else{
575  return &fXaxis;
576  }
577 }
578 
579 TAxis* TH2DAntarctica::GetYaxis(){
580  AntarcticaBackground* b = getBackground();
581  if(b){
582  return b->GetYaxis();
583  }
584  else{
585  return &fYaxis;
586  }
587 }
AntarcticaBackground::ExecuteEvent
void ExecuteEvent(Int_t event, Int_t x, Int_t y)
Definition: AntarcticaBackground.cxx:721
RampdemReader::SurfaceAboveGeoid
static Double_t SurfaceAboveGeoid(Double_t longitude, Double_t latitude, RampdemReader::dataSet=rampdem)
Definition: RampdemReader.cxx:249
AntarcticaBackground::Draw
virtual void Draw(Option_t *opt="")
Definition: AntarcticaBackground.cxx:494
TGraphAntarctica::getBackground
ClassDef(TGraphAntarctica, 1) private AntarcticaBackground * getBackground() const
Don&#39;t persist.
Definition: TGraphAntarctica.h:89
TH2DAntarctica::FillRandomly
void FillRandomly(Int_t nTimes=5000)
Definition: TH2DAntarctica.cxx:487
TProfile2DAntarctica::FillRandomly
void FillRandomly(Int_t nTimes=5000)
Definition: TH2DAntarctica.cxx:116
TH2DAntarctica::FillWithErrorContours
Int_t FillWithErrorContours(Double_t lon, Double_t lat, Double_t phi, Double_t theta, Double_t snr, Double_t ll_thresh, UsefulAdu5Pat upat, Double_t dist_thresh=0)
Definition: TH2DAntarctica.cxx:340
UsefulAdu5Pat
As UsefulAnitaEvent is to RawAnitaEvent, UsefulAdu5Pat is to Adu5Pat. Well not quite as useful but yo...
Definition: UsefulAdu5Pat.h:39
RampdemReader::LonLatToEastingNorthing
static void LonLatToEastingNorthing(Double_t lon, Double_t lat, Double_t &easting, Double_t &northing)
Definition: RampdemReader.cxx:536
AntarcticCoord::surfaceDistance
static double surfaceDistance(const AntarcticCoord &a, const AntarcticCoord &b)
Definition: AntarcticaGeometry.cxx:937
UsefulAdu5Pat::getThetaAndPhiWave2
void getThetaAndPhiWave2(Double_t sourceLon, Double_t sourceLat, Double_t sourceAlt, Double_t &thetaWave, Double_t &phiWave, TVector3 *sourcePos=NULL) const
Definition: UsefulAdu5Pat.cxx:758
RampdemReader::EastingNorthingToLonLat
static void EastingNorthingToLonLat(Double_t easting, Double_t northing, Double_t &lon, Double_t &lat)
Definition: RampdemReader.cxx:623
Generated on Tue Mar 23 2021 12:52:25 for ANITA Software by   doxygen 1.8.14