coredistribution.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <string>
4 #include <TCanvas.h>
5 #include <TH2F.h>
6 #include <TMath.h>
7 #include <TFile.h>
8 #include <TGraph.h>
9 #include <TGraph2D.h>
10 #include <TROOT.h>
11 #include <TStyle.h>
12 #include <TColor.h>
13 #include <TTimeStamp.h>
14 #include <TVector3.h>
15 #include <TProfile.h>
16 #include <TLine.h>
17 
18 #include "RecEventFile.h"
19 #include "RecEvent.h"
20 #include "RdEvent.h"
21 #include "DetectorGeometry.h"
22 #include "SDEvent.h"
23 #include "DetectorGeometry.h"
24 #include "SdRecStation.h"
25 #include "ShowerRRecDataQuantities.h"
26 #include "StationRRecDataQuantities.h"
27 #include "RdTrace.h"
28 
29 using namespace std;
30 
31 //from https://root-forum.cern.ch/t/set-real-aspect-ratio-of-th2-or-tgraph/20734
45 bool SetRealAspectRatio(TCanvas* const c, const Int_t axis = 1)
46 {
47  if(!c)
48  {
49  cout << "Error in SetRealAspectRatio: canvas is NULL";
50  return false;
51  }
52 
53  {
54  //Get the current min-max values if SetUserRange has been called
55  c->Update();
56  const Double_t xmin = c->GetUxmin();
57  const Double_t xmax = c->GetUxmax();
58  const Double_t ymin = c->GetUymin();
59  const Double_t ymax = c->GetUymax();
60 
61  //Get the length of zoomed x and y axes
62  const Double_t xlength = xmax - xmin;
63  const Double_t ylength = ymax - ymin;
64  const Double_t ratio = xlength/ylength;
65 
66  //Get how many pixels are occupied by the canvas
67  const Int_t npx = c->GetWw();
68  const Int_t npy = c->GetWh();
69 
70  //Get x-y coordinates at the edges of the canvas (extrapolating outside the axes, NOT at the edges of the histogram)
71  const Double_t x1 = c->GetX1();
72  const Double_t y1 = c->GetY1();
73  const Double_t x2 = c->GetX2();
74  const Double_t y2 = c->GetY2();
75 
76  //Get the length of extrapolated x and y axes
77  const Double_t xlength2 = x2 - x1;
78  const Double_t ylength2 = y2 - y1;
79  const Double_t ratio2 = xlength2/ylength2;
80 
81  //Get now number of pixels including window borders
82  const Int_t bnpx = c->GetWindowWidth();
83  const Int_t bnpy = c->GetWindowHeight();
84 
85  if(axis==1)
86  {
87  c->SetCanvasSize(TMath::Nint(npy*ratio2), npy);
88  c->SetWindowSize((bnpx-npx)+TMath::Nint(npy*ratio2), bnpy);
89  }
90  else if(axis==2)
91  {
92  c->SetCanvasSize(npx, TMath::Nint(npx/ratio2));
93  c->SetWindowSize(bnpx, (bnpy-npy)+TMath::Nint(npx/ratio2));
94  }
95  else
96  {
97  cout << "Error in SetRealAspectRatio: axis value " << axis << " is neither 1 (resize along x-axis) nor 2 (resize along y-axis).";
98  return false;
99  }
100  }
101 
102  //Check now that resizing has worked
103  {
104  //Get the current min-max values if SetUserRange has been called
105  c->Update();
106  const Double_t xmin = c->GetUxmin();
107  const Double_t xmax = c->GetUxmax();
108  const Double_t ymin = c->GetUymin();
109  const Double_t ymax = c->GetUymax();
110 
111  //Get the length of zoomed x and y axes
112  const Double_t xlength = xmax - xmin;
113  const Double_t ylength = ymax - ymin;
114  const Double_t ratio = xlength/ylength;
115 
116  //Get how many pixels are occupied by the canvas
117  const Int_t npx = c->GetWw();
118  const Int_t npy = c->GetWh();
119 
120  //Get x-y coordinates at the edges of the canvas (extrapolating outside the axes, NOT at the edges of the histogram)
121  const Double_t x1 = c->GetX1();
122  const Double_t y1 = c->GetY1();
123  const Double_t x2 = c->GetX2();
124  const Double_t y2 = c->GetY2();
125 
126  //Get the length of extrapolated x and y axes
127  const Double_t xlength2 = x2 - x1;
128  const Double_t ylength2 = y2 - y1;
129  const Double_t ratio2 = xlength2/ylength2;
130 
131  //Get now number of pixels including window borders
132  const Int_t bnpx = c->GetWindowWidth();
133  const Int_t bnpy = c->GetWindowHeight();
134 
135  //Check accuracy +/-1 pixel due to rounding
136  if(abs(TMath::Nint(npy*ratio2) - npx)<2)
137  {
138  cout << "Resizing finished successfully." << endl;
139  return true;
140  }
141  else
142  {
143  cout << "Resizing failed." << endl;
144  return false;
145  }
146  }
147  // References:
148  // https://root.cern.ch/root/roottalk/roottalk01/3676.html
149  // https://root-forum.cern.ch/t/making-the-both-axes-square-on-the-pad/4325/1
150 }
151 
152 int main() {
153  vector<string> v_filename;
154  ifstream fileList ("./filelist.txt");
155  char filename[10000] ;
156  while (fileList >> filename) {
157  v_filename.push_back(filename);
158  }
159 
160  RecEventFile *file = new RecEventFile(v_filename);
161  RecEvent* theevent= new RecEvent();
162  int nevents=file->GetNEvents();
163  cout << nevents << endl;
164  file->SetBuffers(&theevent);
165  DetectorGeometry *geo =new DetectorGeometry();
166  file->ReadDetectorGeometry(*geo);
167 
168  float distance=0;
169  float distance_before=0;
170  int distance_stadtion_id=0;
171 
172  double station_x[] = {-26348.3,-26493.3,-26607,-26347.7,-26478.8,-26225.9,-26361.8,-26109.3,-26485.5,-26229.1,
173  -25984.7,-26365.4,-26104.7,-26236.9,-26356.2,-26109.9,-26094.2,-25848.3,-25983.3,-25856.8,
174  -26233.4,-25989.5,-27108,-26857.8,-26605.2,-26357.5,-27491,-27223.5,-26981.8,-26737.3,-26484.5,
175  -26231.9,-27608.2,-27358.9,-27108.8,-26859.7,-28479.5,-28231.6,-27980.1,-27728.9,-27488.2,-27229.2,
176  -26980.3,-26733.5,-28357.1,-28113.1,-27850.6,-27609.5,-27356.8,-27113.6,-26853.6,-26609.5,-28464.4,
177  -28238.8,-27983.1,-27739.2,-27481.5,-27232.9,-26985.7,-26735.6,-26478.5,-26233.1,-25986.4,-25733.6,
178  -26044.2,-26603.2,-26241.8,-27548.6,-27168.3,-26797.2,-26416.9,-28481.7,-28104.3,-27737.5,-27368.3,
179  -28296.5,-27929};
180 
181  double station_y[] = {15594.1,15519.2,15456.6,15450,15381.8,15378.2,15309.3,15314.2,15230.4,15236.5,15236.1,
182  15167.3,15135.9,15093.5,15020.5,15017.8,15427.1,15157.7,15088.1,15016.8,14946.4,14945.7,
183  15886.9,15898.5,15887.8,15887.6,15669.3,15673.4,15669.7,15659.2,15677.4,15669.4,15444.9,
184  15449.4,15457,15450.6,15244.7,15240.4,15234.7,15242.5,15234.9,15241.5,15229.4,15236,15016.9,
185  15005.2,15018.3,15017.4,15018.5,15021.3,15013.7,15013,14876.1,14856,14833.2,14809.8,14806.9,
186  14799.1,14799.4,14793.5,14798.5,14798.3,14798.2,14796.2,16750.8,16424.2,16432.8,16097,16108,
187  16101.6,16100.4,15781.6,15778,15780.8,15771.9,15451.2,15452.9};
188 
189  for(int i=0; i<77; i++){
190  station_x[i]/=1000;
191  station_y[i]/=1000;
192  }
193 
194  TGraph* aera2stations = new TGraph(77, station_x, station_y);
195  TGraph* corepos_withCuts = new TGraph();
196  TGraph* corepos_withoutCuts = new TGraph();
197  TGraph* example_event = new TGraph();
198 
199  //Loop over events
200  while(file->ReadNextEvent()==RecEventFile::eSuccess) {
201  const RdEvent& theradio= theevent->GetRdEvent();
202  const RdRecShower& rshow=theradio.GetRdRecShower();
203 
204  if (!theradio.HasRdStations()) continue ;
205  const std::vector<RdRecStation>& vecRstat=theradio.GetRdStationVector();
206  const SDEvent& sdev=theevent->GetSDEvent();
207 
208  const SdRecShower& sshow= sdev.GetSdRecShower();
209  const TVector3& score=sshow.GetCoreSiteCS();
210  const TVector3& rcore=rshow.GetCoreSiteCS();
211  const TVector3& saxis=sshow.GetAxisSiteCS();
212 
213  double rtheta=TMath::RadToDeg()*rshow.GetZenith();
214  double rphi=TMath::RadToDeg()*rshow.GetAzimuth();
215  double stheta=TMath::RadToDeg()*sshow.GetZenith();
216  double sphi=TMath::RadToDeg()*sshow.GetAzimuth();
217 
218  if (stheta<60 || stheta>84)
219  continue;
220 
221  if(theradio.GetRdRunNumber() == 100785 && theradio.GetRdEventId() == 137009){
222  example_event->SetPoint(1, score.X()/1000, score.Y()/1000);
223  continue;
224  }
225 
226  //SD quality check
227  const int minRecLevel=3; // has reconstructed energy
228  const int minCandidateStations=4;
229  const int T4Trigger=2; // has T4 and SD triggered independently of FD
230  const int T5Trigger=3; //has T5 Posterior, which is the T5 HAS in this case
231 
232  const float minZenithSD=60;
233  const float maxZenithSD=84;
234 
235  if(sdev.GetRecLevel()<minRecLevel ||
236  sdev.GetNumberOfCandidates()<minCandidateStations ||
237  sdev.GetT4Trigger()<T4Trigger || sdev.GetT5Trigger()<T5Trigger ||
238  stheta>maxZenithSD || stheta<minZenithSD ||
239  log10(sshow.GetEnergy()) < 18.5){
240  corepos_withoutCuts->SetPoint(corepos_withoutCuts->GetN(),
241  score.X()/1000, score.Y()/1000);
242  } else {
243  corepos_withCuts->SetPoint(corepos_withCuts->GetN(),
244  score.X()/1000, score.Y()/1000);
245  }
246 
247  if (score.Y()/1000 < 1){
248  cout << "far event id: Run number, event id" << endl;
249  cout << theradio.GetRdRunNumber() << " " << theradio.GetRdEventId() << endl;
250  }
251 
252  /*if (sdev.GetRecLevel()<minRecLevel) continue;
253  if (sdev.GetNumberOfCandidates()<minCandidateStations) continue;
254  if (sdev.GetT4Trigger()<T4Trigger) continue;
255  if (sdev.GetT5Trigger()<T5Trigger) continue;
256  if (stheta>maxZenithSD) continue;
257  if (stheta<minZenithSD) continue;*/
258 
259  }
260 
261  TFile *f = new TFile("corepositions.root","RECREATE");
262  corepos_withoutCuts->Write();
263  corepos_withCuts->Write();
264  aera2stations->Write();
265  example_event->Write();
266  f->Close();
267 
268  TCanvas canvas("c","c",-1800,1800);
269  canvas.SetFixedAspectRatio();
270  canvas.SetLeftMargin(0.14);
271  //canvas.SetRightMargin(0.02);
272  canvas.SetRightMargin(0.14);
273  canvas.SetTopMargin(0.03);
274  canvas.SetBottomMargin(0.1);
275 
276  /*Double_t stops[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000,
277  0.6250, 0.7500, 0.8750, 1.0000};
278  double alpha = 1.;
279  Double_t red[9] = { 0./255., 48./255., 119./255., 173./255., 212./255.,
280  224./255., 228./255., 228./255., 245./255.};
281  Double_t green[9] = { 0./255., 13./255., 30./255., 47./255., 79./255.,
282  127./255., 167./255., 205./255., 245./255.};
283  Double_t blue[9] = { 0./255., 68./255., 75./255., 43./255., 16./255.,
284  22./255., 55./255., 128./255., 245./255.};
285 
286  //reverse(red, 9);
287  //reverse(green, 9);
288  //reverse(blue, 9);
289 
290  TColor::CreateGradientColorTable(9, stops, red, green, blue, 255, alpha);*/
291 
292  corepos_withoutCuts->SetMarkerStyle(33);
293  corepos_withoutCuts->SetMarkerColor(38);
294  corepos_withoutCuts->SetMarkerSize(2.);
295  corepos_withoutCuts->SetTitle(";x / km;y / km");
296  corepos_withoutCuts->GetXaxis()->SetLabelFont(42);
297  corepos_withoutCuts->GetXaxis()->SetLabelOffset(0);
298  corepos_withoutCuts->GetXaxis()->SetLabelSize(0.05);
299  corepos_withoutCuts->GetXaxis()->SetTitleSize(0.055);
300  corepos_withoutCuts->GetXaxis()->SetTitleOffset(0.75);
301  corepos_withoutCuts->GetXaxis()->SetNdivisions(505);
302  //corepos_withoutCuts->GetXaxis()->SetRangeUser(-32, -17);
303  //corepos_withoutCuts->GetYaxis()->SetRangeUser(-1, 25);
304  corepos_withoutCuts->GetYaxis()->SetLabelFont(42);
305  corepos_withoutCuts->GetYaxis()->SetLabelOffset(0.015);
306  corepos_withoutCuts->GetYaxis()->SetLabelSize(0.05);
307  corepos_withoutCuts->GetYaxis()->SetTitleSize(0.055);
308  corepos_withoutCuts->GetYaxis()->SetTitleOffset(1.15);
309  corepos_withoutCuts->Draw("AP");
310 
311  gPad->SetTheta(90.0);
312  gPad->SetPhi(0.0);
313 
314  corepos_withCuts->SetMarkerStyle(20);
315  corepos_withCuts->SetMarkerColor(kBlack);
316  corepos_withCuts->SetMarkerSize(1.5);
317  corepos_withCuts->Draw("p SAMES");
318 
319  aera2stations->SetMarkerStyle(22);
320  aera2stations->SetMarkerColor(TColor::GetColor("#d62728"));
321  aera2stations->SetMarkerSize(1.6);
322  aera2stations->Draw("p SAMES");
323 
324  example_event->SetMarkerStyle(kFullCross);
325  example_event->SetMarkerColor(TColor::GetColor("#ff7f0e"));
326  example_event->SetMarkerSize(3.);
327  example_event->Draw("p SAMES");
328 
329  gPad->Update();
330 
331  cout << "min x and max y for calc. of end points SD border" << endl;
332  cout << gPad->GetUymax() << " " << gPad->GetUxmin() << endl;
333 
334  //SD border line
335  //TLine* line = new TLine(-37.13878153, 3.116284307, -33.40331677, -3.366836367);
336  //line->Draw();
337  //line = new TLine(-37.13878153, 3.116284307, -29.62013503, 16.10875224);
338  //TLine* line = new TLine(-33.37943414, 9.627520781, -29.62013503, 16.10875224);
339  TLine* line = new TLine(-32.801, 10.624771980074598, -29.62013503, 16.10875224);
340  line->SetLineStyle(9);
341  line->SetLineColor(16);
342  line->SetLineWidth(2);
343  line->Draw();
344  line = new TLine(-29.62013503, 16.10875224, -30.37781708, 17.31638591);
345  line->SetLineStyle(9);
346  line->SetLineColor(16);
347  line->SetLineWidth(2);
348  line->Draw();
349  //line = new TLine(-30.37781708, 17.31638591, -25.06886873, 26.46995979);
350  line = new TLine(-30.37781708, 17.31638591, -25.145403613376814, 26.338);
351  line->SetLineStyle(9);
352  line->SetLineColor(16);
353  line->SetLineWidth(2);
354  line->Draw();
355 
356  SetRealAspectRatio(&canvas, 1);
357 
358  canvas.SaveAs("coredistribution.C");
359  canvas.SaveAs("coredistribution.pdf");
360 }
SetRealAspectRatio
bool SetRealAspectRatio(TCanvas *const c, const Int_t axis=1)
Function to resize a canvas so that a 2D histogram or TGraph / TGraph2D are shown in real aspect rati...
Definition: coredistribution.cc:45
main
int main(int argc, char *argv[])
Definition: DBSync.cc:58
utl::abs
double abs(const SVector< n, T > &v)
Definition: ValidationTests/SdCalibration/Validatrix.cc:32
file
const string file
Definition: testReadStream.cc:38
filename
char * filename
Definition: dump1090.h:266
, generated on Tue Sep 26 2023.