make_pi0_xcheck.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TH1.h"
3 #include "TF1.h"
4 #include "TFile.h"
5 #include "TLegend.h"
6 #include "TLegendEntry.h"
7 #include "TArrayD.h"
8 #include "TAxis.h"
9 #include "TLatex.h"
10 
11 #include "CAFAna/Core/Cut.h"
12 #include "CAFAna/Cuts/NueCutsFirstAna.h"
13 #include "CAFAna/Analysis/Plots.h"
14 #include "CAFAna/Analysis/Style.h"
16 #include "CAFAna/Core/Spectrum.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Vars/Vars.h"
21 
22 #include "CAFAna/nue/pi0xcheck/make_pi0_xcheck.h"
23 
25 
26 using namespace ana;
27 
28 void make_pi0_xcheck(bool isFHC = 1)
29 {
30  const std::string tag = isFHC ? "FHC" : "RHC";
31 
32  std::string fData = isFHC ? "prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns" :
33  "prod_caf_R17-11-14-prod4reco.neutron-respin.b_nd_numi_rhc_full_v1_goodruns" ;
34  std::string fMC = isFHC ? "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_fhc_nova_v08_full_v1" :
35  "prod_caf_R17-11-14-prod4reco.neutron-respin.c_nd_genie_nonswap_rhc_nova_v08_full_v1" ;
36 
38  SpectrumLoaderBase* loaderMC = new SpectrumLoader(fMC);
39 
40  std::cout << "MC: " << fMC << std::endl;
41  std::cout << "Data: " << fData << std::endl;
42 
44  loaderMC->SetSpillCut(kStandardDQCuts);
45 
47 
48  Spectrum* dataSpects[kNSels][kNPlots];
49  Spectrum* mcSpects[kNSels][kNPlots];
50  Spectrum* bkgSpects[kNSels][kNPlots];
51 
52  Spectrum* dataSpects2[kNSels][kN2D];
53  Spectrum* mcSpects2[kNSels][kN2D];
54  Spectrum* bkgSpects2[kNSels][kN2D];
55 
56  Spectrum* dataSpectsMulti[kNSels][kNMulti];
57  Spectrum* mcSpectsMulti[kNSels][kNMulti];
58  Spectrum* bkgSpectsMulti[kNSels][kNMulti];
59 
60  for(unsigned int i = 0; i < kNSels; ++i){
61  for(unsigned int j = 0; j < kNPlots; ++j){
62  dataSpects[i][j] = new Spectrum(plots[j].label, plots[j].bins, *loader, plots[j].var, sels[i].cut);
63  mcSpects[i][j] = new Spectrum(plots[j].label, plots[j].bins, *loaderMC, plots[j].var, sels[i].cut,
64  kNoShift, weight);
65  bkgSpects[i][j] = new Spectrum(plots[j].label, plots[j].bins, *loaderMC, plots[j].var, sels[i].cut && !kTruePi0,
66  kNoShift, weight);
67  }
68  }
69 
70  for(unsigned int i = 0; i < kNSels; ++i){
71  for(unsigned int j = 0; j < kNMulti; ++j){
72  dataSpectsMulti[i][j] = new Spectrum(plotsMulti[j].label, plotsMulti[j].bins, *loader, plotsMulti[j].var, sels[i].cut);
73  mcSpectsMulti[i][j] = new Spectrum(plotsMulti[j].label, plotsMulti[j].bins, *loaderMC, plotsMulti[j].var, sels[i].cut,
74  kNoShift, weight);
75  bkgSpectsMulti[i][j] = new Spectrum(plotsMulti[j].label, plotsMulti[j].bins, *loaderMC, plotsMulti[j].var, sels[i].cut && !kTruePi0,
76  kNoShift, weight);
77  }
78  }
79 
80  for(unsigned int i = 0; i < kNSels; ++i){
81  for(unsigned int j = 0; j < kN2D; ++j){
82  dataSpects2[i][j] = new Spectrum(plots2d[j].label1, plots2d[j].label2, *loader, plots2d[j].bins1, plots2d[j].var1,
83  plots2d[j].bins2, plots2d[j].var2, sels[i].cut);
84  mcSpects2[i][j] = new Spectrum(plots2d[j].label1, plots2d[j].label2, *loaderMC, plots2d[j].bins1, plots2d[j].var1,
85  plots2d[j].bins2, plots2d[j].var2, sels[i].cut,
86  kNoShift, weight);
87  bkgSpects2[i][j] = new Spectrum(plots2d[j].label1, plots2d[j].label2, *loaderMC, plots2d[j].bins1, plots2d[j].var1,
88  plots2d[j].bins2, plots2d[j].var2, sels[i].cut && !kTruePi0,
89  kNoShift, weight);
90  }
91  }
92 
93  loader->Go();
94  loaderMC->Go();
95 
96  TFile fout(("make_pi0_xcheck_output_"+tag+".root").c_str(), "RECREATE");
97 
98  for(unsigned int i=0;i<kNSels;++i){
99  for(unsigned int j=0;j<kNPlots;++j){
100  std::string mysuffix=plots[j].suffix + "_" + sels[i].suffix;
101 
102  std::string d = tag+"_data_"+mysuffix;
103  std::string m = tag+"_mc_"+mysuffix;
104  std::string b = tag+"_bg_"+mysuffix;
105 
106  dataSpects[i][j]->SaveTo(&fout, d.c_str());
107  mcSpects[i][j]->SaveTo(&fout, m.c_str());
108  bkgSpects[i][j]->SaveTo(&fout, b.c_str());
109  }
110  }
111 
112  for(unsigned int i=0;i<kNSels;++i){
113  for(unsigned int j=0;j<kN2D;++j){
114  std::string mysuffix=plots2d[j].suffix + "_" + sels[i].suffix;
115 
116  std::string d = tag+"_data_"+mysuffix;
117  std::string m = tag+"_mc_"+mysuffix;
118  std::string b = tag+"_bg_"+mysuffix;
119 
120  dataSpects2[i][j]->SaveTo(&fout, d.c_str());
121  mcSpects2[i][j]->SaveTo(&fout, m.c_str());
122  bkgSpects2[i][j]->SaveTo(&fout, b.c_str());
123  }
124  }
125 
126  for(unsigned int i=0;i<kNSels;++i){
127  for(unsigned int j=0;j<kNMulti;++j){
128  std::string mysuffix=plotsMulti[j].suffix + "_" + sels[i].suffix;
129 
130  std::string d = tag+"_data_"+mysuffix;
131  std::string m = tag+"_mc_"+mysuffix;
132  std::string b = tag+"_bg_"+mysuffix;
133 
134  dataSpectsMulti[i][j]->SaveTo(&fout, d.c_str());
135  mcSpectsMulti[i][j]->SaveTo(&fout, m.c_str());
136  bkgSpectsMulti[i][j]->SaveTo(&fout, b.c_str());
137  }
138  }
139 }
std::vector< PlotDef2D > plots2d
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void make_pi0_xcheck(bool isFHC=1)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
const unsigned int kN2D
const unsigned int kNSels
void SetSpillCut(const SpillCut &cut)
const unsigned int kNPlots
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
virtual void Go()=0
Load all the registered spectra.
const Cut sels[kNumSels]
Definition: vars.h:44
const std::vector< Plot > plots
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
Float_t d
Definition: plot.C:236
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
const double j
Definition: BetheBloch.cxx:29
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:570
const Cut kTruePi0([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);if(sr->mc.nu[0].prim.size()==0) return false;float Pi0Mass=0.134977;float KE=-10.0;float TE=-99.0;int nprim=sr->mc.nu[0].prim.size();int npi0=0;for(int i=0;i< nprim;i++){if(sr->mc.nu[0].prim[i].pdg==111){TE=sr->mc.nu[0].prim[i].p.E;KE=pow(pow(TE, 2.0)-pow(Pi0Mass, 2.0), 0.5);if(KE >=0.0) npi0++;}}if(npi0){return true;}return false;})
static bool isFHC
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const std::vector< PlotDefMulti > plotsMulti
Base class for the various types of spectrum loader.
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const unsigned int kNMulti
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
const hit & b
Definition: hits.cxx:21
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49