ND_DataMC.C
Go to the documentation of this file.
1 // ND data/MC comparison plots for 3 flavor group.
2 // Makes POT norm with full syst band and area norm with shape only syst band
3 
7 #include "CAFAna/Cuts/Cuts.h"
13 #include "CAFAna/Cuts/SpillCuts.h"
14 #include "CAFAna/Core/Var.h"
17 #include "CAFAna/Systs/Systs.h"
22 #include "CAFAna/Vars/HistAxes.h"
25 #include "CAFAna/Vars/XsecTunes.h"
26 
27 #include "OscLib/IOscCalc.h"
28 
29 
30 #include "TFile.h"
31 #include "TSystem.h"
32 
33 using namespace ana;
34 
35 void ND_DataMC(TString sFHC_RHC, TString sNumu_Nue){
36  bool isRHC = sFHC_RHC.Contains("rhc") ? true : false;
37 
38  // --- datasets ---
39  TString fnameMC = "prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_"+sNumu_Nue+"2020";
40  TString fnameData = "prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_"+sNumu_Nue+"2020";
41 
42  if(isRHC){
43  fnameMC = "prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_"+sNumu_Nue+"2020";
44  fnameData = "prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_"+sNumu_Nue+"2020";
45  }
46 
47  SpectrumLoader ldrMC((std::string)fnameMC);
49 
50  SpectrumLoader ldrData((std::string)fnameData);
52 
53  // --- weights, shifts ---
55 
56  // --- output file ---
57  TString fileOut = sFHC_RHC+"_"+sNumu_Nue+".root";
58 
59  // --- numu-specific cuts, plots ---
60  if(sNumu_Nue.Contains("numu")){
61 
62  // selection cutflow
63  std::vector<Cut> selCuts = {kNumuQuality&&kNumuContainND2020}; //, kNumu2020ND};
64  std::vector<std::string> selNames = {"QualCont"}; //, "Full"};
65 
66 // std::vector<Cut> selCuts = {kNumu2020ND};
67 // std::vector<std::string> selNames = {"Full"};
68  unsigned int nSels = selNames.size();
69 
70  // quartiles
71  std::vector<Cut> hadEfracCuts = GetNumuEhadFracQuantCuts2020(isRHC, 4);
72  std::vector<std::string> hadEfracNames = {"hadEfracQ1", "hadEfracQ2", "hadEfracQ3", "hadEfracQ4"};
73  unsigned int nhadEfracQuants = hadEfracNames.size();
74 
75  std::vector<Cut> quantCuts;
76  std::vector<std::string> quantNames;
77  for(unsigned int q1 = 0; q1<nhadEfracQuants; ++q1){
78  std::vector<Cut> ptCuts = GetNumuPtQuantCuts2020(isRHC, q1+1, 3);
79  ptCuts.push_back(kNoCut);
80 
81  std::vector<std::string> ptNames = {"pTfracQ1", "pTfracQ2", "pTfracQ3", "pTfracAllQs"};
82  unsigned int nptQuants = ptNames.size();
83 
84  for(unsigned int q2 = 0; q2<nptQuants; ++q2){
85  quantCuts.push_back(hadEfracCuts[q1]&&ptCuts[q2]);
86  quantNames.push_back(hadEfracNames[q1]+"_"+ptNames[q2]);
87  }
88  }
89 
90  quantCuts.push_back(kNoCut);
91  quantNames.push_back("hadEfracAllQs_pTfracAllQs");
92  unsigned int nQuants = quantNames.size();
93 
94  // variables
95  std::vector<HistDef> defs;
97  unsigned int nVars = defs.size();
98  std::vector<std::string> varNames; std::vector<HistAxis> axis;
99  for(unsigned int var = 0; var<nVars; ++var){
100  varNames.emplace_back(defs[var].name);
101  axis.emplace_back(defs[var].axis);
102  }
103 
104  // systematics
105  std::vector<const ISyst*> systs = get3FlavorAna2020AllSysts(EAnaType2020::kNumuAna, true, BeamType2020::kFHC);
107  DataMCPair* numuDataMC[nVars][nSels][nQuants];
108 
109  Cut WrongSign = kIsAntiNu;
110  if(isRHC) WrongSign = kIsNu;
111 
112  for(unsigned int var = 0; var<nVars; ++var){
113  for(unsigned int sel = 0; sel<nSels; ++sel){
114  for(unsigned int quant = 0; quant<nQuants; ++quant){
115  Cut thisCut = selCuts[sel] && quantCuts[quant];
116  numuDataMC[var][sel][quant] = new DataMCPair(thisCut, axis[var], ldrData, ldrMC, weight, systs, {{"Total Background", !kIsNumuCC},{"Wrong Sign", WrongSign}});
117  }
118  }
119  }
120 
121  ldrMC.Go();
122  ldrData.Go();
123 
124  TFile* fout = new TFile(fileOut, "RECREATE");
125 
126  for(unsigned int var = 0; var<nVars; ++var){
127  for(unsigned int sel = 0; sel<nSels; ++sel){
128  for(unsigned int quant = 0; quant<nQuants; ++quant){
129  TString name = varNames[var]+"_"+selNames[sel]+"_"+quantNames[quant];
130  numuDataMC[var][sel][quant]->SaveTo(fout, name);
131  }
132  }
133  }
134  }
135 
136  // --- nue-specific cuts, plots ---
137  else{
138  // selection cutflow
139  double highedge;
140  if(sFHC_RHC.Contains("fhc")) highedge = kNue2020CVNFHCHighEdge;
141  else highedge = kNue2020CVNRHCHighEdge;
142  Cut CVNeHigh = kCVNe_looseptp > highedge;
143  std::vector<std::string> selNames = {"Presel", "CVNeLow", "CVNeHigh"};
144  std::vector<Cut> selCuts = {kNue2020NDPresel, kNue2020ND, kNue2020ND&&CVNeHigh};
145  unsigned int nSels = selNames.size();
146 
147  std::string decompFName = "nue2020_rhc_prop_decomp_ratios.root";
148  std::string decompShortName = ""; //"prop";
149  if(sFHC_RHC.Contains("fhc")) decompFName = "nue2020_fhc_combo_decomp_ratios.root";
150 
151  const Var kDecompWeight(DecompWeightFunc(FindCAFAnaDir()+
152  "/nue/decomp/" + decompFName,
153  decompShortName,
154  kNue2020AnaBin));
155  // variables
156  std::vector<HistDef> defs;
157  AddHistDefRecoND(defs);
158  AddHistDefSlice(defs);
159  AddHistDefNueND(defs);
160  AddHistDefNueShower(defs);
161  AddHistDefNuePID(defs);
163  AddHistDefProngCVN(defs);
164  AddHistDefNueEnergy(defs);
165  unsigned int nVars = defs.size();
166  std::vector<std::string> varNames; std::vector<HistAxis> axis;
167  for(unsigned int var = 0; var<nVars; ++var){
168  varNames.emplace_back(defs[var].name);
169  axis.emplace_back(defs[var].axis);
170  }
171 
172  Spectrum* specData[nVars][nSels];
173  IPrediction* predMC[nVars][nSels];
174  IPrediction* predMCdecomp[nVars][nSels];
175 
176  for(unsigned int var = 0; var<nVars; ++var){
177  for(unsigned int sel = 0; sel<nSels; ++sel){
178  Cut thisCut = selCuts[sel];
179  predMC[var][sel] = new PredictionNoExtrap(ldrMC, kNullLoader, kNullLoader, axis[var], thisCut, kNoShift, weight);
180  predMCdecomp[var][sel] = new PredictionNoExtrap(ldrMC, kNullLoader, kNullLoader, axis[var], thisCut, kNoShift, weight*kDecompWeight);
181  specData[var][sel] = new Spectrum(ldrData, axis[var], thisCut, kNoShift);
182  }
183  }
184 
185  ldrMC.Go();
186  ldrData.Go();
187 
188  TString fileOut1(fileOut); fileOut1.ReplaceAll(".root", "_mc.root");
189  TString fileOut2(fileOut); fileOut2.ReplaceAll(".root", "_data.root");
190  TString fileOut3(fileOut); fileOut3.ReplaceAll(".root", "_mc_decomp.root");
191 
192  TFile* fout_mc = new TFile(fileOut1, "RECREATE");
193  TFile* fout_data = new TFile(fileOut2, "RECREATE");
194  TFile* fout_mcdecomp = new TFile(fileOut3, "RECREATE");
195 
196  std::vector<TString> compNames = {"nue", "numu", "antinue", "antinumu", "nc"};
197  unsigned int nComp = compNames.size();
198 
199  for(unsigned int var = 0; var<nVars; ++var){
200  for(unsigned int sel = 0; sel<nSels; ++sel){
201  TString name = varNames[var]+"_"+selNames[sel];
202  predMC[var][sel]->ComponentCC(12,12).Unoscillated().SaveTo(fout_mc, name+"_nue");
203  predMC[var][sel]->ComponentCC(14,14).Unoscillated().SaveTo(fout_mc, name+"_numu");
204  predMC[var][sel]->ComponentCC(-12,-12).Unoscillated().SaveTo(fout_mc, name+"_antinue");
205  predMC[var][sel]->ComponentCC(-14,-14).Unoscillated().SaveTo(fout_mc, name+"_antinumu");
206  predMC[var][sel]->ComponentNCTotal().SaveTo(fout_mc, name+"_nc");
207 
208  predMCdecomp[var][sel]->ComponentCC(12,12).Unoscillated().SaveTo(fout_mcdecomp, name+"_nue");
209  predMCdecomp[var][sel]->ComponentCC(14,14).Unoscillated().SaveTo(fout_mcdecomp, name+"_numu");
210  predMCdecomp[var][sel]->ComponentCC(-12,-12).Unoscillated().SaveTo(fout_mcdecomp, name+"_antinue");
211  predMCdecomp[var][sel]->ComponentCC(-14,-14).Unoscillated().SaveTo(fout_mcdecomp, name+"_antinumu");
212  predMCdecomp[var][sel]->ComponentNCTotal().SaveTo(fout_mcdecomp, name+"_nc");
213 
214  specData[var][sel]->SaveTo(fout_data, name);
215  }
216  }
217  }
218 }
void AddHistDefNumuNDDataMC(std::vector< HistDef > &hd)
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void AddHistDefNueShower(std::vector< HistDef > &hd)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
void AddHistDefSlice(std::vector< HistDef > &hd)
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: DataMCPair.cxx:282
const double kNue2020CVNRHCHighEdge
Definition: NueCuts2020.h:53
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
void SetSpillCut(const SpillCut &cut)
void AddHistDefNueSelectionExtras(std::vector< HistDef > &hd)
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
const Cut kNue2020ND
Definition: NueCuts2020.h:178
const Cut kIsNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg > 0;})
Definition: TruthCuts.h:54
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Double_t q2[12][num]
Definition: f2_nu.C:137
std::vector< Cut > GetNumuPtQuantCuts2020(const bool isRHC, const unsigned int ehad_quant, const unsigned int nquants)
const double kNue2020CVNFHCHighEdge
Definition: NueCuts2020.h:48
const Cut kNue2020NDPresel
Definition: NueCuts2020.h:169
void AddHistDefRecoND(std::vector< HistDef > &hd)
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
virtual OscillatableSpectrum ComponentCC(int from, int to) const
const HistDef defs[kNumVars]
Definition: vars.h:15
void AddHistDefProngCVN(std::vector< HistDef > &hd)
std::vector< float > Spectrum
Definition: Constants.h:759
void AddHistDefNuePID(std::vector< HistDef > &hd)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
void ND_DataMC(TString sFHC_RHC, TString sNumu_Nue)
Definition: ND_DataMC.C:35
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
virtual Spectrum ComponentNCTotal() const
void AddHistDefNueND(std::vector< HistDef > &hd)
const Var kNue2020AnaBin([](const caf::SRProxy *sr){int selBin=kNue2020SelectionBin(sr);float nuE=kNueEnergy2020(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2020, official Binning.
Definition: NueCuts2020.h:191
std::vector< const ISyst * > get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap)
TString sNumu_Nue
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
const Cut kNumuQuality
Definition: NumuCuts.h:18
const std::string selNames[kNumSels]
Definition: vars.h:46
void AddHistDefNueEnergy(std::vector< HistDef > &hd)
const Var kCVNe_looseptp
Definition: Vars.cxx:36
const Cut kNumuContainND2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const caf::SRVector3DProxy &start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const caf::SRVector3DProxy &stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 40.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2020.h:31
std::vector< Cut > GetNumuEhadFracQuantCuts2020(const bool isRHC, const unsigned int nquants)
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105
enum BeamMode string