test_genieweights.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Binning.h"
2 #include "CAFAna/Cuts/Cuts.h"
3 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Systs/Systs.h"
10 
11 #include "TCanvas.h"
12 #include "TH1.h"
13 #include "TFile.h"
14 #include "TLegend.h"
15 #include "TLine.h"
16 #include "TLatex.h"
17 #include "THStack.h"
19 
20 using namespace ana;
21 
23 
24  TString filename = "test_genieweight_v0.root";
25  TString imgname ="plot_test_genieweight.png";
26  //Options from Vars/GenieVars.h:
27  //kTuftsWeightCC CCQE CCnonQE; kUMNWeightCC CCQE CCnonQE; kRPAWeightCCQE
28  const Var kGenieWeight = kTuftsWeightCC;
29 
30  SpectrumLoader loaderNDMC
31  ("/nova/prod/mc/S15-05-22/genie/nd/decaf/numu_contain/000108/**/neardet_genie_fhc_nonswap_nogenierw_2000_r000**_s*_c000*");
32 
33  loaderNDMC.SetSpillCut(kStandardSpillCuts);
34 
35  const Cut numuCC([](const caf::SRProxy* sr)
36  {
37  if(sr->mc.nnu == 0) return false;
38  assert(sr->mc.nnu == 1);
39  return (sr->mc.nu[0].iscc && abs(sr->mc.nu[0].pdg)==14);
40  });
41  //*******************
42  const HistAxis axisEhad(
43  "E_{had} (GeV)", Binning::Simple(50, 0, 1), kHadE);
44  //*******************
45  Spectrum sp_QE_noW (loaderNDMC, axisEhad, numuCC && kIsQE, kNoShift);
46  Spectrum sp_RES_noW (loaderNDMC, axisEhad, numuCC && kIsRes,kNoShift);
47  Spectrum sp_DIS_noW (loaderNDMC, axisEhad, numuCC && kIsDIS,kNoShift);
48  Spectrum sp_Coh_noW (loaderNDMC, axisEhad, numuCC && kIsCoh,kNoShift);
49 
50  Spectrum sp_QE_Wei (loaderNDMC, axisEhad, numuCC && kIsQE, kNoShift,kGenieWeight);
51  Spectrum sp_RES_Wei (loaderNDMC, axisEhad, numuCC && kIsRes,kNoShift,kGenieWeight);
52  Spectrum sp_DIS_Wei (loaderNDMC, axisEhad, numuCC && kIsDIS,kNoShift,kGenieWeight);
53  Spectrum sp_Coh_Wei (loaderNDMC, axisEhad, numuCC && kIsCoh,kNoShift,kGenieWeight);
54 
55  loaderNDMC.Go();
56 
57  TFile * file = new TFile (filename,"recreate");
58  sp_QE_noW. SaveTo(file, "sp_QE_noweight");
59  sp_RES_noW.SaveTo(file, "sp_RES_noweight");
60  sp_DIS_noW.SaveTo(file, "sp_DIS_noweight");
61  sp_Coh_noW.SaveTo(file, "sp_Coh_noweight");
62 
63  sp_QE_Wei. SaveTo(file, "sp_QE_weight");
64  sp_RES_Wei.SaveTo(file, "sp_RES_weight");
65  sp_DIS_Wei.SaveTo(file, "sp_DIS_weight");
66  sp_Coh_Wei.SaveTo(file, "sp_Coh_weight");
67 
68 
69  //return;
70 
71  //****************************************************
72  double thePOT = 1E20;
73 
74  auto hQEno = sp_QE_noW.ToTH1(thePOT);
75  auto hRESno = sp_RES_noW.ToTH1(thePOT);
76  auto hDISno = sp_DIS_noW.ToTH1(thePOT);
77  auto hCohno = sp_Coh_noW.ToTH1(thePOT);
78 
79  hQEno->SetLineColor(kMagenta);
80  hRESno->SetLineColor(kGreen+2);
81  hDISno->SetLineColor(kRed);
82  hCohno->SetLineColor(kBlue);
83 
84  hQEno->SetLineStyle(2);
85  hRESno->SetLineStyle(2);
86  hDISno->SetLineStyle(2);
87  hCohno->SetLineStyle(2);
88 
89  THStack * hsno = new THStack ("noW","Stacked histograms; E_{had} (GeV); Events/1E20 POT");
90  hsno->Add(hCohno,"histo");
91  hsno->Add(hQEno,"histo");
92  hsno->Add(hRESno,"histo");
93  hsno->Add(hDISno,"histo");
94 
95  auto hQE = sp_QE_Wei.ToTH1(thePOT);
96  auto hRES = sp_RES_Wei.ToTH1(thePOT);
97  auto hDIS = sp_DIS_Wei.ToTH1(thePOT);
98  auto hCoh = sp_Coh_Wei.ToTH1(thePOT);
99 
100  hQE->SetMarkerColor(kMagenta);
101  hRES->SetMarkerColor(kGreen+2);
102  hDIS->SetMarkerColor(kRed);
103  hCoh->SetMarkerColor(kBlue);
104 
105  hQE->SetMarkerStyle(29);
106  hRES->SetMarkerStyle(29);
107  hDIS->SetMarkerStyle(29);
108  hCoh->SetMarkerStyle(29);
109 
110  hQE->SetMarkerSize(2);
111  hRES->SetMarkerSize(2);
112  hDIS->SetMarkerSize(2);
113  hCoh->SetMarkerSize(2);
114 
115 
116  THStack * hsW = new THStack ("W","Stacked histograms; E_{had} (GeV); Events/1E20 POT");
117  hsW->Add(hCoh,"hist p");
118  hsW->Add(hQE,"hist p");
119  hsW->Add(hRES,"hist p");
120  hsW->Add(hDIS,"hist p");
121 
122  hsno->SetMaximum(max(hsno->GetMaximum(),hsW->GetMaximum()));
123  hsno->Draw();
124  hsW->Draw("same");
125 
126  TLegend* leg = new TLegend(0.6,0.6,0.89,0.89);
127  leg->AddEntry(hDIS,"DIS (weighed)", "p");
128  leg->AddEntry(hRES,"RES (weighed)", "p");
129  leg->AddEntry(hQE, " QE (weighed)", "p");
130  leg->AddEntry(hCoh,"Coh (weighed)", "p");
131  leg->Draw("same");
132  gPad->Print(imgname);
133 
134  return;
135 }
const Var kHadE
Definition: NumuVars.h:23
const Cut kIsQE
Definition: TruthCuts.cxx:104
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
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 Cut kIsRes
Definition: TruthCuts.cxx:111
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
void abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Cut kIsDIS
Definition: TruthCuts.cxx:118
string filename
Definition: shutoffs.py:106
void SetSpillCut(const SpillCut &cut)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
const SystShifts kNoShift
Definition: SystShifts.cxx:21
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TFile * file
Definition: cellShifts.C:17
assert(nhit_max >=nhit_nbins)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kGreen
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kIsCoh
Definition: TruthCuts.cxx:133
void test_genieweights()