ppfx_make_systs.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/HistAxis.h"
7 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Core/Binning.h"
11 #include "CAFAna/Core/Cut.h"
12 #include "CAFAna/Core/SystShifts.h"
13 #include "CAFAna/Cuts/Cuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "CAFAna/Cuts/NueCutsSecondAna.h"
18 #include "CAFAna/Cuts/TimingCuts.h"
19 #include "CAFAna/Systs/Systs.h"
20 
21 #include "TStopwatch.h"
22 #include "CAFAna/Core/Utilities.h"
23 #include "Utilities/func/MathUtil.h"
24 
25 #include "TFile.h"
26 #include "TRandom.h"
27 
28 #include <fstream>
29 #include <iostream>
30 #include <cmath>
31 
32 
33 using namespace ana;
34  //----------------------------------------------------------------------
35 
37 {
38 
39  const Cut kIsNumu(
40  [](const caf::SRProxy* sr)
41  {
42  if(sr->mc.nnu == 0) return false;
43  assert(sr->mc.nnu == 1);
44  return (abs(sr->mc.nu[0].pdg) == 14);
45  });
46 
47  const Cut kIsNue(
48  [](const caf::SRProxy* sr)
49  {
50  if(sr->mc.nnu == 0) return false;
51  assert(sr->mc.nnu == 1);
52  return (abs(sr->mc.nu[0].pdg) == 12);
53  });
54 
55  const Cut kIsRock(
56  [](const caf::SRProxy* sr)
57  {
58  if (sr->mc.nnu == 0)
59  return false;
60  if(fabs(sr->mc.nu[0].vtx.X())>180||fabs(sr->mc.nu[0].vtx.Y())>180||sr->mc.nu[0].vtx.Z()<0||sr->mc.nu[0].vtx.Z()>1250)return true;
61  return false;
62  });
63 
64 
65  const Var kTrueE(
66  [](const caf::SRProxy* sr){
67  if(sr->mc.nnu!=1)return (float)-5.0;
68  return (float)sr->mc.nu[0].E;
69  });
70 
71 
72  std::vector<Var> kPPFXFluxUnivWgt;
73 
74  for(unsigned int UnivIdx = 0; UnivIdx < 100; UnivIdx++){
75 
76  const Var tempPPFXWgt(
77  [UnivIdx](const caf::SRProxy* sr){
78  if (!sr->hdr.ismc)return 1.f;
79  if (sr->mc.nnu != 1)return 1.f;
80  if(sr->mc.nu[0].rwgt.ppfx.vuniv[UnivIdx] <= 0)return 1.f;
81  return (float)sr->mc.nu[0].rwgt.ppfx.vuniv[UnivIdx];
82  });
83 
84  kPPFXFluxUnivWgt.push_back(tempPPFXWgt);
85  }
86 
87  const Binning bins = Binning::Simple(400, 0, 10);
88  const HistAxis TrueEAxis("True Energy (GeV)", bins, kTrueE);
89 
90 
91  struct VarDef
92  {
94  HistAxis axis;
95  };
96 
97  struct CompDef
98  {
100  Cut cut;
101  };
102 
103  std::vector<const ISyst*> systs = {&kBeamHornCurrent, &kBeamSpotSize, &kBeamPosX, &kBeamPosY, &kBeamH1PosX, &kBeamH1PosY,
105 
106  SystShifts shifts;
107 
108  std::vector <VarDef> vars = {
109  {"trueE", TrueEAxis},
110  };
111 
112  std::vector <CompDef> comps = {
113  {"numu", kIsNumu && !kIsAntiNu},
114  {"nue", kIsNue && !kIsAntiNu},
115  {"anti_numu", kIsNumu && kIsAntiNu},
116  {"anti_nue", kIsNue && kIsAntiNu},
117  };
118 
119  std::vector <std::vector<Spectrum*>> spects_RHC_nd;
120  std::vector <std::vector<std::vector<Spectrum*>>> spects_RHC_univ_nd;
121  std::vector <std::vector<Spectrum*>> spects_RHC_fd_ns;
122  std::vector <std::vector<std::vector<Spectrum*>>> spects_RHC_univ_fd_ns;
123  std::vector <std::vector<Spectrum*>> spects_FHC_nd;
124  std::vector <std::vector<std::vector<Spectrum*>>> spects_FHC_univ_nd;
125  std::vector <std::vector<Spectrum*>> spects_FHC_fd_ns;
126  std::vector <std::vector<std::vector<Spectrum*>>> spects_FHC_univ_fd_ns;
127  std::string FileName = "ppfx_flux_smooth.root";
128 
129  const std::string nearNonSFHC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
130  const std::string farNonSFHC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1";
131  const std::string nearNonSRHC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
132  const std::string farNonSRHC = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1";
133 
134  SpectrumLoader loaderNDMCFHC(nearNonSFHC);
135  SpectrumLoader loaderFDMCFHC(farNonSFHC);
136  SpectrumLoader loaderNDMCRHC(nearNonSRHC);
137  SpectrumLoader loaderFDMCRHC(farNonSRHC);
138 
139  loaderNDMCFHC.SetSpillCut(kStandardDQCuts);
140  loaderFDMCFHC.SetSpillCut(kStandardDQCuts);
141  loaderNDMCRHC.SetSpillCut(kStandardDQCuts);
142  loaderFDMCRHC.SetSpillCut(kStandardDQCuts);
143 
144  // Define MC weight
145  const Var kMCWeight = kPPFXFluxCVWgtSmooth*kXSecCVWgt2018;
146 
147  int N_PPFX = 100;
148 
149  for(int varIdx = 0; varIdx < (int) vars.size(); ++varIdx){
150  std::vector<Spectrum*> temp_spec_RHC_nd;
151  std::vector<std::vector<Spectrum*>> spec_RHC_univ_nd;
152  std::vector<Spectrum*> temp_spec_RHC_fd_ns;
153  std::vector<std::vector<Spectrum*>> spec_RHC_univ_fd_ns;
154  std::vector<Spectrum*> temp_spec_FHC_nd;
155  std::vector<std::vector<Spectrum*>> spec_FHC_univ_nd;
156  std::vector<Spectrum*> temp_spec_FHC_fd_ns;
157  std::vector<std::vector<Spectrum*>> spec_FHC_univ_fd_ns;
158  const HistAxis& axis = vars[varIdx].axis;
159 
160  for(int compIdx = 0; compIdx < (int) comps.size(); ++compIdx){
161  std::vector<Spectrum*> temp_spec_RHC_univ_nd;
162  std::vector<Spectrum*> temp_spec_RHC_univ_fd_ns;
163  std::vector<Spectrum*> temp_spec_FHC_univ_nd;
164  std::vector<Spectrum*> temp_spec_FHC_univ_fd_ns;
165  temp_spec_RHC_nd.emplace_back(new Spectrum(loaderNDMCRHC, axis, comps[compIdx].cut && !kIsRock, kNoShift, kMCWeight));
166  temp_spec_RHC_fd_ns.emplace_back(new Spectrum(loaderFDMCRHC, axis, comps[compIdx].cut, kNoShift, kMCWeight));
167  temp_spec_FHC_nd.emplace_back(new Spectrum(loaderNDMCFHC, axis, comps[compIdx].cut && !kIsRock, kNoShift, kMCWeight));
168  temp_spec_FHC_fd_ns.emplace_back(new Spectrum(loaderFDMCFHC, axis, comps[compIdx].cut, kNoShift, kMCWeight));
169 
170  for(int UnivIndx = 0; UnivIndx < N_PPFX; UnivIndx++){
171 
172  const Var kPPFXFluxWgt = GetPPFXFluxUnivWgtSmooth(UnivIndx);
173  const Var kMCUnivWeight = kPPFXFluxWgt*kXSecCVWgt2018;
174  temp_spec_RHC_univ_nd.emplace_back(new Spectrum(loaderNDMCRHC, axis, comps[compIdx].cut && !kIsRock, kNoShift, kMCUnivWeight));
175  temp_spec_RHC_univ_fd_ns.emplace_back(new Spectrum(loaderFDMCRHC, axis, comps[compIdx].cut, kNoShift, kMCUnivWeight));
176  temp_spec_FHC_univ_nd.emplace_back(new Spectrum(loaderNDMCFHC, axis, comps[compIdx].cut && !kIsRock, kNoShift, kMCUnivWeight));
177  temp_spec_FHC_univ_fd_ns.emplace_back(new Spectrum(loaderFDMCFHC, axis, comps[compIdx].cut, kNoShift, kMCUnivWeight));
178  }
179  spec_RHC_univ_nd.push_back(temp_spec_RHC_univ_nd);
180  spec_RHC_univ_fd_ns.push_back(temp_spec_RHC_univ_fd_ns);
181  spec_FHC_univ_nd.push_back(temp_spec_FHC_univ_nd);
182  spec_FHC_univ_fd_ns.push_back(temp_spec_FHC_univ_fd_ns);
183  }
184  spects_RHC_nd.push_back(temp_spec_RHC_nd);
185  spects_RHC_univ_nd.push_back(spec_RHC_univ_nd);
186  spects_RHC_fd_ns.push_back(temp_spec_RHC_fd_ns);
187  spects_RHC_univ_fd_ns.push_back(spec_RHC_univ_fd_ns);
188  spects_FHC_nd.push_back(temp_spec_FHC_nd);
189  spects_FHC_univ_nd.push_back(spec_FHC_univ_nd);
190  spects_FHC_fd_ns.push_back(temp_spec_FHC_fd_ns);
191  spects_FHC_univ_fd_ns.push_back(spec_FHC_univ_fd_ns);
192  }
193 
194  //Load files
195  TStopwatch gtimer;
196  gtimer.Start();
197  loaderNDMCRHC.Go();
198  loaderFDMCRHC.Go();
199  loaderNDMCFHC.Go();
200  loaderFDMCFHC.Go();
201  gtimer.Stop();
202  std::cout<<"loading time "<<gtimer.CpuTime()<<std::endl;
203 
204 
205  TFile* fout = new TFile(FileName.c_str(), "RECREATE");
206  for(int varIdx = 0; varIdx < (int)vars.size(); ++varIdx){
207  const char* var = vars[varIdx].name.c_str();
208  for(int compIdx = 0; compIdx < (int)comps.size(); ++compIdx){
209  const char* comp = comps[compIdx].name.c_str();
210  spects_RHC_nd[varIdx][compIdx]->SaveTo(fout, TString::Format("spect_nd_RHC_%s_%s", var, comp));
211  spects_RHC_fd_ns[varIdx][compIdx]->SaveTo(fout, TString::Format("spect_fd_RHC_%s_%s_ns", var, comp));
212  spects_FHC_nd[varIdx][compIdx]->SaveTo(fout, TString::Format("spect_nd_FHC_%s_%s", var, comp));
213  spects_FHC_fd_ns[varIdx][compIdx]->SaveTo(fout, TString::Format("spect_fd_FHC_%s_%s_ns", var, comp));
214 
215  for(int UnivIdx = 0; UnivIdx < N_PPFX; UnivIdx++){
216  spects_RHC_univ_nd[varIdx][compIdx][UnivIdx]->SaveTo(fout, TString::Format("spect_RHC_univ_nd_%s_%s_%d", var, comp, UnivIdx));
217  spects_RHC_univ_fd_ns[varIdx][compIdx][UnivIdx]->SaveTo(fout, TString::Format("spect_RHC_univ_fd_%s_%s_ns_%d", var, comp, UnivIdx));
218  spects_FHC_univ_nd[varIdx][compIdx][UnivIdx]->SaveTo(fout, TString::Format("spect_FHC_univ_nd_%s_%s_%d", var, comp, UnivIdx));
219  spects_FHC_univ_fd_ns[varIdx][compIdx][UnivIdx]->SaveTo(fout, TString::Format("spect_FHC_univ_fd_%s_%s_ns_%d", var, comp, UnivIdx));
220  }
221 
222  }
223 
224  }
225  delete fout;
226 
227 }
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const BeamSyst kBeamTarget((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"7mmTargetZ","+/- 7mm Target z Position")
Target z position shift +/-7mm.
Definition: BeamSysts.h:113
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const Cut kIsNumu([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==14);})
Definition: TruthCuts.h:10
const BeamSyst kBeamHornCurrent((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"2kA","2kAHornCurrent","+/- 2kA Horn Current")
Horn Current +/-2kA.
Definition: BeamSysts.h:97
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
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 abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void SetSpillCut(const SpillCut &cut)
const BeamSyst kBeamGeomWater((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmHornWater","+/- 1mm water on Horn 1")
Water layer on horn 1: +/- 1mm.
Definition: BeamSysts.h:128
const unsigned int N_PPFX
Definition: CalculateXSec.C:77
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
const BeamSyst kBeamH1PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1X","+/- 3mm Horn 1 X Position")
Horn 1 and 2 position +-3mm in X and Y separately.
Definition: BeamSysts.h:107
if(dump)
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
Struct to hold Var/plot information.
const TString name
const Cut kIsNue([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==12);})
Definition: TruthCuts.h:9
static const BeamSyst & kBeamAll
Definition: BeamSysts.h:133
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const std::map< std::pair< std::string, std::string >, Variable > vars
std::vector< std::string > comps
std::vector< float > Spectrum
Definition: Constants.h:610
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const BeamSyst kBeamPosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftY","Beam Position Y")
Definition: BeamSysts.h:104
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Var var
const Cut cut
Definition: exporter_fd.C:30
const BeamSyst kBeamH1PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1Y","+/- 3mm Horn 1 Y Position")
Definition: BeamSysts.h:108
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const BeamSyst kBeamSpotSize((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"0.2mmBeamSpotSize","0p2mmBeamSpotSize"," 1.3 +/- 0.2 mm Spot Size")
Beam Spot Size 1.3 +/- 0.2 mm both XY.
Definition: BeamSysts.h:100
const BeamSyst kBeamH2PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2X","+/- 3mm Horn 2 X Position")
Definition: BeamSysts.h:109
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
assert(nhit_max >=nhit_nbins)
const BeamSyst kBeamPosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftX","Beam Position X")
Beam position on target +-1 mm, X/Y separately.
Definition: BeamSysts.h:103
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const BeamSyst kBeamH2PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2Y","+/- 3mm Horn 2 Y Position")
Definition: BeamSysts.h:110
const BeamSyst kBeamMagField((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"Magnetic Field in Decay Pipe","MagnFieldDecayPipe","Magnetic Field in Decay Pipe")
Constant magnetic field in decay pipe.
Definition: BeamSysts.h:116
Var kMCWeight
Definition: syst_header.h:382
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
void ppfx_make_systs()
enum BeamMode string