plot_kinematics_cafana.C
Go to the documentation of this file.
1 #ifdef __CINT__
3 #else
4 
5 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Vars/Vars.h"
8 using namespace ana;
9 
11 
12 #include "TCanvas.h"
13 #include "TH1.h"
14 #include "TLegend.h"
15 
16 void Compare(const Spectrum* gen, const Spectrum* gib, const std::string& elem)
17 {
18  TH1* hgen = gen->ToTH1(1e20, kBlue);
19  TH1* hgib = gib->ToTH1(1e20, kRed);
20 
21  if(hgen->GetBinContent(hgen->GetMaximumBin()) > hgib->GetBinContent(hgib->GetMaximumBin())){
22  hgen->SetTitle(("#nu_{#mu} CC on "+elem).c_str());
23  hgen->GetYaxis()->SetTitle("Events / 10^{20} POT");
24  hgen->Draw("hist");
25  hgib->Draw("hist same");
26  }
27  else{
28  hgib->SetTitle(("#nu_{#mu} CC on "+elem).c_str());
29  hgib->GetYaxis()->SetTitle("Events / 10^{20} POT");
30  hgib->Draw("hist");
31  hgen->Draw("hist same");
32  }
33 }
34 
36 {
37  SpectrumLoader lgenie("/nfs/raid5/bckhouse/gibuu/reprocess/*.orig.caf.root");
38  SpectrumLoader lgibuu("/nfs/raid5/bckhouse/gibuu/reprocess/*.gibuu.caf.root");
39 
40  const Cut kPreselAny({"mc.nnu", "mc.nu.iscc", "mc.nu.pdg",
41  "mc.nu.vtx.x", "mc.nu.vtx.y", "mc.nu.vtx.z"},
42  [](const caf::StandardRecord* sr)
43  {
44  if(sr->mc.nnu == 0) return false;
45  if(!sr->mc.nu[0].iscc || !sr->mc.nu[0].pdg == 14) return false;
46  if(sr->mc.nu[0].vtx.x < -200 || sr->mc.nu[0].vtx.x > +200 ||
47  sr->mc.nu[0].vtx.y < -200 || sr->mc.nu[0].vtx.y > +200 ||
48  sr->mc.nu[0].vtx.z < -100 || sr->mc.nu[0].vtx.z > 1600) return false;
49  return true;
50  });
51 
52  const Var kGiBUUWgt({"mc.nnu", "mc.nu.mode", "mc.nu.genweight"},
53  [](const caf::StandardRecord* sr)
54  {
55  if(sr->mc.nnu == 0) return 0.f;
56  if(sr->mc.nu[0].mode >= 0) return 1.f; // GENIE
57  return sr->mc.nu[0].genweight;
58  });
59 
60  const Var kLogWgt({"mc.nnu", "mc.nu.mode", "mc.nu.genweight"},
61  [](const caf::StandardRecord* sr)
62  {
63  if(sr->mc.nnu == 0) return -100.f;
64  if(sr->mc.nu[0].mode >= 0) return -100.f; // GENIE
65  return log10(sr->mc.nu[0].genweight);
66  });
67 
68  const Var kTgtZ({"mc.nnu", "mc.nu.tgtZ"},
69  [](const caf::StandardRecord* sr)
70  {
71  if(sr->mc.nnu == 0) return -1;
72  return sr->mc.nu[0].tgtZ;
73  });
74 
75  const Var kTrueX({"mc.nnu", "mc.nu.x"},
76  [](const caf::StandardRecord* sr)
77  {
78  if(sr->mc.nnu == 0) return -1.f;
79  return sr->mc.nu[0].x;
80  });
81 
82  const Var kTrueY({"mc.nnu", "mc.nu.y"},
83  [](const caf::StandardRecord* sr)
84  {
85  if(sr->mc.nnu == 0) return -1.f;
86  return sr->mc.nu[0].y;
87  });
88 
89  const Var kTrueQ2({"mc.nnu", "mc.nu.q2"},
90  [](const caf::StandardRecord* sr)
91  {
92  if(sr->mc.nnu == 0) return -1.f;
93  return sr->mc.nu[0].q2;
94  });
95 
96  const Var kTrueW2({"mc.nnu", "mc.nu.W2"},
97  [](const caf::StandardRecord* sr)
98  {
99  if(sr->mc.nnu == 0) return -1.f;
100  return sr->mc.nu[0].W2;
101  });
102 
103  struct Spects
104  {
105  Spectrum *E, *x, *y, *q2, *W2;
106  };
107 
108  std::map<int, Spects> ssgen, ssgib;
109 
110  Spectrum sZgen("Target Z", Binning::Simple(30,0,30), lgenie, kTgtZ, kPreselAny);
111  Spectrum sZgib("Target Z", Binning::Simple(30,0,30), lgibuu, kTgtZ, kPreselAny, kNoShift, kGiBUUWgt);
112 
113  Spectrum sWtgib("log_{10}(weight)", Binning::Simple(60,-2,+4), lgibuu, kLogWgt, kPreselAny, kNoShift, kGiBUUWgt);
114 
115  std::map<int, std::string> names = {{1, "H"},
116  {6, "C"},
117  {8, "O"},
118  {17, "Cl"},
119  {22, "Ti"},
120  {26, "Fe"}};
121  for(int Z: {1, 6, 8, 17, 22, 26}){
122  const Cut kPresel = kPreselAny && (kTgtZ == Z);
123 
124  ssgen[Z].E = new Spectrum("True neutrino energy", Binning::Simple(40,0,10), lgenie, kTrueE, kPresel);
125  ssgib[Z].E = new Spectrum("True neutrino energy", Binning::Simple(40,0,10), lgibuu, kTrueE, kPresel, kNoShift, kGiBUUWgt);
126 
127  ssgen[Z].x = new Spectrum("x", Binning::Simple(40,0,2), lgenie, kTrueX, kPresel);
128  ssgib[Z].x = new Spectrum("x", Binning::Simple(40,0,2), lgibuu, kTrueX, kPresel, kNoShift, kGiBUUWgt);
129 
130  ssgen[Z].y = new Spectrum("y", Binning::Simple(50,0,1), lgenie, kTrueY, kPresel);
131  ssgib[Z].y = new Spectrum("y", Binning::Simple(50,0,1), lgibuu, kTrueY, kPresel, kNoShift, kGiBUUWgt);
132 
133  ssgen[Z].q2 = new Spectrum("Q^{2}", Binning::Simple(40,0,2), lgenie, kTrueQ2, kPresel);
134  ssgib[Z].q2 = new Spectrum("Q^{2}", Binning::Simple(40,0,2), lgibuu, kTrueQ2, kPresel, kNoShift, kGiBUUWgt);
135 
136  ssgen[Z].W2 = new Spectrum("W^{2}", Binning::Simple(40,0,2), lgenie, kTrueW2, kPresel);
137  ssgib[Z].W2 = new Spectrum("W^{2}", Binning::Simple(40,0,2), lgibuu, kTrueW2, kPresel, kNoShift, kGiBUUWgt);
138  } // end for Z
139 
140  lgenie.Go();
141  lgibuu.Go();
142 
143  sZgen.ToTH1(1e20, kBlue)->Draw("hist");
144  sZgib.ToTH1(1e20, kRed)->Draw("hist same");
145 
146  TLegend* leg = new TLegend(.725, .725, .875, .875);
147  leg->SetFillStyle(0);
148  leg->AddEntry(sZgen.ToTH1(1e20, kBlue), "GENIE", "l");
149  leg->AddEntry(sZgib.ToTH1(1e20, kRed), "GiBUU", "l");
150  leg->Draw();
151 
152  gPad->Print("tgtZ.pdf");
153 
154  new TCanvas;
155 
156  sWtgib.ToTH1(1e20)->Draw("hist");
157  gPad->Print("genweight.pdf");
158 
159  for(int Z: {1, 6, 8, 17, 22, 26}){
160  new TCanvas;
161  Compare(ssgen[Z].E, ssgib[Z].E, names[Z]);
162  leg->Draw();
163  gPad->Print(("Enu_"+names[Z]+".pdf").c_str());
164 
165  new TCanvas;
166  Compare(ssgen[Z].x, ssgib[Z].x, names[Z]);
167  leg->Draw();
168  gPad->Print(("x_"+names[Z]+".pdf").c_str());
169 
170  new TCanvas;
171  Compare(ssgen[Z].y, ssgib[Z].y, names[Z]);
172  leg->Draw();
173  gPad->Print(("y_"+names[Z]+".pdf").c_str());
174 
175  new TCanvas;
176  Compare(ssgen[Z].q2, ssgib[Z].q2, names[Z]);
177  leg->Draw();
178  gPad->Print(("q2_"+names[Z]+".pdf").c_str());
179 
180  new TCanvas;
181 
182  Compare(ssgen[Z].W2, ssgib[Z].W2, names[Z]);
183  leg->Draw();
184  gPad->Print(("W2_"+names[Z]+".pdf").c_str());
185  }
186 }
187 
188 #endif
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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:148
const Var kTrueQ2
Definition: TruthVars.h:27
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Float_t Z
Definition: plot.C:38
Double_t q2[12][num]
Definition: f2_nu.C:137
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
Float_t E
Definition: plot.C:20
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
void plot_kinematics_cafana()
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
std::vector< float > Spectrum
Definition: Constants.h:570
const SystShifts kNoShift
Definition: SystShifts.cxx:21
const Var kTrueY
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The StandardRecord is the primary top-level object in the Common Analysis File trees.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
T log10(T number)
Definition: d0nt_math.hpp:120
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
void Compare(const Spectrum *gen, const Spectrum *gib, const std::string &elem)