multiplicity.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 <TROOT.h>
10 #include <TStyle.h>
11 #include <TLegend.h>
12 #include <TColor.h>
13 #include <TTimeStamp.h>
14 #include <TEllipse.h>
15 #include <TVector3.h>
16 #include <TProfile.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 int main() {
32  TH2F* Multiplicity = new TH2F("Muliplicity", "Multiplicity", 15 , 60, 90, 75, 0, 75);
33  TProfile* Multiplicity_prof = new TProfile("Multiplicity_prof", "Multiplicity_prof", 15, 60, 90, 0, 75);
34 
35  TGraph* FurthestStation = new TGraph();
36  TGraph* FurthestStationWithCuts = new TGraph();
37  TProfile* FurthestStation_prof = new TProfile("FurthestStation_prof", "FurthestStation_prof", 15, 60, 90, 0, 3000,"S");
38 
39  TGraph* FurthestStationEnergy = new TGraph();
40  TGraph* FurthestStationEnergyWithCuts = new TGraph();
41  TProfile* FurthestStationEnergy_prof = new TProfile("FurthestStationEnergy_prof", "FurthestStationEnergy_prof", 10, 17.5, 19.5, 0, 2000,"S");
42 
43  TH2F* FurthestStationEnergyTheta = new TH2F("FurthestStationEnergyTheta", "FurthestStationEnergyTheta", 20, 10, 360, 2000, 0, 2000);
44  TProfile* FurthestStationEnergyTheta_prof = new TProfile("FurthestStationEnergyTheta_prof", "FurthestStationEnergyTheta_prof", 20, 10, 360, 0, 2000,"S");
45 
46  TH1F* NumberOfStations = new TH1F("NumberOfStations", "NumberOfStations", 75 , 0, 75);
47 
48  int dutch_stations[46]= {102, 104, 105, 108, 109, 112, 134, 135, 136, 143,
49  144, 149, 158, 159, 160, 169, 170, 183, 184, 185,
50  188, 192, 193, 194, 200, 201, 202, 203, 210, 211,
51  212, 213, 222, 223, 224, 225, 226, 235, 236, 237,
52  244, 245, 246, 251, 252, 257};
53 
54  vector<string> v_filename;
55  ifstream fileList ("./filelist.txt");
56  char filename[10000] ;
57  while (fileList >> filename) {
58  v_filename.push_back(filename);
59  }
60 
61  RecEventFile *file = new RecEventFile(v_filename);
62  RecEvent* theevent= new RecEvent();
63  int nevents=file->GetNEvents();
64  cout << nevents << endl;
65  file->SetBuffers(&theevent);
66  DetectorGeometry *geo =new DetectorGeometry();
67  file->ReadDetectorGeometry(*geo);
68 
69  //Loop over events
70  while(file->ReadNextEvent()==RecEventFile::eSuccess) {
71  const RdEvent& theradio= theevent->GetRdEvent();
72  const RdRecShower& rshow=theradio.GetRdRecShower();
73 
74  if (!theradio.HasRdStations()) continue ;
75  const std::vector<RdRecStation>& vecRstat=theradio.GetRdStationVector();
76  const SDEvent& sdev=theevent->GetSDEvent();
77  const SdRecShower& sshow= sdev.GetSdRecShower();
78 
79  const TVector3& score=sshow.GetCoreSiteCS();
80  const TVector3& rcore=rshow.GetCoreSiteCS();
81 
82  const TVector3& saxis=sshow.GetAxisSiteCS();
83 
84  double rtheta=TMath::RadToDeg()*rshow.GetZenith();
85  double rphi=TMath::RadToDeg()*rshow.GetAzimuth();
86  double stheta=TMath::RadToDeg()*sshow.GetZenith();
87  double sphi=TMath::RadToDeg()*sshow.GetAzimuth();
88 
89  if (stheta<60 || stheta>84)
90  continue;
91 
92  double RdZenith = rshow.GetZenith();
93  double RdAzimuth = rshow.GetAzimuth();
94  double SdZenith = sshow.GetZenith();
95  double SdAzimuth = sshow.GetAzimuth();
96  const double xr = sin(RdZenith)*cos(RdAzimuth);
97  const double yr = sin(RdZenith)*sin(RdAzimuth);
98  const double zr = cos(RdZenith);
99  const double xs = sin(SdZenith)*cos(SdAzimuth);
100  const double ys = sin(SdZenith)*sin(SdAzimuth);
101  const double zs = cos(SdZenith);
102  const double phi = 180./TMath::Pi()*acos((xr*xs+yr*ys+zr*zs)/(
103  sqrt(xr*xr+yr*yr+zr*zr)*sqrt(xs*xs+ys*ys+zs*zs)));
104  if (phi > 20){
105  cout << "opening angle SD RD greater 20 degree. Skipping event" << endl;
106  continue;
107  }
108 
109  //Multiplicity
110  if(rshow.HasParameter(revt::eNumberOfStationsWithPulseFound)) {
111  Multiplicity->Fill(stheta,rshow.GetParameter(
112  revt::eNumberOfStationsWithPulseFound));
113  Multiplicity_prof->Fill(stheta,rshow.GetParameter(
114  revt::eNumberOfStationsWithPulseFound), 1);
115  }
116 
117  //Furthest station with signal
118  double furthestDistance = 0.0;
119  for (vector<RdRecStation>::const_iterator stiter=vecRstat.begin();
120  stiter!=vecRstat.end(); ++stiter) {
121  /*
122  for(int i=0; i<46; i++)
123  if(((stiter->GetId())-100)==dutch_stations[i])
124  continue;
125  */
126  if (stiter->IsRejected())
127  continue; // Do not draw stations which are rejected
128 
129  double axis_distance_to_sd=geo->GetStationAxisDistance(
130  geo->GetRdStationPosition(stiter->GetId()), saxis , score);
131 
132  if(stiter->GetParameter(revt::eSignalToNoise)>10.0)
133  if (axis_distance_to_sd > furthestDistance)
134  furthestDistance = axis_distance_to_sd;
135  }
136 
137  //if (furthestDistance > 3000){
138  // cout << theradio.GetRdEventId() << " " << theradio.GetRdRunNumber()
139  // << " " << stheta << " " << furthestDistance << endl;
140  //}
141 
142  if(theradio.GetRdRunNumber() == 100809 && theradio.GetRdEventId() == 114253){
143  continue; //skip unrealistic event, prop. noise?
144  }
145 
146  FurthestStation_prof->Fill(stheta,furthestDistance,1);
147 
148  double energy_sd = sshow.GetEnergy();
149  FurthestStationEnergy_prof->Fill(log10(energy_sd),furthestDistance,1);
150 
151  FurthestStationEnergyTheta->Fill(log10(energy_sd)/cos(stheta*TMath::DegToRad()),
152  furthestDistance);
153  FurthestStationEnergyTheta_prof->Fill(log10(energy_sd)/cos(stheta*TMath::DegToRad()),
154  furthestDistance,1);
155 
156  NumberOfStations->Fill(rshow.GetParameter(revt::eNumberOfStationsWithPulseFound));
157 
158  //SD quality check
159  const int minRecLevel = 3; // has reconstructed energy
160  const int minCandidateStations = 4;
161  const int T4Trigger = 2; // has T4 and SD triggered independently of FD
162  const int T5Trigger = 3; //has T5 Posterior, which is the T5 HAS in this case
163 
164  const float minZenithSD = 60;
165  const float maxZenithSD = 84;
166 
167  if(sdev.GetRecLevel()<minRecLevel ||
168  sdev.GetNumberOfCandidates()<minCandidateStations ||
169  sdev.GetT4Trigger()<T4Trigger || sdev.GetT5Trigger()<T5Trigger ||
170  stheta>maxZenithSD || stheta<minZenithSD ||
171  log10(sshow.GetEnergy()) < 18.5){
172  FurthestStation->SetPoint(FurthestStation->GetN(),
173  stheta,furthestDistance);
174  FurthestStationEnergy->SetPoint(FurthestStation->GetN(),
175  log10(energy_sd),furthestDistance);
176  } else {
177  FurthestStationWithCuts->SetPoint(FurthestStationWithCuts->GetN(),
178  stheta,furthestDistance);
179  FurthestStationEnergyWithCuts->SetPoint(FurthestStation->GetN(),
180  log10(energy_sd),furthestDistance);
181  }
182  }
183 
184  cout << "number of entries, min, max" << endl;
185  cout << FurthestStation_prof->GetEntries() << endl;
186  cout << TMath::MinElement(FurthestStation->GetN(),FurthestStation->GetX()) << endl;
187  cout << TMath::MaxElement(FurthestStation->GetN(),FurthestStation->GetX()) << endl;
188 
189  //make plots
190  TCanvas Histos("Histos","Histos",1800,1000);
191  Histos.SetLeftMargin(0.14);
192  Histos.SetRightMargin(0.05);
193  Histos.SetTopMargin(0.06);
194  Histos.SetBottomMargin(0.14);
195 
196  FurthestStation->SetTitle(";#theta_{SD} / #circ;Axis distance of furthest station / m");
197  FurthestStation->SetMarkerStyle(33);
198  FurthestStation->SetMarkerSize(2.);
199  FurthestStation->SetMarkerColor(38);
200  FurthestStation->GetXaxis()->SetTitleSize(0.055);
201  FurthestStation->GetXaxis()->SetTitleOffset(1.2);
202  FurthestStation->GetXaxis()->SetRangeUser(59, 90);
203  FurthestStation->GetXaxis()->SetLabelFont(42);
204  FurthestStation->GetXaxis()->SetLabelSize(0.05);
205  FurthestStation->GetXaxis()->SetLabelOffset(0.01);
206  FurthestStation->GetYaxis()->SetTitleSize(0.055);
207  FurthestStation->GetYaxis()->SetTitleOffset(1.1);
208  FurthestStation->GetYaxis()->SetLabelFont(42);
209  FurthestStation->GetYaxis()->SetLabelSize(0.05);
210  FurthestStation->GetYaxis()->SetLabelOffset(0.01);
211  FurthestStation->Draw("AP");
212 
213  FurthestStationWithCuts->SetMarkerStyle(20);
214  FurthestStationWithCuts->SetMarkerSize(1.5);
215  FurthestStationWithCuts->SetMarkerColor(kBlack);
216  FurthestStationWithCuts->Draw("P SAMES");
217 
218  FurthestStation_prof->SetStats(kFALSE);
219  FurthestStation_prof->SetMarkerStyle(34);
220  FurthestStation_prof->SetMarkerSize(1.8);
221  FurthestStation_prof->SetLineWidth(2);
222  FurthestStation_prof->SetLineColor(kRed);
223  FurthestStation_prof->SetMarkerColor(kRed);
224  FurthestStation_prof->Draw("P SAMES");
225 
226  Histos.SaveAs("furtheststation_noCuts.pdf");
227 
228  TCanvas Histos2("Histos","Histos",1800,1000);
229  Histos2.SetLeftMargin(0.14);
230  Histos2.SetRightMargin(0.05);
231  Histos2.SetTopMargin(0.06);
232  Histos2.SetBottomMargin(0.14);
233 
234  FurthestStationEnergy->SetTitle(";log_{10}(E_{SD} / eV);Axis distance of furthest station / m");
235  FurthestStationEnergy->SetMarkerStyle(33);
236  FurthestStationEnergy->SetMarkerSize(2.);
237  FurthestStationEnergy->SetMarkerColor(38);
238  FurthestStationEnergy->GetXaxis()->SetTitleSize(0.055);
239  FurthestStationEnergy->GetXaxis()->SetTitleOffset(1.2);
240  FurthestStationEnergy->GetXaxis()->SetRangeUser(17.5, 19.5);
241  FurthestStationEnergy->GetXaxis()->SetLabelFont(42);
242  FurthestStationEnergy->GetXaxis()->SetLabelSize(0.05);
243  FurthestStationEnergy->GetXaxis()->SetLabelOffset(0.01);
244  FurthestStationEnergy->GetYaxis()->SetTitleSize(0.055);
245  FurthestStationEnergy->GetYaxis()->SetTitleOffset(1.1);
246  FurthestStationEnergy->GetYaxis()->SetLabelFont(42);
247  FurthestStationEnergy->GetYaxis()->SetLabelSize(0.05);
248  FurthestStationEnergy->GetYaxis()->SetLabelOffset(0.01);
249  FurthestStationEnergy->Draw("AP");
250 
251  FurthestStationEnergyWithCuts->SetMarkerStyle(20);
252  FurthestStationEnergyWithCuts->SetMarkerSize(1.5);
253  FurthestStationEnergyWithCuts->SetMarkerColor(kBlack);
254  FurthestStationEnergyWithCuts->Draw("P SAMES");
255 
256  FurthestStationEnergy_prof->SetStats(kFALSE);
257  FurthestStationEnergy_prof->SetMarkerStyle(34);
258  FurthestStationEnergy_prof->SetMarkerSize(1.8);
259  FurthestStationEnergy_prof->SetLineWidth(2);
260  FurthestStationEnergy_prof->SetLineColor(kRed);
261  FurthestStationEnergy_prof->SetMarkerColor(kRed);
262  FurthestStationEnergy_prof->Draw("P SAMES");
263 
264  Histos2.SaveAs("FurthestStationEnergy_noCuts.pdf");
265 
266  FurthestStationEnergyTheta->GetXaxis()->SetTitle("log10(Energy/eV)/cos(#theta)");
267  FurthestStationEnergyTheta->GetYaxis()->SetTitle("Axis distance of furthest station [m]");
268  FurthestStationEnergyTheta->SetMarkerStyle(20);
269  FurthestStationEnergyTheta->SetMarkerSize(0.4);
270  FurthestStationEnergyTheta->SetMarkerColor(kBlack);
271  FurthestStationEnergyTheta->GetXaxis()->SetTitleOffset(1.2);
272  FurthestStationEnergyTheta->GetYaxis()->SetTitleOffset(1.2);
273  FurthestStationEnergyTheta->Draw("P");
274 
275  FurthestStationEnergyTheta_prof->SetStats(kFALSE);
276  FurthestStationEnergyTheta_prof->SetMarkerStyle(34);
277  FurthestStationEnergyTheta_prof->SetMarkerSize(1.4);
278  FurthestStationEnergyTheta_prof->SetLineColor(kRed);
279  FurthestStationEnergyTheta_prof->SetMarkerColor(kBlue);
280  FurthestStationEnergyTheta_prof->Draw("P SAME");
281 
282  Histos.SaveAs("furtheststationenergytheta_noCuts.pdf");
283 }
main
int main(int argc, char *argv[])
Definition: DBSync.cc:58
file
const string file
Definition: testReadStream.cc:38
filename
char * filename
Definition: dump1090.h:266
, generated on Tue Sep 26 2023.