KinematicsPlots.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void KinematicsPlots(std::string hornCurrent, std::string prodName)
3 {
4  std::cout << "CINT not supported (or tolerated) here" << std::endl;
5 }
6 #else
7 
8 #include <iostream>
9 #include <map>
10 #include <memory>
11 #include <string>
12 #include <typeinfo>
13 #include <unordered_map>
14 #include <unordered_set>
15 
16 #include "boost/algorithm/string.hpp"
17 #include "boost/tokenizer.hpp"
18 
19 #include "TCanvas.h"
20 #include "TError.h"
21 #include "TFile.h"
22 #include "TH1.h"
23 #include "TH2.h"
24 #include "THStack.h"
25 #include "TKey.h"
26 #include "TLegend.h"
27 #include "TLegendEntry.h"
28 #include "TList.h"
29 #include "TString.h"
30 #include "TLatex.h"
31 #include "TROOT.h"
32 
33 #include "CAFAna/Core/Spectrum.h"
34 #include "CAFAna/Core/Utilities.h"
35 
36 // /////////////////////////////////////
37 //
38 // custom types for use below
39 
40 const std::map<std::string, int> INT_COLORS = {
41  { "QE", kAzure+7 },
42  { "RES", kGreen+2 },
43  { "DIS", kGray },
44  { "MEC", kOrange-1 },
45  { "Other", kBlack }
46 };
47 
48 const std::vector<std::string> INT_ORDER = {
49  "Other",
50  "DIS",
51  "RES",
52  "QE",
53  "MEC"
54 };
55 
56 const std::map<std::string, std::string> INT_LEGEND = {
57  { "QE", "QE" },
58  { "RES", "RES" },
59  { "DIS", "DIS" },
60  { "MEC", "NOvA 2p2h" },
61  { "Other", "Other" }
62 };
63 
64 // would probably be better to determine all of these dynamically from the keys in the file...
65 // but I don't feel like it right now
66 const std::vector<std::string> MECQE_WGTS = {
67  "CV2018",
68  "CV2018-NoMEC",
69  "NOvA-MECdown",
70  "NOvA-MECup",
71  "MINERvA-1p1h",
72  "MINERvA-2p2h",
73  "MINERvA-2p2hnp",
74  "MINERvA-2p2hpp",
75 };
76 
77 // again, would be nice to determine dynamically, but this is easier
78 const std::unordered_map<std::string, double> Q0Q3_MAX = {
79  {"FHC", 8e4},
80 // {"FHC", 2e3},
81  {"RHC", 3e4},
82 };
83 
84 // scale down by this much. hard-code for now...
85 const double Q0Q3_SCALE = 1e-4;
86 
87 const double MC_NOMINAL_POT = 9e20; // for MC-only plots
88 
89 const double Q3_SLICE_SIZE = 0.1;
90 const unsigned int Q0Q3_SLICE_OFFSET = 1;
91 
92 class PlotKey
93 {
94  public:
95  PlotKey() = default;
96  PlotKey(const std::string & wgtnm, const std::string & intnm)
97  : _mecWgt(wgtnm), _interactionType(intnm),
98  _combination(Form("{%s}{%s}", wgtnm.c_str(), intnm.c_str()))
99  {};
100 
101  PlotKey(const PlotKey & other) = default;
102 
103  bool operator<(const PlotKey& other) const {
104  if (this->_mecWgt == other.mecWgt())
105  return this->_interactionType < other.interactionType();
106  else
107  return this->_mecWgt < other.mecWgt();
108  };
109 
110  std::string mecWgt() const { return this->_mecWgt; };
111  std::string interactionType() const { return this->_interactionType; };
112 
113  const std::string & combination() const { return this->_combination; };
114  private:
117 
118  // need this for ordering in the map... sigh.
120 };
121 
122 std::ostream & operator<<(std::ostream &os, PlotKey const &pk)
123 {
124  return os << pk.combination();
125 }
126 
127 // /////////////////////////////////////
128 //
129 // helper obj
130 
131 
132 class Plotter
133 {
134  public:
135  Plotter(const std::string & hornCurrent, const std::string & prodName, const std::string & outDir)
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" ))
140  {
141  this->GetDataPlots();
142  this->GetMCPlots();
143  }
144 
146  {
147  fHistFileData->Close();
148  fHistFileMC->Close();
149 
150  delete fHistFileData;
151  delete fHistFileMC;
152  }
153 
154  private:
158 
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;
161 
162  std::map<PlotKey, std::vector<std::unique_ptr<TCanvas>>> fq0q3Canvases; //ugh
163 
164  TFile * fHistFileData;
165  TFile * fHistFileMC;
166 
167  std::unordered_set<std::string> fdataPlotsWritten;
168  std::unordered_set<std::string> fMCPlotsWritten;
169 
170  void WritePlotToHistFile(const std::string & plotName, TH1 * hist, bool isData=true)
171  {
172  auto name = plotName.c_str();
173 
174  TFile * f = (isData) ? fHistFileData : fHistFileMC;
175  decltype(fdataPlotsWritten) & cache = isData ? fdataPlotsWritten : fMCPlotsWritten;
176 
177  if (cache.find(name) != cache.end())
178  return;
179 
180  f->cd();
181  hist->SetName(name);
182  hist->Write(name);
183  cache.insert(name);
184  }
185 
187  {
188  ana::DontAddDirectory inclGuard;
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();
191 
192  fDataPlots.clear();
193  TIter it(keys);
194  std::for_each(it.Begin(), TIter::End(), [this](TObject * listItem){
195  TKey * key = dynamic_cast<TKey*>(listItem);
196 
197  TDirectory * dir = dynamic_cast<TDirectory*>(key->ReadObj());
198  if (!dir)
199  return;
200 
201  fDataPlots.emplace(key->GetName(), std::move(ana::Spectrum::LoadFrom(dir)));
202  });
203  }
204 
205  ////////
206 
207  bool SplitMCPlotName(const std::string& keyedName, std::string & name, PlotKey & plotkey)
208  {
209  std::vector<std::string> fields;
210 
211  typedef boost::tokenizer<boost::char_separator<char> >
212  tokenizer;
213  boost::char_separator<char> sep("{}");
214  tokenizer tokens(keyedName, sep);
215 
216 // std::cout << "Parsing keyed name: " << keyedName << std::endl;
217  for (const auto & f : tokens)
218  {
219  // std::cout << " raw = " << f;
220  if (std::count(f.begin(), f.end(), '=') < 2) // if it's a named category, it has a 'group=' and a 'cat='
221  fields.push_back(f);
222  else
223  {
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());
227  }
228 
229  // std::cout << "; parsed = " << fields.back() << std::endl;
230  }
231 
232  // truth variables have only a var name and an interaction class...
233  if (fields.size() == 2)
234  plotkey = PlotKey("", fields[1]);
235  // ... while other plots have all the relevant fields: mec weight name, interaction class...
236  else if (fields.size() == 3)
237  plotkey = PlotKey(fields[1], fields[2]);
238  else
239  {
240  std::cerr << Form("Warning: Can't parse Spectrum name: '%s'", keyedName.c_str()) << std::endl;
241  return false;
242  }
243 
244  name = fields[0];
245 
246 // std::cout << " name, computed key = " << name << ", " << plotkey << std::endl;
247 
248  return true;
249  }
250 
251  ////////
252 
253  void GetMCPlots()
254  {
255  ana::DontAddDirectory inclGuard;
256  fMCPlots.clear();
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());
259  TFile f(fileName);
260  if (f.IsZombie())
261  return;
262 
263  TList * keys = f.GetListOfKeys();
264 
265  TIter it(keys);
266  std::for_each(it.Begin(), TIter::End(), [this](TObject * listItem){
267  TKey * key = dynamic_cast<TKey*>(listItem);
268 
269  TDirectory * dir = dynamic_cast<TDirectory*>(key->ReadObj());
270  if (!dir)
271  return;
272 
274  PlotKey plotkey;
275  bool success = SplitMCPlotName(key->GetName(), name, plotkey);
276  if (!success)
277  return;
278 
279  // if slot not already occupied OR preceding entry has nothing in it
280  if (fMCPlots[name].find(plotkey) == fMCPlots[name].end() || fMCPlots[name].at(plotkey)->ToTH1(fMCPlots[name].at(plotkey)->POT())->GetEntries() == 0)
281  fMCPlots[name][plotkey] = std::move(ana::Spectrum::LoadFrom(dir));
282  });
283 
284  std::cout << "MC plots: ";
285  for (const auto & plotPair : fMCPlots)
286  {
287  std::cout << plotPair.first << ":\n";
288  // for (const auto & keyPair : plotPair.second)
289  // std::cout << " " << keyPair.first << "\n";
290  std::cout << std::endl;
291  }
292  // throw std::runtime_error("look these over first");
293 
294  }
295 
296  ////////
297 
299  const PlotKey & key,
300  ana::Spectrum * dataPlot,
301  const std::map<PlotKey,
302  const ana::Spectrum*> & mcplots,
303  const PlotKey & plotKey
304  )
305  {
306  // these guys are unfortunately going to get leaked.
307  // it's too much work otherwise to make sure they persist
308  // through this function's return until the Print()ing
309  // of the canvas. this macro is small anyway
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);
314 
315  // this one, on the other hand, is returned and stuck
316  // into a unique_ptr, so he's ok
317  TCanvas * canvas = new TCanvas;
318  double dataPOT = dataPlot->POT();
319  static bool announcedPOT = false;
320  if (!announcedPOT)
321  {
322  std::cout << " data POT = " << dataPOT << std::endl;
323  announcedPOT = true;
324  }
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) // this is a ridiculous way to decide if this is one of the q0,q3 slices, but the histogram name isn't available here
329  dataTH1->SetMarkerSize(1.2);
330  else
331  dataTH1->SetMarkerSize(0.8);
332  dataTH1->Draw("pe"); // draw the data first so that the labeled axes show up. don't feel like copying them to the THStack
333  this->WritePlotToHistFile(plotName.c_str(), dataTH1, true);
334 
335  for (const auto & intClass : INT_ORDER)
336  {
337  PlotKey key(plotKey.mecWgt(), intClass);
338  TH1 * mcTH1 = mcplots.at(key)->ToTH1(dataPOT);
339  mcTH1->SetLineWidth(1);
340 // std::cout << " integral = " << mcTH1->Integral() << std::endl;
341  mcTH1->SetFillColor(INT_COLORS.at(intClass));
342  this->WritePlotToHistFile(plotName + key.combination(), mcTH1, false);
343  stack->Add(mcTH1, "hist");
344  // do this the hard way to get the legend order the same as the stack instead of reversed
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"));
347  }
348  stack->Draw("same");
349 
350  TH1 * mcsum = dynamic_cast<TH1*>((*stack->GetStack())[stack->GetStack()->LastIndex()]);
351 
352  auto maxVal = std::max(dataTH1->GetMaximum(), mcsum->GetMaximum());
353 
354  maxVal *= 1.2;
355  dataTH1->SetMaximum(maxVal);
356  dataTH1->GetYaxis()->SetRange(0, maxVal);
357 
358  dataTH1->Draw("pe same"); // yep, do it again, to move the points back up to the top
359 
360  legend->Draw();
361 
362  return canvas;
363  }
364 
365  ////////
366 
367  public:
368 
370  {
371  for (const auto & dataPair : fDataPlots)
372  {
373  const std::string & plotName = dataPair.first;
374  const std::unique_ptr<ana::Spectrum> & dataPlot = dataPair.second;
375 
376  for (const auto & mecWgtName : MECQE_WGTS)
377  {
378  if (fMCPlots.find(plotName) == fMCPlots.end())
379  {
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;
381  continue;
382  }
383 
384  // std::cout << " MC plots " << (MCPlots.find(plotName) == MCPlots.end() ? "DO NOT have" : "HAVE") << " plot with key '" << plotName << "' in them" << std::endl;
385  PlotKey key(mecWgtName, "");
386  const auto & plotCollection = fMCPlots.at(plotName);
387 
388  // std::cout << " MC plots with name '" << plotName << "' are: " << std::endl;
389  // for (const auto & plotPair : plotCollection)
390  // std::cout << " " << plotPair.first << std::endl;;
391 
392  // find the set of keys with the same weights but possibly different interaction types
393  // typename std::remove_reference<decltype(cont)>::type
394  std::map<PlotKey, const ana::Spectrum*> plotsByInteraction;
395  for (const auto & plotPair : plotCollection)
396  {
397  // std::cout << " candidate: " << plotPair.first;
398  if ( plotPair.first.mecWgt() != key.mecWgt())
399  {
400  // std::cout << " ... NOPE" << std::endl;
401  continue;
402  }
403  // std::cout << " YES" << std::endl;
404  plotsByInteraction.insert( std::make_pair(plotPair.first, plotPair.second.get()) );
405 // std::cout << plotPair.first << " " << plotPair.second->Integral(plotPair.second->POT()) << std::endl;;
406  }
407 
408  try
409  {
410  std::unique_ptr<TCanvas> dataMCComparison(DataMCComparison(plotName, key, dataPlot.get(), plotsByInteraction, key));
411 
412  std::size_t q3Loc = plotName.find("q3=");
413  if (q3Loc != std::string::npos)
414  {
415  TIter iter(dataMCComparison->GetListOfPrimitives());
416  TObject * obj;
417  while ( (obj = iter()) )
418  {
419  TH1 * h = dynamic_cast<TH1*>(obj);
420  if (!h)
421  continue;
422  h->GetXaxis()->SetRangeUser(0, 1);
423  h->GetYaxis()->SetRangeUser(0, Q0Q3_MAX.at(this->fHornCurrent));
424  }
425  }
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) );
428 
429  // stash the q0q3 plots away for later
430  if (q3Loc != std::string::npos)
431  {
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) // gotta watch out for tricky integer overflows
436  canvasColl.resize(q3SliceIdx+1);
437 // std::cout << " for plot " << plotName << ", concluded q3 start = " << lb << ", slice number = " << q3SliceIdx << std::endl;
438  canvasColl[q3SliceIdx] = std::move(dataMCComparison);
439  }
440  }
441  catch (std::out_of_range & e)
442  {
443  std::cout << " couldn't find plot '" << plotName << "' with key '" << key << "'. Skipping data-MC comparison." << std::endl;
444  continue;
445  }
446  } // for (mecWgtName)
447  } // for (dataPair)
448  } // DataMCComparisons()
449 
450  ////////
451 
452  void q0q3Panel()
453  {
454  if (this->fq0q3Canvases.empty())
455  return;
456 
457  double q0max = Q0Q3_MAX.at(this->fHornCurrent) * Q0Q3_SCALE;
458  q0max *= 0.98; // for axis autolabeling purposes
459 
460  for (const auto & keyPair : this->fq0q3Canvases)
461  {
462  auto & key = keyPair.first;
463  auto & canvases = keyPair.second;
464  std::cout << " making q0q3 panel for key: " << key << std::endl;
465 
466  const double bottomMarginMaster = 0.12;
467  const double leftMarginMaster = 0.1;
468  const double rightMarginMaster = 0.05;
469  const double topMarginMaster = 0.05;
470 
471  TCanvas masterC("q0q3panel", "(q0, |q|) panel", 700, 700);
472  auto nCanvases = canvases.size();
473 // masterC.Divide(3, nCanvases / 3, 0, 0);
474  auto nRows = nCanvases / 3;
475  double padWidth = (1. - leftMarginMaster - rightMarginMaster) / 3;
476  double padHeight = (1. - topMarginMaster - bottomMarginMaster) / nRows;
477 
478  // yes, skip 0. nothing interesting in that bin anywho
479  for (unsigned int cIdx = Q0Q3_SLICE_OFFSET; cIdx < nCanvases; cIdx++)
480  {
481  if (!canvases[cIdx])
482  {
483  std::cout << " null" << std::endl;
484  continue;
485  }
486 
487  TIter it_prim(canvases[cIdx]->GetListOfPrimitives());
488  bool adjustedData = false; // the data histogram sometimes is in there twice
489  while (TObject * obj = it_prim())
490  {
491  if (auto l = dynamic_cast<TLegend*>(obj))
492  canvases[cIdx]->GetListOfPrimitives()->Remove(l);
493  if ( auto h = dynamic_cast<TH1*>(obj) )
494  {
495 
496  h->SetTitle(";;");
497  if (!adjustedData)
498  {
499  adjustedData = true;
500  h->Scale(Q0Q3_SCALE);
501  }
502  h->GetYaxis()->SetRangeUser(0, q0max);
503  h->GetYaxis()->SetDecimals();
504 
505  if ((cIdx - Q0Q3_SLICE_OFFSET) % 3 == 0)
506  {
507  h->SetLabelSize(0.03, "Y");
508  h->GetYaxis()->ChangeLabel(-1, 0, 0); // erase the last axis label so it doesn't clash with the 0 above it
509  }
510  else
511  h->SetLabelSize(0, "Y");
512 
513  if ((cIdx-Q0Q3_SLICE_OFFSET) / 3 >= (nCanvases-Q0Q3_SLICE_OFFSET) / 3 - 1)
514  {
515  h->SetLabelSize(0.03, "X");
516  h->GetXaxis()->ChangeLabel(-1, -1, 0); // erase the last axis label so it doesn't clash with the 0 next to it
517  }
518  else
519  h->SetLabelSize(0, "X");
520 
521  h->GetXaxis()->SetRangeUser(0, 0.48); // not 0.5 so the last label isn't shown
522  h->GetXaxis()->SetDecimals();
523  }
524  else if ( auto st = dynamic_cast<THStack*>(obj) )
525  {
526  for (auto o : *st->GetHists())
527  {
528  auto h = dynamic_cast<TH1*>(o);
529  h->Scale(Q0Q3_SCALE);
530  h->GetYaxis()->SetRangeUser(0, q0max);
531  }
532  st->Modified();
533  }
534  }
535 
536  auto idx = cIdx - Q0Q3_SLICE_OFFSET;
537  auto column = idx % 3;
538  double leftMargin = leftMarginMaster + (column * padWidth);
539  double rightMargin = rightMarginMaster + (3 - column - 1) * padWidth;
540  auto row = idx / 3;
541  double topMargin = topMarginMaster + row * padHeight;
542  double bottomMargin = bottomMarginMaster + (nRows - row - 1) * padHeight;
543 
544  masterC.cd();
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));
546  l->SetTextAlign(22); // middle-middle
547  l->SetTextSize(0.025);
548  canvases[cIdx]->GetListOfPrimitives()->Add(l);
549 // canvases[cIdx]->GetListOfPrimitives()->Print();
550 
551  canvases[cIdx]->Update();
552 
553  auto p = new TPad(std::to_string(cIdx).c_str(), "", 0, 0, 1, 1);
554  p->cd();
555 // gROOT->SetSelectedPad(canvases[cIdx].get());
556  auto c = canvases[cIdx]->DrawClonePad();
557 // gROOT->GetListOfCanvases()->Remove(c);
558  p->SetMargin(leftMargin, rightMargin, bottomMargin, topMargin);
559 // std::cout << " left right top bottom margins are " << leftMargin << " " << rightMargin << " " << topMargin << " " << bottomMargin << std::endl;
560  p->SetFillStyle(0);
561  p->RedrawAxis();
562 
563  masterC.cd();
564  p->Draw();
565 
566  // canvases[cIdx]->Print();
567 
568  }
569 
570  masterC.cd();
571 
572  // rebuild the legend now
573  TLegend leg(leftMarginMaster + 0.3*padWidth, 1 - (topMarginMaster+0.95*padHeight),
574  leftMarginMaster + padWidth, 1 - (topMarginMaster+0.2*padHeight));
575  leg.SetFillStyle(0);
576  auto pad = dynamic_cast<TPad*>(masterC.GetPrimitive("1"));
577  auto stack = dynamic_cast<THStack*>(pad->GetPrimitive("stack"));
578  for (auto o : *stack->GetHists())
579  {
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"));
583  }
584  TH1 * h;
585  for (auto o : *pad->GetListOfPrimitives())
586  {
587  if ( (h = dynamic_cast<TH1*>(o)) )
588  break;
589  }
590  leg.GetListOfPrimitives()->AddFirst(new TLegendEntry(h, "ND Data", "lpe"));
591  leg.Draw();
592 
593  TLatex xaxisTitle(0.5, 0.03, "Visible E_{had} (GeV)");
594  xaxisTitle.SetNDC();
595  xaxisTitle.SetTextAlign(22);
596  xaxisTitle.SetTextSize(0.04);
597  xaxisTitle.Draw();
598  TLatex yaxisTitle(0.025, 0.5, Form("10^{%d} Events", int(-log10(Q0Q3_SCALE))));
599  yaxisTitle.SetTextAngle(90);
600  yaxisTitle.SetNDC();
601  yaxisTitle.SetTextAlign(22);
602  yaxisTitle.SetTextSize(0.04);
603  yaxisTitle.Draw();
604 
605  TLatex beamText(leftMarginMaster + 0.02, 1 - topMarginMaster + 0.01,
606  (fHornCurrent == "FHC" ? "Neutrino beam" : "Antineutrino beam"));
607  beamText.SetNDC();
608  beamText.SetTextAlign(11); // align left-bottom
609  beamText.SetTextSize(0.03);
610  beamText.Draw();
611 
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));
614  }
615  }
616 
617  ////////
618 
619  void TruthPlots()
620  {
621  TCanvas c;
622  // migration matrices: don't stack these
623  for (const auto & intName : INT_ORDER)
624  {
625  c.Clear();
626  if (fMCPlots.find("Q3Migration") != fMCPlots.end())
627  {
628  auto & spectrum = fMCPlots.at("Q3Migration").at(PlotKey("NoMECWgt", intName));
629  TH2 * mcTH2 = spectrum->ToTH2(MC_NOMINAL_POT);
630  mcTH2->Draw("colz");
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) );
633  }
634 
635  // these only exist for QE & Dytman-MEC
636  if (intName != "QE" && intName != "MEC")
637  continue;
638  if (fMCPlots.find("TrueQ0VsTrueQmag") == fMCPlots.end())
639  continue;
640  for (const auto & mecWgtName : MECQE_WGTS)
641  {
642  c.Clear();
643  PlotKey pk(mecWgtName, intName);
644 // std::cout << " trying plot: TrueQ0VsTrueQmag/" << pk << std::endl;
645  auto & spectrum = fMCPlots.at("TrueQ0VsTrueQmag").at(PlotKey(pk));
646  TH2 * mcTH2 = spectrum->ToTH2(MC_NOMINAL_POT);
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);
649  mcTH2->Draw("colz");
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) );
652 
653  }
654  }
655 
656  for (const auto & plotListPair : fMCPlots)
657  {
658  const auto & plotName = plotListPair.first;
659  const auto & plots = plotListPair.second;
660 // if (plotName.find("True") == std::string::npos)
661 // continue;
662  if (plotName.find("Vs") != std::string::npos)
663  continue;
664  for (const auto & mecWgtName : MECQE_WGTS)
665  {
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)
672  {
673  PlotKey key(mecWgtName, intClass);
674 // std::cout << " trying plot: " << plotName << "/" << key << " (" << ((plots.find(key) == plots.end()) ? "NOT FOUND" : "found") << ")" << std::endl;
675  mcTH1 = plots.at(key)->ToTH1(MC_NOMINAL_POT);
676  mcTH1->SetLineColor(INT_COLORS.at(intClass));
677  mcTH1Area = dynamic_cast<TH1*>(mcTH1->Clone());
678  mcTH1Area->Scale(1./mcTH1->Integral());
679  mcTH1->SetFillStyle(1001);
680  mcTH1->SetFillColor(INT_COLORS.at(intClass));
681  this->WritePlotToHistFile(plotName + key.combination(), mcTH1, false);
682  shapeStack.Add(mcTH1Area, "hist");
683  sumStack.Add(mcTH1, "hist");
684  // do this the hard way to get the legend order the same as the stack instead of reversed
685  legend.GetListOfPrimitives()->AddFirst(new TLegendEntry(mcTH1, intClass.c_str(), "l"));
686  }
687  c.Clear();
688  sumStack.Draw();
689  sumStack.GetXaxis()->SetTitle(mcTH1->GetXaxis()->GetTitle());
690  sumStack.GetYaxis()->SetTitle("Events");
691  legend.Draw();
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) );
694 
695  shapeStack.Draw("nostack");
696  shapeStack.GetXaxis()->SetTitle(mcTH1->GetXaxis()->GetTitle());
697  shapeStack.GetYaxis()->SetTitle("Fraction of events");
698  legend.Draw();
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) );
701  } // for (mecWgtName)
702  } // for (plotListPair)
703 
704  }
705 
706  ////////
707 
708  void WeightPlots()
709  {
710  const auto & wgtPlots = fMCPlots.at("weight");
711  for (const auto & mecWgt : MECQE_WGTS )
712  {
713  TCanvas c;
714  THStack stack;
715  TLegend legend(0.7, 0.5, 0.9, 0.85);
716  legend.SetFillStyle(0);
717 
718  std::string sameOpt = "";
719  for (const auto & intType : INT_ORDER)
720  {
721  PlotKey key (mecWgt, intType);
722  // std::cout << " looking for XS weight plot: " << xsWgt << " / " << intType << std::endl;
723  if (wgtPlots.find(key) == wgtPlots.end())
724  {
725  // std::cout << " ... but couldn't find it. skip." << std::endl;
726  continue;
727  }
728 
729  const auto & spectrum = wgtPlots.at(key); // too much typing to spell out the type (ha)
730 
731  TH1 * th1 = spectrum->ToTH1(MC_NOMINAL_POT, INT_COLORS.at(intType)); // unit normalized
732  // std::cout << th1->GetEntries() << std::endl;
733  double integral = th1->Integral();
734  if (integral > 0)
735  th1->Scale(1./integral);
736 
737  th1->Draw(("hist" + sameOpt).c_str());
738  if (sameOpt == "")
739  sameOpt = " same";
740 
741  this->WritePlotToHistFile("weight" + key.combination(), th1, false);
742 
743  legend.GetListOfPrimitives()->AddFirst(new TLegendEntry(th1, key.interactionType().c_str(), "l"));
744  }
745 
746  legend.Draw();
747 
748  for (const auto & ext : {"png", "root", "eps", "pdf"})
749  c.Print( Form("%s/weights/%s.%s", fOutDir.c_str(), mecWgt.c_str(), ext) );
750  } // for(xsWgt)
751  } // WeightPlots()
752 };
753 
754 // /////////////////////////////////////
755 //
756 // main
757 
758 
759 void KinematicsPlots(std::string hornCurrent, std::string prodName)
760 {
761  // too many "<image> file has been created" messages
762  gErrorIgnoreLevel = kWarning;
763 
764  TH1::AddDirectory(false);
765  Plotter pl(hornCurrent, prodName, "/nova/ana/users/jwolcott/2p2h_tuning/2018/plots/kinematics");
766  pl.DataMCComparisons();
767  pl.q0q3Panel();
768  pl.WeightPlots();
769  pl.TruthPlots();
770 }
771 
772 #endif
const std::vector< std::string > MECQE_WGTS
T max(const caf::Proxy< T > &a, T b)
void split(double tt, double *fr)
const XML_Char * name
Definition: expat.h:151
void WritePlotToHistFile(const std::string &plotName, TH1 *hist, bool isData=true)
void KinematicsPlots(std::string hornCurrent, std::string prodName)
keys
Reco plots.
Definition: caf_analysis.py:46
fileName
Definition: plotROC.py:78
std::unordered_set< std::string > fMCPlotsWritten
std::map< std::string, std::unique_ptr< ana::Spectrum > > fDataPlots
std::unordered_set< std::string > fdataPlotsWritten
set< int >::iterator it
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.
Definition: Spectrum.cxx:209
bool SplitMCPlotName(const std::string &keyedName, std::string &name, PlotKey &plotkey)
PlotKey()=default
PlotKey(const std::string &wgtnm, const std::string &intnm)
const char * p
Definition: xmltok.h:285
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
OStream cerr
Definition: OStream.cxx:7
const std::string & combination() const
void WeightPlots()
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
void GetDataPlots()
void DataMCComparisons()
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
std::string outDir
TFile * fHistFileData
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
void TruthPlots()
std::string fProdName
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
std::string _combination
std::string fHornCurrent
const std::vector< Plot > plots
def success(message)
Definition: log.py:5
void GetMCPlots()
double POT() const
Definition: Spectrum.h:231
std::map< std::string, std::map< PlotKey, std::unique_ptr< ana::Spectrum > > > fMCPlots
const char sep
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void out_of_range(const char *function, int max, int index, const char *msg1="", const char *msg2="")
const double Q0Q3_SCALE
void q0q3Panel()
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const unsigned int Q0Q3_SLICE_OFFSET
std::map< PlotKey, std::vector< std::unique_ptr< TCanvas > > > fq0q3Canvases
T log10(T number)
Definition: d0nt_math.hpp:120
TDirectory * dir
Definition: macro.C:5
std::string _mecWgt
void End(void)
Definition: gXSecComp.cxx:210
std::string fOutDir
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
Float_t e
Definition: plot.C:35
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
TFile * fHistFileMC
std::string _interactionType