preselect.cc
Go to the documentation of this file.
1 #include <TROOT.h>
2 #include <TFile.h>
3 #include <TMath.h>
4 #include <TVector3.h>
5 #include <TSystem.h>
6 #include <TH1D.h>
7 #include <TH2D.h>
8 #include <TProfile.h>
9 
10 #include <iomanip>
11 #include <iostream>
12 #include <string>
13 #include <cmath>
14 #include <sstream>
15 #include <vector>
16 #include <algorithm>
17 #include <cmath>
18 #include <cstdio>
19 #include <iterator>
20 #include <fstream>
21 #include <algorithm>
22 #include <map>
23 
24 #include <RecEvent.h>
25 #include <RecEventFile.h>
26 #include <FileInfo.h>
27 #include <DetectorGeometry.h>
28 #include <EyeGeometry.h>
29 #include <FDEvent.h>
30 #include <FdRecLevel.h>
31 #include <EventInfo.h>
32 #include <Shower.h>
33 #include <GenShower.h>
34 #include <FdGenShower.h>
35 #include <TelescopeGeometry.h>
36 #include <FdTelescopeData.h>
37 #include <FdApertureLight.h>
38 #include <FdGeometry.h>
39 #include <FdApertureLight.h>
40 
41 #include <boost/algorithm/string.hpp>
42 #include <boost/lexical_cast.hpp>
43 #include <boost/filesystem.hpp>
44 
45 #include <chrono>
46 #include <thread>
47 
48 using namespace std;
49 
50 bool doCompress = false;
51 
52 void
53 preselect(boost::filesystem::path steerdir,
54  const vector<string>& inputs)
55 {
56  const double degree = 180 / M_PI;
57 
58 
59  // output file for ROOT
60  TFile* outPlots = TFile::Open((steerdir/"preselect_plots.root").string().c_str(), "update");
61  outPlots->cd();
62 
63 
64  // first try to read histograms from ROOT
65  TH1::SetDefaultSumw2();
66 
67  TH1D* hRecLevel = (TH1D*) outPlots->Get("hRecLevel");
68 
69  TH1D* hEnergyGen = (TH1D*) outPlots->Get("hEnergyGen");
70  TH1D* hEnergyTrg = (TH1D*) outPlots->Get("hEnergyTrg");
71  TH1D* hEnergyRec = (TH1D*) outPlots->Get("hEnergyRec");
72  TH1D* hEnergySel = (TH1D*) outPlots->Get("hEnergySel");
73 
74  TH1D* hXmaxGen = (TH1D*) outPlots->Get("hXmaxGen");
75  TH1D* hXmaxTrg = (TH1D*) outPlots->Get("hXmaxTrg");
76  TH1D* hXmaxRec = (TH1D*) outPlots->Get("hXmaxRec");
77  TH1D* hXmaxSel = (TH1D*) outPlots->Get("hXmaxSel");
78 
79  TH1D* hRpGen = (TH1D*) outPlots->Get("hRpGen");
80  TH1D* hRpTrg = (TH1D*) outPlots->Get("hRpTrg");
81  TH1D* hRpRec = (TH1D*) outPlots->Get("hRpRec");
82  TH1D* hRpSel = (TH1D*) outPlots->Get("hRpSel");
83 
84  TH1D* hZenithGen = (TH1D*) outPlots->Get("hZenithGen");
85  TH1D* hZenithTrg = (TH1D*) outPlots->Get("hZenithTrg");
86  TH1D* hZenithRec = (TH1D*) outPlots->Get("hZenithRec");
87  TH1D* hZenithSel = (TH1D*) outPlots->Get("hZenithSel");
88 
89  TH1D* hAngleGen = (TH1D*) outPlots->Get("hAngleGen");
90  TH1D* hAngleTrg = (TH1D*) outPlots->Get("hAngleTrg");
91  TH1D* hAngleRec = (TH1D*) outPlots->Get("hAngleRec");
92  TH1D* hAngleSel = (TH1D*) outPlots->Get("hAngleSel");
93 
94  TH1D* hLightTotalGen = (TH1D*) outPlots->Get("hLightTotalGen");
95  TH1D* hLightTotalTrg = (TH1D*) outPlots->Get("hLightTotalTrg");
96  TH1D* hLightTotalRec = (TH1D*) outPlots->Get("hLightTotalRec");
97  TH1D* hLightTotalSel = (TH1D*) outPlots->Get("hLightTotalSel");
98 
99  TH1D* hLightCherGen = (TH1D*) outPlots->Get("hLightCherGen");
100  TH1D* hLightCherTrg = (TH1D*) outPlots->Get("hLightCherTrg");
101  TH1D* hLightCherRec = (TH1D*) outPlots->Get("hLightCherRec");
102  TH1D* hLightCherSel = (TH1D*) outPlots->Get("hLightCherSel");
103 
104  TH1D* hCFracGen = (TH1D*) outPlots->Get("hCFracGen");
105  TH1D* hCFracTrg = (TH1D*) outPlots->Get("hCFracTrg");
106  TH1D* hCFracRec = (TH1D*) outPlots->Get("hCFracRec");
107  TH1D* hCFracSel = (TH1D*) outPlots->Get("hCFracSel");
108 
109  TH2D* hPosGen = (TH2D*) outPlots->Get("hPosGen");
110  TH2D* hPosTrg = (TH2D*) outPlots->Get("hPosTrg");
111  TH2D* hPosRec = (TH2D*) outPlots->Get("hPosRec");
112  TH2D* hPosSel = (TH2D*) outPlots->Get("hPosSel");
113 
114  TH2D* hDirGen = (TH2D*) outPlots->Get("hDirGen");
115  TH2D* hDirTrg = (TH2D*) outPlots->Get("hDirTrg");
116  TH2D* hDirRec = (TH2D*) outPlots->Get("hDirRec");
117  TH2D* hDirSel = (TH2D*) outPlots->Get("hDirSel");
118 
119  TH2D* hChFractXmax = (TH2D*) outPlots->Get("hChFractXmax");
120  TH2D* hChFractMean = (TH2D*) outPlots->Get("hChFractMean");
121  TH2D* hChFractMin = (TH2D*) outPlots->Get("hChFractMin");
122 
123  TProfile* proCherXmax = (TProfile*) outPlots->Get("proCherXmax");
124  TProfile* proCherMean = (TProfile*) outPlots->Get("proCherMean");
125 
126  TProfile* proXmax = (TProfile*) outPlots->Get("proXmax");
127  TProfile* proEnergy = (TProfile*) outPlots->Get("proEnergy");
128  TProfile* proChi0 = (TProfile*) outPlots->Get("proChi0");
129  TProfile* proAngle = (TProfile*) outPlots->Get("proAngle");
130 
131 
132  // if histogram reading is not successful, create new ones:
133  TH1::SetDefaultSumw2();
134 
135  if (!hRecLevel) hRecLevel = new TH1D("hRecLevel", "hRecLevel", 16, -1, 15);
136 
137  if (!hEnergyGen) hEnergyGen = new TH1D("hEnergyGen", "hEnergyGen", 100, 15, 18);
138  if (!hEnergyTrg) hEnergyTrg = new TH1D("hEnergyTrg", "hEnergyTrg", 100, 15, 18);
139  if (!hEnergyRec) hEnergyRec = new TH1D("hEnergyRec", "hEnergyRec", 100, 15, 18);
140  if (!hEnergySel) hEnergySel = new TH1D("hEnergySel", "hEnergySel", 100, 15, 18);
141 
142  if (!hXmaxGen) hXmaxGen = new TH1D("hXmaxGen", "hXmaxGen", 100, 300, 1000);
143  if (!hXmaxTrg) hXmaxTrg = new TH1D("hXmaxTrg", "hXmaxTrg", 100, 300, 1000);
144  if (!hXmaxRec) hXmaxRec = new TH1D("hXmaxRec", "hXmaxRec", 100, 300, 1000);
145  if (!hXmaxSel) hXmaxSel = new TH1D("hXmaxSel", "hXmaxSel", 100, 300, 1000);
146 
147  if (!hRpGen) hRpGen = new TH1D("hRpGen", "hRpGen", 100, 0, 5000);
148  if (!hRpTrg) hRpTrg = new TH1D("hRpTrg", "hRpTrg", 100, 0, 5000);
149  if (!hRpRec) hRpRec = new TH1D("hRpRec", "hRpRec", 100, 0, 5000);
150  if (!hRpSel) hRpSel = new TH1D("hRpSel", "hRpSel", 100, 0, 5000);
151 
152  if (!hZenithGen) hZenithGen = new TH1D("hZenithGen", "hZenithGen", 100, 0, 90);
153  if (!hZenithTrg) hZenithTrg = new TH1D("hZenithTrg", "hZenithTrg", 100, 0, 90);
154  if (!hZenithRec) hZenithRec = new TH1D("hZenithRec", "hZenithRec", 100, 0, 90);
155  if (!hZenithSel) hZenithSel = new TH1D("hZenithSel", "hZenithSel", 100, 0, 90);
156 
157  if (!hAngleGen) hAngleGen = new TH1D("hAngleGen", "hAngleGen", 100, -10, 120);
158  if (!hAngleTrg) hAngleTrg = new TH1D("hAngleTrg", "hAngleTrg", 100, -10, 120);
159  if (!hAngleRec) hAngleRec = new TH1D("hAngleRec", "hAngleRec", 100, -10, 120);
160  if (!hAngleSel) hAngleSel = new TH1D("hAngleSel", "hAngleSel", 100, -10, 120);
161 
162  if (!hLightTotalGen) hLightTotalGen = new TH1D("hLightTotalGen", "hLightTotalGen", 100, -2, 6.5);
163  if (!hLightTotalTrg) hLightTotalTrg = new TH1D("hLightTotalTrg", "hLightTotalTrg", 100, -2, 6.5);
164  if (!hLightTotalRec) hLightTotalRec = new TH1D("hLightTotalRec", "hLightTotalRec", 100, -2, 6.5);
165  if (!hLightTotalSel) hLightTotalSel = new TH1D("hLightTotalSel", "hLightTotalSel", 100, -2, 6.5);
166 
167  if (!hLightCherGen) hLightCherGen = new TH1D("hLightCherGen", "hLightCherGen", 100, -3, 6.5);
168  if (!hLightCherTrg) hLightCherTrg = new TH1D("hLightCherTrg", "hLightCherTrg", 100, -3, 6.5);
169  if (!hLightCherRec) hLightCherRec = new TH1D("hLightCherRec", "hLightCherRec", 100, -3, 6.5);
170  if (!hLightCherSel) hLightCherSel = new TH1D("hLightCherSel", "hLightCherSel", 100, -3, 6.5);
171 
172  if (!hCFracGen) hCFracGen = new TH1D("hCFracGen", "hCFracGen", 100, 0, 1);
173  if (!hCFracTrg) hCFracTrg = new TH1D("hCFracTrg", "hCFracTrg", 100, 0, 1);
174  if (!hCFracRec) hCFracRec = new TH1D("hCFracRec", "hCFracRec", 100, 0, 1);
175  if (!hCFracSel) hCFracSel = new TH1D("hCFracSel", "hCFracSel", 100, 0, 1);
176 
177  if (!hPosGen) hPosGen = new TH2D("hPosGen", "hPosGen", 100, -45000, -20000, 100, 3000, 28000);
178  if (!hPosTrg) hPosTrg = new TH2D("hPosTrg", "hPosTrg", 100, -45000, -20000, 100, 3000, 28000);
179  if (!hPosRec) hPosRec = new TH2D("hPosRec", "hPosRec", 100, -45000, -20000, 100, 3000, 28000);
180  if (!hPosSel) hPosSel = new TH2D("hPosSel", "hPosSel", 100, -45000, -20000, 100, 3000, 28000);
181 
182  if (!hDirGen) hDirGen = new TH2D("hDirGen", "hDirGen", 100, 0, 360, 100, 0, 90);
183  if (!hDirTrg) hDirTrg = new TH2D("hDirTrg", "hDirTrg", 100, 0, 360, 100, 0, 90);
184  if (!hDirRec) hDirRec = new TH2D("hDirRec", "hDirRec", 100, 0, 360, 100, 0, 90);
185  if (!hDirSel) hDirSel = new TH2D("hDirSel", "hDirSel", 100, 0, 360, 100, 0, 90);
186 
187  if (!hChFractXmax) hChFractXmax = new TH2D("hChFractXmax", "hChFractXmax", 100, -10, 120, 100, 0, 1);
188  if (!hChFractMean) hChFractMean = new TH2D("hChFractMean", "hChFractMean", 100, -10, 120, 100, 0, 1);
189  if (!hChFractMin) hChFractMin = new TH2D("hChFractMin", "hChFractMin", 100, -10, 120, 100, 0, 1);
190 
191  if (!proCherXmax) proCherXmax = new TProfile("proCherXmax", "proCherXmax", 100, -10, 120);
192  if (!proCherMean) proCherMean = new TProfile("proCherMean", "proCherMean", 100, -10, 120);
193 
194  if (!proEnergy) proEnergy = new TProfile("proEnergy", "proEnergy", 10, 0, 1);
195  if (!proXmax) proXmax = new TProfile("proXmax", "proXmax", 10, 0, 1);
196  if (!proChi0) proChi0 = new TProfile("proChi0", "proChi0", 10, 0, 1);
197  if (!proAngle) proAngle = new TProfile("proAngle", "proAngle", 10, 0, 1);
198 
199 
200  // open input ADST files
201  RecEvent* theRecEvent = new RecEvent();
202  RecEvent* theRecEventOut = new RecEvent();
203  DetectorGeometry* theGeometry = 0;
204  theGeometry = new DetectorGeometry();
205  FileInfo theInfo;
206 
207 
208  // list of counter variables
209  unsigned int countFile = 0;
210  unsigned int countEvent = 0;
211  unsigned int countSelected = 0;
212  unsigned int countSelectedAndReco = 0;
213  unsigned int countHeatFOV = 0;
214  unsigned int countMiss = 0; // count missing Cherenkov events
215 
216 
217  // loop over input files
218  for (const auto& input : inputs) {
219 
220  ++countFile;
221 
222  RecEventFile dataFile(input.c_str());
223  dataFile.SetBuffers(&theRecEvent);
224  dataFile.ReadDetectorGeometry(*theGeometry);
225  dataFile.ReadFileInfo(theInfo);
226 
227  ostringstream outputName;
228  outputName << input.substr(0, input.rfind(".root")) << "_presel.root";
229  if (boost::filesystem::exists(outputName.str())) {
230  cerr << "Output ADST (presel) " << outputName.str() << " does already exist. Do not overwrite" << endl;
231  exit(1);
232  }
233  RecEventFile preselADST(outputName.str().c_str(), RecEventFile::eWrite);
234  preselADST.SetBuffers(&theRecEventOut);
235  preselADST.WriteDetectorGeometry(*theGeometry);
236  preselADST.WriteFileInfo(theInfo);
237 
238 
239  const double HeightOfRefCSOriginMeter = 1400; // from Offline xml
240  double ObservationLevelMeter = 1300; // CORSIKA observation level [m]
241 
242  unsigned int countInFile = 0;
243 
244 
245  // event loop
246  while (dataFile.ReadNextEvent() == RecEventFile::eSuccess){
247 
248  ++countEvent;
249  ++countInFile;
250 
251  //if (EventID_int!=58)
252  //continue;
253 
254 
255  // compression is an experimental feature...
256  if (doCompress) {
257  (*theRecEventOut) = *(new RecEvent()); // reset
258  theRecEventOut->GetGenShower() = theRecEvent->GetGenShower();
259  theRecEventOut->SetDetector(theRecEvent->GetDetector());
260  theRecEventOut->SetEventId(theRecEvent->GetEventId());
261  vector<double> empty;
262  theRecEventOut->GetGenShower().SetDepth(empty);
263  theRecEventOut->GetGenShower().SetElectrons(empty);
264  theRecEventOut->GetGenShower().SetEnergyDeposit(empty);
265  }// end compression
266  else {
267  (*theRecEventOut) = (*theRecEvent);
268  }
269 
270  TVector3 shower_axis = theRecEvent->GetGenShower().GetAxisSiteCS();
271 
272  // this return core at relative height in reference-CS (ePampaAmarilla)
273  TVector3 shower_core = theRecEvent->GetGenShower().GetCoreAtAltitudeSiteCS(ObservationLevelMeter-HeightOfRefCSOriginMeter);
274 
275  double genX = shower_core.X();
276  double genY = shower_core.Y();
277  double genAzimuth = shower_axis.Phi() * degree;
278  if (genAzimuth<0) genAzimuth += 360;
279  double genZenith = shower_axis.Theta() * degree;
280  double genEnergy = theRecEvent->GetGenShower().GetEnergy();
281  double genXmax = theRecEvent->GetGenShower().GetXmaxInterpolated();
282 
283 
284  // find hottest Sim tel over all eyes
285 
286  int hottestEye = -1;
287  int hottestTel = -1;
288  double hottestTelPhotons = -1;
289  for (RecEvent::EyeIterator eyeIt= theRecEvent->EyesBegin();
290  eyeIt !=theRecEvent->EyesEnd();
291  ++eyeIt){
292  const int eyeId = eyeIt->GetEyeId();
293  //if (eyeId != 5) // only heat
294  //continue;
295  for (FDEvent::TelescopeDataConstIterator telIt = eyeIt->TelescopeDataBegin();
296  telIt != eyeIt->TelescopeDataEnd(); ++telIt) {
297  const FdTelescopeData& telData = telIt->second;
298  const FdApertureLight& theGenlight = telData.GetGenApertureLight();
299  const vector<double>& vecTotal = theGenlight.GetTotalLightAtAperture();
300  double totalPhotonsTel = 0;
301  for (const auto& itot : vecTotal) {totalPhotonsTel += itot;}
302  if (totalPhotonsTel > hottestTelPhotons) {
303  hottestEye = eyeId;
304  hottestTel = telIt->first;
305  hottestTelPhotons = totalPhotonsTel;
306  }
307  }
308  }
309 
310 
311  EFdRecLevel fdRec = eNoFdEvent;
312  const int eyeId = hottestEye;
313  bool isSel = false;
314  bool isHeatFOV = false;
315  //bool isStrange = false;
316  bool isTrig = false;
317  bool isReco = false;
318  double genRp = 0;
319  double genCherPhotons = 0;
320  double genTotalPhotons = 0;
321  double genXmaxAngle = -1;
322  double genMeanAngle = -1;
323  double genMinAngle = -1;
324  double genCherFrac = 0;
325 
326  double genChi0 = 0;
327 
328  double recEnergy = 0;
329  double recXmax = 0;
330  double recChi0 = 0;
331  double recAngle = 0;
332 
333  if (hottestTelPhotons>0) { // no HEAT... ?
334 
335  const FDEvent& eye = theRecEvent->GetEye(hottestEye);
336  const FdTelescopeData& telData = eye.GetTelescopeData(hottestTel);
337 
338  const int eyeId = eye.GetEyeId();
339  fdRec = eye.GetRecLevel();
340 
341  isSel = false;
342  isHeatFOV = false;
343  //isStrange = false;
344 
345  isTrig = (eye.HasT3Reconstruction());
346  isReco = (fdRec>=eHasEnergy);
347  genRp = eye.GetGenGeometry().GetRp();
348  genChi0 = eye.GetGenGeometry().GetChi0() * degree;
349 
350  const FdApertureLight& theGenlight = telData.GetGenApertureLight();
351  const vector<double>& vecCher = theGenlight.GetCherLightAtAperture();
352  const vector<double>& vecTotal = theGenlight.GetTotalLightAtAperture();
353  for (const auto& icher : vecCher) {genCherPhotons += icher;}
354  for (const auto& itot : vecTotal) {genTotalPhotons += itot;}
355 
356  genXmaxAngle = theGenlight.GetXmaxAngle() * degree;
357  genMeanAngle = theGenlight.GetMeanAngle() * degree;
358  genMinAngle = theGenlight.GetMinAngle() * degree;
359 
360  if (genTotalPhotons>100 && genMeanAngle<40 && eyeId==5) {
361  isHeatFOV = true;
362  }
363 
364  if (genTotalPhotons>0)
365  genCherFrac = genCherPhotons / genTotalPhotons;
366 
367  if (isReco) {
368  const FdRecShower& theRecShw = eye.GetFdRecShower();
369  const FdRecGeometry& theRecGeo = eye.GetFdRecGeometry();
370  recEnergy = theRecShw.GetEnergy();
371  recXmax = theRecShw.GetXmax();
372  recChi0 = theRecGeo.GetChi0() * degree;
373  const FdApertureLight& theRecLight = telData.GetRecApertureLight();
374  recAngle = theRecLight.GetXmaxAngle() * degree;
375  }
376 
377  //if (thisMeanAngle<genMeanAngle) {// && thisMeanAngle>0) {
378  //if (thisXmaxAngle<genXmaxAngle && thisXmaxAngle>0) {
379  //if (isTrig && genCherFrac > 0.5) {
380  // isSel = true;
381  //}
382 
383 
384  // compression is an experimental feature !
385  if (doCompress) {
386  if (eyeId == 5) {
387 
388  theRecEventOut->AddEye(eye);
389 
390  if (!isTrig) {
391  FDEvent& compressedEye = theRecEventOut->GetEye(5);
392  compressedEye.ClearStationList();
393  for (FDEvent::TelescopeDataIterator telIt = compressedEye.TelescopeDataBegin();
394  telIt != compressedEye.TelescopeDataEnd(); ++telIt) {
395  FdGenApertureLight& genLight = telIt->second.GetGenApertureLight();
396  const int size = genLight.GetTime().size();
397  if (size < 1)
398  continue;
399  const double deltaTsimMC = genLight.GetTimeBinning();
400  const double deltaTsimMCnew = deltaTsimMC * size;
401  genLight.SetTimeBinning(deltaTsimMCnew);
402 
403  {
404  vector<double> vnew;
405  vnew.push_back(genLight.GetTime()[0]);
406  genLight.SetTime(vnew);
407  }
408 
409  {
410  vector<double> vnew(1, 0.0);
411  for (const auto& ival : genLight.GetTotalLightAtAperture()) vnew[0] += ival;
412  genLight.SetTotalLightAtAperture(vnew);
413  }
414 
415  {
416  vector<double> vnew(1, 0.0);
417  for (const auto& ival : genLight.GetFluorLightAtAperture()) vnew[0] += ival;
418  genLight.SetFluorLightAtAperture(vnew);
419  }
420 
421  {
422  vector<double> vnew(1, 0.0);
423  for (const auto& ival : genLight.GetCherLightAtAperture()) vnew[0] += ival;
424  genLight.SetCherLightAtAperture(vnew);
425  }
426 
427  {
428  vector<double> vnew(1, 0.0);
429  for (const auto& ival : genLight.GetMieCherLightAtAperture()) vnew[0] += ival;
430  genLight.SetMieCherLightAtAperture(vnew);
431  }
432 
433  {
434  vector<double> vnew(1, 0.0);
435  for (const auto& ival : genLight.GetRayCherLightAtAperture()) vnew[0] += ival;
436  genLight.SetRayCherLightAtAperture(vnew);
437  }
438 
439  /*
440  if (HasMultipleScatteredLight()) {
441  const std::vector<Double_t>& GetMultScatFluorLightAtAperture() const { return fDia_multScat_fluor; }
442  const std::vector<Double_t>& GetMultScatCherLightAtAperture() const { return fDia_multScat_cher; }
443  }
444  */
445  //void SetMultScatFluorLightAtAperture(const std::vector<Double_t>& phot);
446  //void SetMultScatCherLightAtAperture(const std::vector<Double_t>& phot);
447  }
448  }
449  }
450  if (isHeatFOV) {
451  theRecEventOut->SetSDEvent(theRecEvent->GetSDEvent());
452  }
453  } // end compress!
454  }
455 
456  /*isStrange = log10(genTotalPhotons)<0 && isTrig;
457  if (isStrange) {
458  }*/
459 
460  //fill a text file with shower info to fed to CORSIKA steering cards. For showers with cherenkov fraction more than 50%
461  //if (genXmaxAngle<=15 && cherFrac>0.5){
462 
463  //if (genMeanAngle<=30 && genTotalPhotons>1){
464  if (genCherFrac>0.5 && isTrig && eyeId == 5) {
465 
466  isSel = true;
467 
468  ++countSelected;
469 
470  if (isReco)
471  ++countSelectedAndReco;
472 
473  preselADST.WriteEvent();
474 
475  } // selection
476  else {
477  if (genCherFrac>0.5)
478  ++countMiss;
479  }
480 
481  // ********************************************
482  // PLOTS
483 
484  if (isReco) { // add quality checks, too
485  proXmax->Fill(genCherFrac, recXmax-genXmax);
486  proEnergy->Fill(genCherFrac, recEnergy/genEnergy);
487  proChi0->Fill(genCherFrac, recChi0-genChi0);
488  proAngle->Fill(genCherFrac, recAngle-genXmaxAngle);
489  }
490 
491 
492  hRecLevel->Fill(fdRec);
493 
494  proCherXmax->Fill(genXmaxAngle, genCherFrac);
495  proCherMean->Fill(genMeanAngle, genCherFrac);
496 
497  hChFractXmax->Fill(genXmaxAngle, genCherFrac);
498  hChFractMin->Fill(genMinAngle, genCherFrac);
499  hChFractMean->Fill(genMeanAngle, genCherFrac);
500 
501  hEnergyGen->Fill(log10(genEnergy));
502  hCFracGen->Fill(genCherFrac);
503  hXmaxGen->Fill(genXmax);
504  hRpGen->Fill(genRp);
505  hZenithGen->Fill(genZenith);
506  hAngleGen->Fill(genXmaxAngle);
507  hLightTotalGen->Fill(log10(genTotalPhotons));
508  hLightCherGen->Fill(log10(genCherPhotons));
509  hPosGen->Fill(genX, genY);
510  hDirGen->Fill(genAzimuth, genZenith);
511 
512  if (isTrig) {
513  hEnergyTrg->Fill(log10(genEnergy));
514  hCFracTrg->Fill(genCherFrac);
515  hXmaxTrg->Fill(genXmax);
516  hRpTrg->Fill(genRp);
517  hZenithTrg->Fill(genZenith);
518  hAngleTrg->Fill(genXmaxAngle);
519  hLightTotalTrg->Fill(log10(genTotalPhotons));
520  hLightCherTrg->Fill(log10(genCherPhotons));
521  hPosTrg->Fill(genX, genY);
522  hDirTrg->Fill(genAzimuth, genZenith);
523  }
524 
525  if (isReco) {
526  hEnergyRec->Fill(log10(genEnergy));
527  hCFracRec->Fill(genCherFrac);
528  hXmaxRec->Fill(genXmax);
529  hRpRec->Fill(genRp);
530  hZenithRec->Fill(genZenith);
531  hAngleRec->Fill(genXmaxAngle);
532  hLightTotalRec->Fill(log10(genTotalPhotons));
533  hLightCherRec->Fill(log10(genCherPhotons));
534  hPosRec->Fill(genX, genY);
535  hDirRec->Fill(genAzimuth, genZenith);
536  }
537 
538  if (isSel) {
539  //if (isStrange) {
540  hEnergySel->Fill(log10(genEnergy));
541  hCFracSel->Fill(genCherFrac);
542  hXmaxSel->Fill(genXmax);
543  hRpSel->Fill(genRp);
544  hZenithSel->Fill(genZenith);
545  hAngleSel->Fill(genXmaxAngle);
546  hLightTotalSel->Fill(log10(genTotalPhotons));
547  hLightCherSel->Fill(log10(genCherPhotons));
548  hPosSel->Fill(genX, genY);
549  hDirSel->Fill(genAzimuth, genZenith);
550  }
551 
552  if (isHeatFOV) {
553  countHeatFOV++;
554  }
555 
556  } // events loop
557 
558  preselADST.Close();
559 
560  } // file loop
561 
562  // ********************************************
563  /*
564  hEnergyTrg->Divide(hEnergyGen);
565  hXmaxTrg->Divide(hXmaxGen);
566  hRpTrg->Divide(hRpGen);
567  hZenithTrg->Divide(hZenithGen);
568  hAngleTrg->Divide(hAngleGen);
569  hLightTotalTrg->Divide(hLightTotalGen);
570  hLightCherTrg->Divide(hLightCherGen);
571  hPosTrg->Divide(hPosGen);
572  hDirTrg->Divide(hDirGen);
573  hEnergyRec->Divide(hEnergyGen);
574  hXmaxRec->Divide(hXmaxGen);
575  hRpRec->Divide(hRpGen);
576  hZenithRec->Divide(hZenithGen);
577  hAngleRec->Divide(hAngleGen);
578  hLightTotalRec->Divide(hLightTotalGen);
579  hLightCherRec->Divide(hLightCherGen);
580  hPosRec->Divide(hPosGen);
581  hDirRec->Divide(hDirGen);
582  */
583 
584  outPlots->Write();
585  outPlots->Close();
586 
587  cout << "Finished: countFile=" << countFile << " countEvent=" << countEvent
588  << " countSelected=" << countSelected
589  << " countSelectedAndReco=" << countSelectedAndReco
590  << " countHeatFOV=" << countHeatFOV
591  << " countMiss=" << countMiss << " (w/ chFrac>0.5)"
592  << endl;
593 
594  boost::filesystem::remove( steerdir / "preselect.txt-LOCK");
595  cout << "unlocked: " << (steerdir / "preselect.txt-LOCK").string() << endl;
596 }
597 
598 
599 
600 
601 int
602 main(int argc, char** argv)
603 {
604  if (argc<3) {
605  cout << " specify steerdir AND input file" << endl;
606  return 1;
607  }
608  const string steerdir(argv[1]);
609  boost::filesystem::path dir(steerdir);
610  if (!boost::filesystem::is_directory(dir)) {
611  cout << " specify steerdir AND input file" << endl;
612  cout << dir << " is not a directory" << endl;
613  return 1;
614  }
615 
616  vector<string> inputs;
617  for (int i=1; i<argc-1; ++i)
618  inputs.push_back(string(argv[i+1]));
619  preselect(steerdir, inputs);
620  return 0;
621 }
622 
degree
const double degree
Definition: PhysicsConstants.h:3
io::eWrite
Read Only access.
Definition: IoCodes.h:18
exit
int exit
Definition: dump1090.h:237
doCompress
bool doCompress
Definition: preselect.cc:50
exists
char * exists
Definition: XbArray.cc:12
main
int main(int argc, char *argv[])
Definition: DBSync.cc:58
preselect
void preselect(boost::filesystem::path steerdir, const vector< string > &inputs)
Definition: preselect.cc:53
, generated on Tue Sep 26 2023.