13 #include <unordered_map> 14 #include <unordered_set> 16 #include "boost/algorithm/string.hpp" 17 #include "boost/tokenizer.hpp" 27 #include "TLegendEntry.h" 60 {
"MEC",
"NOvA 2p2h" },
78 const std::unordered_map<std::string, double>
Q0Q3_MAX = {
98 _combination(Form(
"{%s}{%s}", wgtnm.c_str(), intnm.c_str()))
136 : fHornCurrent(hornCurrent), fProdName(prodName),
137 fOutDir( Form(
"%s/%s/%s", outDir.c_str(), hornCurrent.c_str(), prodName.c_str()) ),
138 fHistFileData(TFile::Open( Form(
"%s/%s/rw_histos_%s_data.root", outDir.c_str(), hornCurrent.c_str(), fProdName.c_str()),
"recreate" )),
139 fHistFileMC(TFile::Open( Form(
"%s/%s/rw_histos_%s_MC.root", outDir.c_str(), hornCurrent.c_str(), fProdName.c_str()),
"recreate" ))
141 this->GetDataPlots();
147 fHistFileData->Close();
148 fHistFileMC->Close();
150 delete fHistFileData;
159 std::map<std::string, std::unique_ptr<ana::Spectrum>>
fDataPlots;
160 std::map<std::string, std::map<PlotKey,std::unique_ptr<ana::Spectrum>>>
fMCPlots;
172 auto name = plotName.c_str();
174 TFile *
f = (isData) ? fHistFileData : fHistFileMC;
175 decltype(fdataPlotsWritten) & cache = isData ? fdataPlotsWritten : fMCPlotsWritten;
177 if (cache.find(
name) != cache.end())
189 TFile
f( Form(
"/nova/ana/users/jwolcott/2p2h_tuning/2018/hists/rwgt_spectra_%s_%s_Data.root", fProdName.c_str(), fHornCurrent.c_str()) );
190 TList *
keys = f.GetListOfKeys();
194 std::for_each(it.Begin(),
TIter::End(), [
this](TObject * listItem){
195 TKey *
key =
dynamic_cast<TKey*
>(listItem);
197 TDirectory *
dir =
dynamic_cast<TDirectory*
>(key->ReadObj());
209 std::vector<std::string>
fields;
211 typedef boost::tokenizer<boost::char_separator<char> >
213 boost::char_separator<char>
sep(
"{}");
214 tokenizer
tokens(keyedName, sep);
217 for (
const auto &
f : tokens)
224 std::vector<std::string> subfields;
225 boost::split(subfields,
f, boost::is_any_of(
"="), boost::token_compress_on);
226 fields.push_back(subfields.back());
233 if (fields.size() == 2)
234 plotkey =
PlotKey(
"", fields[1]);
236 else if (fields.size() == 3)
237 plotkey =
PlotKey(fields[1], fields[2]);
240 std::cerr << Form(
"Warning: Can't parse Spectrum name: '%s'", keyedName.c_str()) <<
std::endl;
257 auto fileName = Form(
"/nova/ana/users/jwolcott/2p2h_tuning/2018/hists/rwgt_spectra_%s_%s_MC.root",
258 this->fProdName.c_str(), this->fHornCurrent.c_str());
263 TList *
keys = f.GetListOfKeys();
266 std::for_each(it.Begin(),
TIter::End(), [
this](TObject * listItem){
267 TKey *
key =
dynamic_cast<TKey*
>(listItem);
269 TDirectory *
dir =
dynamic_cast<TDirectory*
>(key->ReadObj());
275 bool success = SplitMCPlotName(key->GetName(),
name, plotkey);
280 if (fMCPlots[name].find(plotkey) == fMCPlots[
name].end() || fMCPlots[
name].at(plotkey)->ToTH1(fMCPlots[name].at(plotkey)->
POT())->GetEntries() == 0)
285 for (
const auto & plotPair : fMCPlots)
287 std::cout << plotPair.first <<
":\n";
310 THStack * stack =
new THStack;
311 stack->SetName(
"stack");
312 TLegend *
legend =
new TLegend(0.7, 0.5, 0.9, 0.85);
313 legend->SetFillStyle(0);
317 TCanvas *
canvas =
new TCanvas;
318 double dataPOT = dataPlot->
POT();
319 static bool announcedPOT =
false;
325 TH1 * dataTH1 = dataPlot->
ToTH1(dataPOT);
326 dataTH1->SetLineWidth(2);
327 dataTH1->SetMarkerStyle(8);
328 if (
std::string(dataTH1->GetXaxis()->GetTitle()).find(
"\"q0\"") != std::string::npos)
329 dataTH1->SetMarkerSize(1.2);
331 dataTH1->SetMarkerSize(0.8);
333 this->WritePlotToHistFile(plotName.c_str(), dataTH1,
true);
338 TH1 * mcTH1 = mcplots.at(key)->ToTH1(dataPOT);
339 mcTH1->SetLineWidth(1);
342 this->WritePlotToHistFile(plotName + key.
combination(), mcTH1,
false);
343 stack->Add(mcTH1,
"hist");
345 if (plotKey.
mecWgt().find(
"NoMEC") == std::string::npos || intClass !=
"MEC")
346 legend->GetListOfPrimitives()->AddFirst(
new TLegendEntry(mcTH1,
INT_LEGEND.at(intClass).c_str(),
"f"));
350 TH1 * mcsum =
dynamic_cast<TH1*
>((*stack->GetStack())[stack->GetStack()->LastIndex()]);
352 auto maxVal =
std::max(dataTH1->GetMaximum(), mcsum->GetMaximum());
355 dataTH1->SetMaximum(maxVal);
356 dataTH1->GetYaxis()->SetRange(0, maxVal);
358 dataTH1->Draw(
"pe same");
371 for (
const auto & dataPair : fDataPlots)
374 const std::unique_ptr<ana::Spectrum> & dataPlot = dataPair.second;
378 if (fMCPlots.find(plotName) == fMCPlots.end())
380 std::cerr << Form(
" warning: MC plot with mec wgt '%s' not found for data plot '%s'", mecWgtName.c_str(), plotName.c_str()) <<
std::endl;
386 const auto & plotCollection = fMCPlots.at(plotName);
394 std::map<PlotKey, const ana::Spectrum*> plotsByInteraction;
395 for (
const auto & plotPair : plotCollection)
398 if ( plotPair.first.mecWgt() != key.
mecWgt())
404 plotsByInteraction.insert(
std::make_pair(plotPair.first, plotPair.second.get()) );
410 std::unique_ptr<TCanvas> dataMCComparison(
DataMCComparison(plotName, key, dataPlot.get(), plotsByInteraction,
key));
412 std::size_t q3Loc = plotName.find(
"q3=");
413 if (q3Loc != std::string::npos)
415 TIter iter(dataMCComparison->GetListOfPrimitives());
417 while ( (obj = iter()) )
419 TH1 *
h =
dynamic_cast<TH1*
>(obj);
422 h->GetXaxis()->SetRangeUser(0, 1);
423 h->GetYaxis()->SetRangeUser(0,
Q0Q3_MAX.at(this->fHornCurrent));
426 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
427 dataMCComparison->Print( Form(
"%s/%s/%s.%s", fOutDir.c_str(), mecWgtName.c_str(), plotName.c_str(),
ext) );
430 if (q3Loc != std::string::npos)
432 auto lb = std::atof(plotName.substr(q3Loc+3, plotName.find(
"-")).c_str());
433 auto & canvasColl = this->fq0q3Canvases[
key];
434 int q3SliceIdx = 10*lb;
435 if (q3SliceIdx > (
int)(canvasColl.size())-1)
436 canvasColl.resize(q3SliceIdx+1);
438 canvasColl[q3SliceIdx] = std::move(dataMCComparison);
443 std::cout <<
" couldn't find plot '" << plotName <<
"' with key '" << key <<
"'. Skipping data-MC comparison." <<
std::endl;
454 if (this->fq0q3Canvases.empty())
460 for (
const auto & keyPair : this->fq0q3Canvases)
462 auto &
key = keyPair.first;
463 auto & canvases = keyPair.second;
466 const double bottomMarginMaster = 0.12;
467 const double leftMarginMaster = 0.1;
468 const double rightMarginMaster = 0.05;
469 const double topMarginMaster = 0.05;
471 TCanvas masterC(
"q0q3panel",
"(q0, |q|) panel", 700, 700);
472 auto nCanvases = canvases.size();
474 auto nRows = nCanvases / 3;
475 double padWidth = (1. - leftMarginMaster - rightMarginMaster) / 3;
476 double padHeight = (1. - topMarginMaster - bottomMarginMaster) / nRows;
487 TIter it_prim(canvases[cIdx]->GetListOfPrimitives());
488 bool adjustedData =
false;
489 while (TObject * obj = it_prim())
491 if (
auto l = dynamic_cast<TLegend*>(obj))
492 canvases[cIdx]->GetListOfPrimitives()->Remove(
l);
493 if (
auto h = dynamic_cast<TH1*>(obj) )
502 h->GetYaxis()->SetRangeUser(0, q0max);
503 h->GetYaxis()->SetDecimals();
507 h->SetLabelSize(0.03,
"Y");
508 h->GetYaxis()->ChangeLabel(-1, 0, 0);
511 h->SetLabelSize(0,
"Y");
515 h->SetLabelSize(0.03,
"X");
516 h->GetXaxis()->ChangeLabel(-1, -1, 0);
519 h->SetLabelSize(0,
"X");
521 h->GetXaxis()->SetRangeUser(0, 0.48);
522 h->GetXaxis()->SetDecimals();
524 else if (
auto st = dynamic_cast<THStack*>(obj) )
526 for (
auto o : *st->GetHists())
528 auto h =
dynamic_cast<TH1*
>(
o);
530 h->GetYaxis()->SetRangeUser(0, q0max);
537 auto column =
idx % 3;
538 double leftMargin = leftMarginMaster + (column * padWidth);
539 double rightMargin = rightMarginMaster + (3 - column - 1) * padWidth;
541 double topMargin = topMarginMaster + row * padHeight;
542 double bottomMargin = bottomMarginMaster + (nRows - row - 1) * padHeight;
545 TLatex *
l =
new TLatex(0.25, q0max*0.9, Form(
"%.1f #leq |#vec{q}|/(GeV/c) < %.1f", cIdx *
Q3_SLICE_SIZE, (cIdx+1)*Q3_SLICE_SIZE));
547 l->SetTextSize(0.025);
548 canvases[cIdx]->GetListOfPrimitives()->Add(l);
551 canvases[cIdx]->Update();
556 auto c = canvases[cIdx]->DrawClonePad();
558 p->SetMargin(leftMargin, rightMargin, bottomMargin, topMargin);
573 TLegend
leg(leftMarginMaster + 0.3*padWidth, 1 - (topMarginMaster+0.95*padHeight),
574 leftMarginMaster + padWidth, 1 - (topMarginMaster+0.2*padHeight));
576 auto pad =
dynamic_cast<TPad*
>(masterC.GetPrimitive(
"1"));
577 auto stack =
dynamic_cast<THStack*
>(
pad->GetPrimitive(
"stack"));
578 for (
auto o : *stack->GetHists())
580 auto h =
dynamic_cast<TH1*
>(
o);
581 auto toks = TString(
h->GetName()).Tokenize(
"{}");
582 leg.GetListOfPrimitives()->AddFirst(
new TLegendEntry(
h,
INT_LEGEND.at(toks->Last()->GetName()).c_str(),
"f"));
585 for (
auto o : *
pad->GetListOfPrimitives())
587 if ( (h = dynamic_cast<TH1*>(
o)) )
590 leg.GetListOfPrimitives()->AddFirst(
new TLegendEntry(h,
"ND Data",
"lpe"));
593 TLatex xaxisTitle(0.5, 0.03,
"Visible E_{had} (GeV)");
595 xaxisTitle.SetTextAlign(22);
596 xaxisTitle.SetTextSize(0.04);
598 TLatex yaxisTitle(0.025, 0.5, Form(
"10^{%d} Events",
int(-
log10(
Q0Q3_SCALE))));
599 yaxisTitle.SetTextAngle(90);
601 yaxisTitle.SetTextAlign(22);
602 yaxisTitle.SetTextSize(0.04);
605 TLatex beamText(leftMarginMaster + 0.02, 1 - topMarginMaster + 0.01,
606 (fHornCurrent ==
"FHC" ?
"Neutrino beam" :
"Antineutrino beam"));
608 beamText.SetTextAlign(11);
609 beamText.SetTextSize(0.03);
612 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
613 masterC.Print(Form(
"%s/%s/%s.%s", fOutDir.c_str(),
key.mecWgt().c_str(),
"q0q3panel",
ext));
626 if (fMCPlots.find(
"Q3Migration") != fMCPlots.end())
628 auto & spectrum = fMCPlots.at(
"Q3Migration").at(
PlotKey(
"NoMECWgt", intName));
631 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
632 c.Print( Form(
"%s/NoMECWgt/Q3Migration_%s.%s", fOutDir.c_str(), intName.c_str(),
ext) );
636 if (intName !=
"QE" && intName !=
"MEC")
638 if (fMCPlots.find(
"TrueQ0VsTrueQmag") == fMCPlots.end())
643 PlotKey pk(mecWgtName, intName);
645 auto & spectrum = fMCPlots.at(
"TrueQ0VsTrueQmag").at(
PlotKey(pk));
647 mcTH2->SetName( Form(
"TrueQ0VsTrueQmag{group=mecWgt,cat=%s}{group=interaction,cat=%s}", pk.
mecWgt().c_str(), pk.
interactionType().c_str()) );
648 this->WritePlotToHistFile(
"TrueQ0VsTrueQmag" + pk.
combination(), mcTH2,
false);
650 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
651 c.Print( Form(
"%s/%s/TrueQ0VsTrueQmag_%s.%s", fOutDir.c_str(), mecWgtName.c_str(), intName.c_str(),
ext) );
656 for (
const auto & plotListPair : fMCPlots)
658 const auto &
plotName = plotListPair.first;
659 const auto &
plots = plotListPair.second;
662 if (
plotName.find(
"Vs") != std::string::npos)
666 THStack shapeStack, sumStack;
667 TLegend
legend(0.7, 0.5, 0.9, 0.85);
668 legend.SetFillStyle(0);
669 TH1 * mcTH1 =
nullptr;
670 TH1 * mcTH1Area =
nullptr;
671 for (
const auto & intClass : INT_ORDER)
677 mcTH1Area =
dynamic_cast<TH1*
>(mcTH1->Clone());
678 mcTH1Area->Scale(1./mcTH1->Integral());
679 mcTH1->SetFillStyle(1001);
682 shapeStack.Add(mcTH1Area,
"hist");
683 sumStack.Add(mcTH1,
"hist");
685 legend.GetListOfPrimitives()->AddFirst(
new TLegendEntry(mcTH1, intClass.c_str(),
"l"));
689 sumStack.GetXaxis()->SetTitle(mcTH1->GetXaxis()->GetTitle());
690 sumStack.GetYaxis()->SetTitle(
"Events");
692 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
693 c.Print( Form(
"%s/%s/%s_MCSums.%s", fOutDir.c_str(), mecWgtName.c_str(),
plotName.c_str(),
ext) );
695 shapeStack.Draw(
"nostack");
696 shapeStack.GetXaxis()->SetTitle(mcTH1->GetXaxis()->GetTitle());
697 shapeStack.GetYaxis()->SetTitle(
"Fraction of events");
699 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
700 c.Print( Form(
"%s/%s/%s_MCShapes.%s", fOutDir.c_str(), mecWgtName.c_str(),
plotName.c_str(),
ext) );
710 const auto & wgtPlots = fMCPlots.at(
"weight");
715 TLegend
legend(0.7, 0.5, 0.9, 0.85);
716 legend.SetFillStyle(0);
723 if (wgtPlots.find(key) == wgtPlots.end())
729 const auto & spectrum = wgtPlots.at(key);
735 th1->Scale(1./integral);
737 th1->Draw((
"hist" + sameOpt).c_str());
741 this->WritePlotToHistFile(
"weight" + key.
combination(), th1,
false);
743 legend.GetListOfPrimitives()->AddFirst(
new TLegendEntry(th1, key.
interactionType().c_str(),
"l"));
748 for (
const auto &
ext : {
"png",
"root",
"eps",
"pdf"})
749 c.Print( Form(
"%s/weights/%s.%s", fOutDir.c_str(),
mecWgt.c_str(),
ext) );
764 TH1::AddDirectory(
false);
765 Plotter pl(hornCurrent, prodName,
"/nova/ana/users/jwolcott/2p2h_tuning/2018/plots/kinematics");
const std::vector< std::string > MECQE_WGTS
T max(const caf::Proxy< T > &a, T b)
void split(double tt, double *fr)
void WritePlotToHistFile(const std::string &plotName, TH1 *hist, bool isData=true)
void KinematicsPlots(std::string hornCurrent, std::string prodName)
std::unordered_set< std::string > fMCPlotsWritten
std::map< std::string, std::unique_ptr< ana::Spectrum > > fDataPlots
std::unordered_set< std::string > fdataPlotsWritten
bool operator<(const PlotKey &other) const
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
bool SplitMCPlotName(const std::string &keyedName, std::string &name, PlotKey &plotkey)
PlotKey(const std::string &wgtnm, const std::string &intnm)
const std::unordered_map< std::string, double > Q0Q3_MAX
std::ostream & operator<<(std::ostream &os, PlotKey const &pk)
std::string mecWgt() const
std::string interactionType() const
const double MC_NOMINAL_POT
const std::string & combination() const
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
TCanvas * canvas(const char *nm, const char *ti, const char *a)
TCanvas * DataMCComparison(const std::string &plotName, const PlotKey &key, ana::Spectrum *dataPlot, const std::map< PlotKey, const ana::Spectrum * > &mcplots, const PlotKey &plotKey)
const std::map< std::string, int > INT_COLORS
Representation of a spectrum in any variable, with associated POT.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
const std::vector< Plot > plots
std::map< std::string, std::map< PlotKey, std::unique_ptr< ana::Spectrum > > > fMCPlots
std::vector< double > POT
void out_of_range(const char *function, int max, int index, const char *msg1="", const char *msg2="")
const unsigned int Q0Q3_SLICE_OFFSET
std::map< PlotKey, std::vector< std::unique_ptr< TCanvas > > > fq0q3Canvases
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
std::string to_string(ModuleType mt)
Prevent histograms being added to the current directory.
Plotter(const std::string &hornCurrent, const std::string &prodName, const std::string &outDir)
const double Q3_SLICE_SIZE
const std::vector< std::string > INT_ORDER
const std::map< std::string, std::string > INT_LEGEND
std::string _interactionType