compareEvents.C
Go to the documentation of this file.
1 #include <string>
2 #include <set>
3 #include <list>
4 #include <unordered_map>
5 #include <iostream>
6 #include <fstream>
7 #include <iomanip>
8 #include <math.h>
9 
10 #include "TH1.h"
11 #include "TCanvas.h"
12 #include "TFile.h"
13 #include "TLegend.h"
14 
15 struct EventID{
16 
18  : run(0)
19  , subrun(0)
20  , event(0)
21  , slice(0)
22  , cycle(0)
23  , energy(0.)
24  , cvn(0.)
25  , wgt(0.)
26  , select("")
27  {}
28 
29  EventID(long r,
30  long sr,
31  long ev,
32  long sl,
33  long c,
34  double en,
35  double id,
36  double weight,
37  std::string sel)
38  : run(r)
39  , subrun(sr)
40  , event(ev)
41  , slice(sl)
42  , cycle(c)
43  , energy(en)
44  , cvn(id)
45  , wgt(weight)
46  , select(sel)
47  {}
48 
49  long run ;
50  long subrun;
51  long event ;
52  long slice ;
53  long cycle ;
54  double energy;
55  double cvn;
56  double wgt;
58 
59  bool operator<(EventID const& other) const {
60  if(run < other.run ) return true;
61  else{
62  if(run == other.run &&
63  subrun < other.subrun) return true;
64  else{
65  if(run == other.run &&
66  subrun == other.subrun &&
67  event < other.event ) return true;
68  else {
69  if(run == other.run &&
70  subrun == other.subrun &&
71  event == other.event &&
72  slice < other.slice ) return true;
73  else{
74  if(run == other.run &&
75  subrun == other.subrun &&
76  event == other.event &&
77  slice == other.slice &&
78  cycle < other.cycle ) return true;
79  }
80  }
81  }
82  }
83  return false;
84  }
85 
86 };
87 
88 //-----------------------------------------------------------------------------
89 void fillEventIDsFromFile(std::string const& inputName,
90  std::set<EventID> & eventList,
91  std::string const& selection)
92 {
93  std::cout << inputName
94  << std::endl;
95 
96  // read in the event lists list of events
97  std::ifstream eventListFile(inputName);
98 
99  std::cout << inputName
100  << " is open "
101  << eventListFile.is_open()
102  << std::endl;
103 
104  double energy, cvn, wgt;
105  double run, subrun, event, slice, cycle;
106 
107  long ev;
108 
109  int highEventNumCnt = 0;
110 
112  std::string selectToUse = "NuMu";
113  if (selection == "nue") selectToUse = "NuE";
114  else if(selection == "nus") selectToUse = "NC";
115 
116  while( !eventListFile.eof() && eventListFile.is_open() ){
117 
118  if(inputName.find("caf") == std::string::npos){
119  eventListFile >> select >> run >> subrun >> event >> slice >> cycle >> energy >> cvn >> wgt;
120  if(select.find(selectToUse) == std::string::npos) continue;
121  }
122  else{
123  eventListFile >> run >> subrun >> event >> slice >> cycle >> energy >> cvn >> wgt;
124  select = "caf_" + selection;
125  }
126 
127  if(event > 1.e6){
128  ev = event / 10;
129  ev *= 10.;
130  ++highEventNumCnt;
131  }
132  else ev = event;
133 
134  EventID evId((long)run,
135  (long)subrun,
136  (long)ev,
137  (long)slice,
138  (long)cycle,
139  energy,
140  cvn,
141  wgt,
142  select);
143 
144  if(eventList.size() < 10)
145  std::cout << select
146  << " "
147  << evId.run
148  << " "
149  << evId.subrun
150  << " "
151  << evId.event
152  << " / "
153  << event
154  << " "
155  << evId.slice
156  << " "
157  << evId.cycle
158  << " "
159  << eventList.size()
160  << std::endl;
161 
162  eventList.insert(evId);
163 
164  } // end loop to load event lists events
165 
166  std::cout << "There are "
167  << eventList.size()
168  << " events from "
169  << inputName
170  << " and "
171  << highEventNumCnt
172  << " high event numbers"
173  << std::endl;
174 
175 }
176 
177 //-----------------------------------------------------------------------------
178 int findMissingEvents(bool firstCAF,
179  std::string const& missingListName,
180  std::set<EventID> const& firstList,
181  std::set<EventID> const& secondList,
182  TH1D* energyDiff,
183  TH1D* cvnDiff,
184  TH1D* wgtDiff,
185  std::unordered_map<std::string, TH1D*> & weightedSpect,
186  std::unordered_map<std::string, TH1D*> & unweightedSpect)
187 {
188  int cntMissing = 0;
189 
190  std::ofstream missingList(missingListName);
191 
192  for(auto const& itr : firstList){
193  if(secondList.find(itr) == secondList.end()){
194  ++cntMissing;
195 
196  missingList << itr.run
197  << " "
198  << itr.subrun
199  << " "
200  << itr.event
201  << " "
202  << itr.slice
203  << " "
204  << itr.cycle
205  << std::endl;
206  }
207  else{
208  if(!firstCAF){
209  energyDiff ->Fill(itr.energy - secondList.find(itr)->energy);
210  cvnDiff ->Fill(itr.cvn - secondList.find(itr)->cvn);
211  wgtDiff ->Fill(itr.wgt - secondList.find(itr)->wgt);
212  weightedSpect[itr.select] ->Fill(itr.energy, itr.wgt * 1.e-3);
213  unweightedSpect[itr.select]->Fill(itr.energy, 1.e-3);
214  }
215  } // end we found the event in the second list
216  } // end loop to look for missing caf evnets
217 
218  missingList.close();
219 
220  return cntMissing;
221 }
222 
223 //-----------------------------------------------------------------------------
225  TH1D* energyDiff,
226  TH1D* cvnDiff,
227  TH1D* wgtDiff,
228  TH1D* missing,
229  std::unordered_map<std::string, TH1D*> weightedSpect,
230  std::unordered_map<std::string, TH1D*> unweightedSpect,
231  std::unordered_map<std::string, TH1D*> spectRatio)
232 {
233  // set how many pads and which pad is which
234  // plot now
235  int xPads = 4;
236  int yPads = 3;
237  int countPad = 1;
238  int deltaEPad = 2;
239  int deltaWgtPad = 3;
240  int deltaCVNPad = 4;
241  int spectPadBeg = 5;
242  int ratioPadBeg = 9;
243 
244  if(fileName.find("nus") != std::string::npos){
245  xPads = 3;
246  yPads = 2;
247  spectPadBeg = 4;
248  ratioPadBeg = 5;
249  }
250 
251  TCanvas *canv = new TCanvas("compareCanv",
252  "Compare CMF vs CAF",
253  300 * xPads,
254  300 * yPads);
255 
256  // the number of divisions is equal to
257  // the number of selections in the x and 4 in y
258  canv->Divide(xPads, yPads);
259 
260  canv->cd(countPad)->SetLogy();
261  missing->Draw();
262 
263  canv->cd(deltaEPad)->SetLogy();
264  energyDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
265  energyDiff->Draw();
266 
267  canv->cd(deltaWgtPad)->SetLogy();
268  wgtDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
269  wgtDiff->Draw();
270 
271  if(fileName.find("nue") != std::string::npos){
272  canv->cd(deltaCVNPad)->SetLogy();
273  cvnDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
274  cvnDiff->Draw();
275  }
276 
277  size_t i = 0;
278  for(auto const& itr : weightedSpect){
279  canv->cd(i + spectPadBeg);
280 
281  itr.second->SetMaximum(std::max(itr.second->GetMaximum(), unweightedSpect[itr.first]->GetMaximum()) * 1.75);
282  itr.second->Draw();
283  unweightedSpect[itr.first]->Draw("same");
284  TLegend *leg = new TLegend(0.5, 0.5, 0.8, 0.8);
285  leg->AddEntry(itr.second, itr.first.c_str(), "");
286  leg->AddEntry(itr.second, "Weighted Spectrum", "l");
287  leg->AddEntry(unweightedSpect[itr.first], "Unweighted Spectrum", "l");
288  leg->SetFillStyle(0);
289  leg->SetBorderSize(0);
290  leg->Draw();
291 
292  canv->cd(i + ratioPadBeg)->SetGridy();
293  spectRatio[itr.first]->GetYaxis()->SetRangeUser(0.5, 1.5);
294  spectRatio[itr.first]->Draw();
295 
296  ++i;
297  } // end loop to draw spectra
298 
299  canv->Update();
300  canv->Write();
301  canv->Print(fileName.c_str(), "pdf");
302 
303  energyDiff->Write();
304  cvnDiff ->Write();
305  missing ->Write();
306 
307  for(auto const& itr : weightedSpect){
308  itr.second ->Write();
309  unweightedSpect[itr.first]->Write();
310  spectRatio [itr.first]->Write();
311  }
312 
313 }
314 
315 //-----------------------------------------------------------------------------
316 void compareEvents(std::string det="DETECTOR",
317  std::string sel="SELECT",
318  std::string sample="FILETYPE",
319  std::string horncur="HORNCUR",
320  std::string tag="TAG",
321  std::string syst="SYST")
322 {
323  if( !syst.empty() ) syst.insert(0, "_");
324 
325  std::set<EventID> eventListEvents;
326  std::set<EventID> cafEvents;
327 
328  TFile *file = new TFile(("compare_hist_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".root").c_str(), "RECREATE");
329 
330  TH1D* energyDiff = new TH1D("energyDiff", ";#Delta E (GeV);Events", 2000, -1.e-4, 1.e-4);
331  energyDiff->GetXaxis()->CenterTitle();
332  energyDiff->GetYaxis()->CenterTitle();
333 
334  TH1D* cvnDiff = new TH1D("cvnDiff", ";#Delta CVN;Events", 2000, -1.e-4, 1.e-4);
335  cvnDiff->GetXaxis()->CenterTitle();
336  cvnDiff->GetYaxis()->CenterTitle();
337 
338  TH1D* wgtDiff = new TH1D("wgtDiff", ";#Delta (PPFX #times XSecCV);Events", 2000, -1.e-4, 1.e-4);
339  wgtDiff->GetXaxis()->CenterTitle();
340  wgtDiff->GetYaxis()->CenterTitle();
341 
342  // Histogram to visualize the number of missing events and totals
343  TH1D* missing = new TH1D("missing", ";Missing in CAF or CMF;Events", 4, 0., 4.);
344  missing->GetXaxis()->SetBinLabel(1, "Total CMF");
345  missing->GetXaxis()->SetBinLabel(2, "CMF not in CAF");
346  missing->GetXaxis()->SetBinLabel(3, "Total CAF");
347  missing->GetXaxis()->SetBinLabel(4, "CAF not in CMF");
348  missing->GetXaxis()->CenterTitle();
349  missing->GetYaxis()->CenterTitle();
350 
351  std::list<std::string> sels;
352  std::vector<double> bins;
353  std::vector<double> nueBins ({1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5.});
354  std::vector<double> numuBins({0., 0.75, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 5.});
355  std::vector<double> nusBins;
356  for(double b = 0.75; b < 20.25; b += 0.25) nusBins.push_back(b);
357 
358  if(sel == "numu"){
359  bins.swap(numuBins);
360  sels.emplace_back("NuMuSelQ1");
361  sels.emplace_back("NuMuSelQ2");
362  sels.emplace_back("NuMuSelQ3");
363  sels.emplace_back("NuMuSelQ4");
364  }
365  else if(sel == "nue"){
366  bins.swap(nueBins);
367  sels.emplace_back("NuESel_LowPID");
368  sels.emplace_back("NuESel_HighPID");
369  sels.emplace_back("NuESel_Peripheral");
370  }
371  else if(sel == "nus"){
372  bins.swap(nusBins);
373  sels.emplace_back("NCSel");
374  }
375 
376  std::unordered_map<std::string, TH1D*> weightedSpect;
377  std::unordered_map<std::string, TH1D*> unweightedSpect;
378  std::unordered_map<std::string, TH1D*> spectRatio;
379 
380  for(auto const& itr : sels){
381  //std::cout << itr << std::endl;
382 
383  weightedSpect[itr] = new TH1D(("weighted" + itr).c_str(),
384  ";E (GeV);Events / GeV (#times 1000)",
385  bins.size() - 1,
386  &bins[0]);
387  weightedSpect[itr]->GetXaxis()->CenterTitle();
388  weightedSpect[itr]->GetYaxis()->CenterTitle();
389  weightedSpect[itr]->GetXaxis()->SetDecimals();
390  weightedSpect[itr]->GetYaxis()->SetDecimals();
391  weightedSpect[itr]->Sumw2();
392 
393  unweightedSpect[itr] = new TH1D(("unweighted" + itr).c_str(),
394  ";E (GeV);Events / GeV (#times 1000)",
395  bins.size() - 1,
396  &bins[0]);
397  unweightedSpect[itr]->GetXaxis()->CenterTitle();
398  unweightedSpect[itr]->GetYaxis()->CenterTitle();
399  unweightedSpect[itr]->GetXaxis()->SetDecimals();
400  unweightedSpect[itr]->GetYaxis()->SetDecimals();
401  unweightedSpect[itr]->Sumw2();
402  unweightedSpect[itr]->SetLineColor(2);
403  unweightedSpect[itr]->SetMarkerColor(2);
404 
405  spectRatio[itr] = new TH1D(("spectRatio" + itr).c_str(),
406  ";E (GeV);Weighted / Unweighted",
407  bins.size() - 1,
408  &bins[0]);
409  spectRatio[itr]->GetXaxis()->CenterTitle();
410  spectRatio[itr]->GetYaxis()->CenterTitle();
411  spectRatio[itr]->GetXaxis()->SetDecimals();
412  spectRatio[itr]->GetYaxis()->SetDecimals();
413  spectRatio[itr]->SetMarkerStyle(22);
414  spectRatio[itr]->SetMarkerSize(0.8);
415  } // end loop to make spectra for each selection
416 
417  std::string inputName("selectedEvents_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt");
418 
419  fillEventIDsFromFile(inputName,
420  eventListEvents,
421  sel);
422 
423  inputName = "cafEvents_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt";
424 
425  size_t pos = inputName.find("-");
426  while( pos != std::string::npos){
427  inputName.replace(pos, 1, "_");
428  pos = inputName.find("-");
429  }
430 
431  fillEventIDsFromFile(inputName, cafEvents, sel);
432 
433  // loop over the event list and look for missing caf events
434  int cntMissingCAF = findMissingEvents(false,
435  "missing_caf_list_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt",
436  eventListEvents,
437  cafEvents,
438  energyDiff,
439  cvnDiff,
440  wgtDiff,
441  weightedSpect,
442  unweightedSpect);
443 
444  int cntMissingEventList = findMissingEvents(true,
445  "missing_event_list_" + det + "_" + sel + "_" + horncur + "_" + sample + syst + "_" + tag + ".txt",
446  cafEvents,
447  eventListEvents,
448  energyDiff,
449  cvnDiff,
450  wgtDiff,
451  weightedSpect,
452  unweightedSpect);
453 
454  std::cout << "There are "
455  << cntMissingCAF
456  << " event list events not in the caf events and "
457  << cntMissingEventList
458  << " caf events not in the event list events"
459  << std::endl;
460 
461  // make the ratios of the weighted and unweighted histograms
462  double norm;
463  for(auto const& itr : weightedSpect){
464  spectRatio[itr.first]->Divide(weightedSpect[itr.first], unweightedSpect[itr.first]);
465 
466  // normalize the bin contents to be per 0.1 GeV
467  for(int i = 1; i < itr.second->GetNbinsX() + 1; ++i){
468  norm = 1. / itr.second->GetBinWidth(i);
469  itr.second->SetBinContent(i, itr.second->GetBinContent(i) * norm);
470  itr.second->SetBinError(i, itr.second->GetBinError(i) * norm);
471 
472  unweightedSpect[itr.first]->SetBinContent(i, unweightedSpect[itr.first]->GetBinContent(i) * norm);
473  unweightedSpect[itr.first]->SetBinError(i, unweightedSpect[itr.first]->GetBinError(i) * norm);
474  }
475  }
476 
477  missing->Fill(0.5, eventListEvents.size());
478  missing->Fill(1.5, cntMissingCAF);
479  missing->Fill(2.5, cafEvents.size());
480  missing->Fill(3.5, cntMissingEventList);
481 
482  drawHistograms("compare_" + det + "_" + sel + "_" + horncur + "_" + sample + "_" + tag + ".pdf",
483  energyDiff,
484  cvnDiff,
485  wgtDiff,
486  missing,
487  weightedSpect,
488  unweightedSpect,
489  spectRatio);
490 
491  file->Write();
492  file->Close();
493 }
T max(const caf::Proxy< T > &a, T b)
fileName
Definition: plotROC.py:78
const Var weight
TCanvas * canv
double cvn
Definition: compareEvents.C:55
Defines an enumeration for prong classification.
void drawHistograms(std::string const &fileName, TH1D *energyDiff, TH1D *cvnDiff, TH1D *wgtDiff, TH1D *missing, std::unordered_map< std::string, TH1D * > weightedSpect, std::unordered_map< std::string, TH1D * > unweightedSpect, std::unordered_map< std::string, TH1D * > spectRatio)
void fillEventIDsFromFile(std::string const &inputName, std::set< EventID > &eventList, std::string const &selection)
Definition: compareEvents.C:89
double energy
Definition: compareEvents.C:54
const Cut sels[kNumSels]
Definition: vars.h:44
bool operator<(EventID const &other) const
Definition: compareEvents.C:59
caf::StandardRecord * sr
const Binning bins
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
Float_t norm
const hit & b
Definition: hits.cxx:21
TFile * file
Definition: cellShifts.C:17
TRandom3 r(0)
double wgt
Definition: compareEvents.C:56
Float_t e
Definition: plot.C:35
int findMissingEvents(bool firstCAF, std::string const &missingListName, std::set< EventID > const &firstList, std::set< EventID > const &secondList, TH1D *energyDiff, TH1D *cvnDiff, TH1D *wgtDiff, std::unordered_map< std::string, TH1D * > &weightedSpect, std::unordered_map< std::string, TH1D * > &unweightedSpect)
void compareEvents(std::string det="DETECTOR", std::string sel="SELECT", std::string sample="FILETYPE", std::string horncur="HORNCUR", std::string tag="TAG", std::string syst="SYST")
EventID(long r, long sr, long ev, long sl, long c, double en, double id, double weight, std::string sel)
Definition: compareEvents.C:29
std::string select
Definition: compareEvents.C:57
enum BeamMode string