4 #include <unordered_map> 60 if(run < other.
run )
return true;
62 if(run == other.
run &&
63 subrun < other.
subrun)
return true;
65 if(run == other.
run &&
67 event < other.
event )
return true;
69 if(run == other.
run &&
71 event == other.
event &&
72 slice < other.
slice )
return true;
74 if(run == other.
run &&
76 event == other.
event &&
77 slice == other.
slice &&
78 cycle < other.
cycle )
return true;
90 std::set<EventID> & eventList,
97 std::ifstream eventListFile(inputName);
101 << eventListFile.is_open()
109 int highEventNumCnt = 0;
113 if (selection ==
"nue") selectToUse =
"NuE";
114 else if(selection ==
"nus") selectToUse =
"NC";
116 while( !eventListFile.eof() && eventListFile.is_open() ){
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;
123 eventListFile >> run >> subrun >>
event >> slice >> cycle >> energy >> cvn >>
wgt;
124 select =
"caf_" + selection;
144 if(eventList.size() < 10)
162 eventList.insert(evId);
172 <<
" high event numbers" 180 std::set<EventID>
const& firstList,
181 std::set<EventID>
const& secondList,
185 std::unordered_map<std::string, TH1D*> & weightedSpect,
186 std::unordered_map<std::string, TH1D*> & unweightedSpect)
190 std::ofstream missingList(missingListName);
192 for(
auto const& itr : firstList){
193 if(secondList.find(itr) == secondList.end()){
196 missingList << itr.run
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);
229 std::unordered_map<std::string, TH1D*> weightedSpect,
230 std::unordered_map<std::string, TH1D*> unweightedSpect,
231 std::unordered_map<std::string, TH1D*> spectRatio)
244 if(fileName.find(
"nus") != std::string::npos){
251 TCanvas *
canv =
new TCanvas(
"compareCanv",
252 "Compare CMF vs CAF",
258 canv->Divide(xPads, yPads);
260 canv->cd(countPad)->SetLogy();
263 canv->cd(deltaEPad)->SetLogy();
264 energyDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
267 canv->cd(deltaWgtPad)->SetLogy();
268 wgtDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
271 if(fileName.find(
"nue") != std::string::npos){
272 canv->cd(deltaCVNPad)->SetLogy();
273 cvnDiff->GetXaxis()->SetRangeUser(-0.5, 0.5);
278 for(
auto const& itr : weightedSpect){
279 canv->cd(i + spectPadBeg);
281 itr.second->SetMaximum(
std::max(itr.second->GetMaximum(), unweightedSpect[itr.first]->GetMaximum()) * 1.75);
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);
292 canv->cd(i + ratioPadBeg)->SetGridy();
293 spectRatio[itr.first]->GetYaxis()->SetRangeUser(0.5, 1.5);
294 spectRatio[itr.first]->Draw();
301 canv->Print(fileName.c_str(),
"pdf");
307 for(
auto const& itr : weightedSpect){
308 itr.second ->Write();
309 unweightedSpect[itr.first]->Write();
310 spectRatio [itr.first]->Write();
323 if( !syst.empty() ) syst.insert(0,
"_");
325 std::set<EventID> eventListEvents;
326 std::set<EventID> cafEvents;
328 TFile *
file =
new TFile((
"compare_hist_" +
det +
"_" + sel +
"_" + horncur +
"_" + sample + syst +
"_" +
tag +
".root").c_str(),
"RECREATE");
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();
334 TH1D* cvnDiff =
new TH1D(
"cvnDiff",
";#Delta CVN;Events", 2000, -1.
e-4, 1.
e-4);
335 cvnDiff->GetXaxis()->CenterTitle();
336 cvnDiff->GetYaxis()->CenterTitle();
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();
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();
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);
360 sels.emplace_back(
"NuMuSelQ1");
361 sels.emplace_back(
"NuMuSelQ2");
362 sels.emplace_back(
"NuMuSelQ3");
363 sels.emplace_back(
"NuMuSelQ4");
365 else if(sel ==
"nue"){
367 sels.emplace_back(
"NuESel_LowPID");
368 sels.emplace_back(
"NuESel_HighPID");
369 sels.emplace_back(
"NuESel_Peripheral");
371 else if(sel ==
"nus"){
373 sels.emplace_back(
"NCSel");
376 std::unordered_map<std::string, TH1D*> weightedSpect;
377 std::unordered_map<std::string, TH1D*> unweightedSpect;
378 std::unordered_map<std::string, TH1D*> spectRatio;
380 for(
auto const& itr : sels){
383 weightedSpect[itr] =
new TH1D((
"weighted" + itr).c_str(),
384 ";E (GeV);Events / GeV (#times 1000)",
387 weightedSpect[itr]->GetXaxis()->CenterTitle();
388 weightedSpect[itr]->GetYaxis()->CenterTitle();
389 weightedSpect[itr]->GetXaxis()->SetDecimals();
390 weightedSpect[itr]->GetYaxis()->SetDecimals();
391 weightedSpect[itr]->Sumw2();
393 unweightedSpect[itr] =
new TH1D((
"unweighted" + itr).c_str(),
394 ";E (GeV);Events / GeV (#times 1000)",
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);
405 spectRatio[itr] =
new TH1D((
"spectRatio" + itr).c_str(),
406 ";E (GeV);Weighted / Unweighted",
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);
417 std::string inputName(
"selectedEvents_" +
det +
"_" + sel +
"_" + horncur +
"_" + sample + syst +
"_" +
tag +
".txt");
423 inputName =
"cafEvents_" +
det +
"_" + sel +
"_" + horncur +
"_" + sample + syst +
"_" +
tag +
".txt";
425 size_t pos = inputName.find(
"-");
426 while( pos != std::string::npos){
427 inputName.replace(pos, 1,
"_");
428 pos = inputName.find(
"-");
435 "missing_caf_list_" +
det +
"_" + sel +
"_" + horncur +
"_" + sample + syst +
"_" + tag +
".txt",
445 "missing_event_list_" +
det +
"_" + sel +
"_" + horncur +
"_" + sample + syst +
"_" + tag +
".txt",
456 <<
" event list events not in the caf events and " 457 << cntMissingEventList
458 <<
" caf events not in the event list events" 463 for(
auto const& itr : weightedSpect){
464 spectRatio[itr.first]->Divide(weightedSpect[itr.first], unweightedSpect[itr.first]);
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);
472 unweightedSpect[itr.first]->SetBinContent(
i, unweightedSpect[itr.first]->GetBinContent(
i) *
norm);
473 unweightedSpect[itr.first]->SetBinError(
i, unweightedSpect[itr.first]->GetBinError(
i) *
norm);
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);
482 drawHistograms(
"compare_" +
det +
"_" + sel +
"_" + horncur +
"_" + sample +
"_" + tag +
".pdf",
T max(const caf::Proxy< T > &a, T b)
dictionary missing
Analysis.
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)
bool operator<(EventID const &other) const
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)