ppfx_smooth_weights_save.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  const Cut kIsNumu(
39  [](const caf::SRProxy* sr)
40  {
41  if(sr->mc.nnu == 0) return false;
42  assert(sr->mc.nnu == 1);
43  return (abs(sr->mc.nu[0].pdg) == 14);
44  });
45 
46  const Cut kIsNue(
47  [](const caf::SRProxy* sr)
48  {
49  if(sr->mc.nnu == 0) return false;
50  assert(sr->mc.nnu == 1);
51  return (abs(sr->mc.nu[0].pdg) == 12);
52  });
53 
54  const Cut kIsRock(
55  [](const caf::SRProxy* sr)
56  {
57  if (sr->mc.nnu == 0)
58  return false;
59  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;
60  return false;
61  });
62 
63 
64  const Var kTrueE(
65  [](const caf::SRProxy* sr){
66  if(sr->mc.nnu!=1)return (float)-5.0;
67  return (float)sr->mc.nu[0].E;
68  });
69 
70 
71  std::vector<Var> kPPFXFluxUnivWgt;
72 
73  for(unsigned int UnivIdx = 0; UnivIdx < 100; UnivIdx++){
74 
75  const Var tempPPFXWgt(
76  [UnivIdx](const caf::SRProxy* sr){
77  if (!sr->hdr.ismc)return 1.f;
78  if (sr->mc.nnu != 1)return 1.f;
79  if(sr->mc.nu[0].rwgt.ppfx.vuniv[UnivIdx] <= 0)return 1.f;
80  return (float)sr->mc.nu[0].rwgt.ppfx.vuniv[UnivIdx];
81  });
82 
83  kPPFXFluxUnivWgt.push_back(tempPPFXWgt);
84  }
85 
86  const Binning bins = Binning::Simple(200, 0, 10);
87  const HistAxis TrueEAxis("True Energy (GeV)", bins, kTrueE);
88 
89  struct CompDef
90  {
92  Cut cut;
93  };
94 
95  std::vector <CompDef> comps = {
96  {"numu", kIsNumu && !kIsAntiNu},
97  {"nue", kIsNue && !kIsAntiNu},
98  {"anumu", kIsNumu && kIsAntiNu},
99  {"anue", kIsNue && kIsAntiNu},
100  };
101 
102  std::map<std::string, Spectrum*> spects;
103 
104  const std::string nearNonSFHC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
105  const std::string farNonSFHC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1";
106  const std::string nearNonSRHC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
107  const std::string farNonSRHC = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1";
108 
109  SpectrumLoader loaderNDMCFHC(nearNonSFHC);
110  SpectrumLoader loaderFDMCFHC(farNonSFHC);
111  SpectrumLoader loaderNDMCRHC(nearNonSRHC);
112  SpectrumLoader loaderFDMCRHC(farNonSRHC);
113 
114  loaderNDMCFHC.SetSpillCut(kStandardDQCuts);
115  loaderFDMCFHC.SetSpillCut(kStandardDQCuts);
116  loaderNDMCRHC.SetSpillCut(kStandardDQCuts);
117  loaderFDMCRHC.SetSpillCut(kStandardDQCuts);
118 
119  int N_PPFX = 100;
120  for(int compIdx = 0; compIdx < (int) comps.size(); ++compIdx){
121  Spectrum* temp_unweighted_ndFHC = new Spectrum(loaderNDMCFHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kUnweighted);
122  Spectrum* temp_unweighted_fdFHC = new Spectrum(loaderFDMCFHC, TrueEAxis, comps[compIdx].cut, kNoShift, kUnweighted);
123  Spectrum* temp_unweighted_ndRHC = new Spectrum(loaderNDMCRHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kUnweighted);
124  Spectrum* temp_unweighted_fdRHC = new Spectrum(loaderFDMCRHC, TrueEAxis, comps[compIdx].cut, kNoShift, kUnweighted);
125 
126  spects.insert(std::make_pair("unweighted_FHC_ND_"+comps[compIdx].name, temp_unweighted_ndFHC));
127  spects.insert(std::make_pair("unweighted_FHC_FD_"+comps[compIdx].name, temp_unweighted_fdFHC));
128  spects.insert(std::make_pair("unweighted_RHC_ND_"+comps[compIdx].name, temp_unweighted_ndRHC));
129  spects.insert(std::make_pair("unweighted_RHC_FD_"+comps[compIdx].name, temp_unweighted_fdRHC));
130 
131  Spectrum* temp_ppfx_cv_ndFHC = new Spectrum(loaderNDMCFHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxCVWgt);
132  Spectrum* temp_ppfx_cv_fdFHC = new Spectrum(loaderFDMCFHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxCVWgt);
133  Spectrum* temp_ppfx_cv_ndRHC = new Spectrum(loaderNDMCRHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxCVWgt);
134  Spectrum* temp_ppfx_cv_fdRHC = new Spectrum(loaderFDMCRHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxCVWgt);
135 
136  spects.insert(std::make_pair("ppfx_cv_FHC_ND_"+comps[compIdx].name, temp_ppfx_cv_ndFHC));
137  spects.insert(std::make_pair("ppfx_cv_FHC_FD_"+comps[compIdx].name, temp_ppfx_cv_fdFHC));
138  spects.insert(std::make_pair("ppfx_cv_RHC_ND_"+comps[compIdx].name, temp_ppfx_cv_ndRHC));
139  spects.insert(std::make_pair("ppfx_cv_RHC_FD_"+comps[compIdx].name, temp_ppfx_cv_fdRHC));
140 
141  for(int UnivIdx = 0; UnivIdx < N_PPFX; UnivIdx++){
142  const Var kPPFXFluxWgt = kPPFXFluxUnivWgt[UnivIdx];
143  Spectrum* temp_ppfx_univ_ndFHC = new Spectrum(loaderNDMCFHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxWgt);
144  Spectrum* temp_ppfx_univ_fdFHC = new Spectrum(loaderFDMCFHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxWgt);
145  Spectrum* temp_ppfx_univ_ndRHC = new Spectrum(loaderNDMCRHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxWgt);
146  Spectrum* temp_ppfx_univ_fdRHC = new Spectrum(loaderFDMCRHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxWgt);
147 
148  int dig0 = UnivIdx % 10;
149  int dig1 = (UnivIdx/10) % 10;
150  spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_FHC_ND_"+comps[compIdx].name, temp_ppfx_univ_ndFHC));
151  spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_FHC_FD_"+comps[compIdx].name, temp_ppfx_univ_fdFHC));
152  spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[compIdx].name, temp_ppfx_univ_ndRHC));
153  spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[compIdx].name, temp_ppfx_univ_fdRHC));
154 
155  }
156  }
157 
158  TStopwatch gtimer;
159  gtimer.Start();
160  loaderNDMCFHC.Go();
161  loaderFDMCFHC.Go();
162  loaderNDMCRHC.Go();
163  loaderFDMCRHC.Go();
164  gtimer.Stop();
165  std::cout<<"loading time "<<gtimer.CpuTime()<<std::endl;
166 
167 
168  std::string filename = "ppfx_weights.root";
169  TFile* fout = new TFile(filename.c_str(), "RECREATE");
170  std::for_each(spects.begin(), spects.end(),
171  [&](std::pair<std::string, Spectrum*> spect){
172  spect.second->SaveTo(fout, spect.first.c_str());
173  });
174  delete fout;
175 
176 
177 }
const XML_Char * name
Definition: expat.h:151
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2136
const Cut kIsNumu([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==14);})
Definition: TruthCuts.h:10
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:617
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:616
string filename
Definition: shutoffs.py:106
void SetSpillCut(const SpillCut &cut)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
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
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
const Cut kIsNue([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==12);})
Definition: TruthCuts.h:9
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
std::vector< std::string > comps
std::vector< float > Spectrum
Definition: Constants.h:610
void ppfx_smooth_weights_save()
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2137
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< bool > ismc
Definition: SRProxy.h:241
assert(nhit_max >=nhit_nbins)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
enum BeamMode string