saveSpectraForUnf.C
Go to the documentation of this file.
1 //Save spectra to do unfolding
2 
3 #ifdef __CINT__
4 void saveSpectraForUnf()
5 {
6  std::cout << "Sorry, you must run in compiled mode" << std::endl;
7 }
8 #else
9 
10 
11 
12 #include <string>
13 #include "CAFAna/Core/SpectrumLoader.h"
14 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Core/ISyst.h"
16 #include "CAFAna/Core/SystShifts.h"
17 #include "CAFAna/Systs/XSecSystLists.h"
18 #include "CAFAna/Cuts/SpillCuts.h"
19 #include "CAFAna/Cuts/TruthCuts.h"
20 
21 #include "CAFAna/Vars/GenieWeights.h"
22 #include "CAFAna/Vars/PPFXWeights.h"
23 #include "CAFAna/XSec/GenieMultiverseSyst.h"
24 #include "CAFAna/XSec/FluxMultiverseSyst.h"
25 
26 #include "NDAna/ncpi0_SemiInclusive/NDNCPi0Xsec.h"
27 #include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncCuts.h"
28 #include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncVars.h"
29 #include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncBins.h"
30 
31 #include <vector>
32 
33 using namespace ana;
34 
35 void saveSpectraForUnf()
36 {
37 
38 
39  const NuTruthVar weightST = kPPFXFluxCVWgtST*kXSecCVWgt2017ST;
40  const Var weight = VarFromNuTruthVar(weightST, 1);
41 
42  //signal up by 10% (defined in vars.h)
43  const NuTruthVar weightSigUpST = kPPFXFluxCVWgtST * kXSecCVWgt2017ST * kSigUpST;
44  const Var weightSigUp = VarFromNuTruthVar(weightSigUpST, 1);
45 
46  //signal weighted by trueKE (defined in vars.h)
47  const NuTruthVar weightSigHighST = kPPFXFluxCVWgtST * kXSecCVWgt2017ST * kSigHighWtST;
48  const Var weightSigHigh = VarFromNuTruthVar(weightSigHighST, 1);
49 
50 
51  const std::string fGENIE = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
52 
53 
54  const std::string fFakeData = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1 with stride 4 with offset 3";
55 
56 
57  SpectrumLoader lNDDA(fGENIE);
58  SpectrumLoader lNDMC(fFakeData);
59 
60  lNDDA.SetSpillCut(kStandardDQCuts && kTightBeamQualityCuts);
61  lNDMC.SetSpillCut(kStandardDQCuts && kTightBeamQualityCuts);
62 
63 
64  //Save Response Matrices for MonteCarlo,
65  Spectrum* response_unf_RecopiKE =
66  new Spectrum(lNDMC, analysis_vars_unf[0].axis,
67  analysis_vars_unf[2].axis,
68  kIsPreSel && kIsSignal,
69  kNoShift,
70  kPPFXFluxCVWgt * kXSecCVWgt2017);
71 
72  Spectrum* response_unf_RecopiCOS =
73  new Spectrum(lNDMC, analysis_vars_unf[1].axis,
74  analysis_vars_unf[3].axis,
75  kIsPreSel && kIsSignal,
76  kNoShift,
77  kPPFXFluxCVWgt * kXSecCVWgt2017);
78 
79 
80 
81  //Save spectrum for fake data where fake data is weighted by only ppfx and xseccv weights
82  Spectrum* NrmlWt_unf_RecopiKE =
83  new Spectrum(lNDDA, analysis_vars_unf[0].axis,
84  kIsPreSel && kIsSignal,
85  kNoShift,
86  weight);
87  Spectrum* NrmlWt_unf_RecopiCOS =
88  new Spectrum(lNDDA, analysis_vars_unf[1].axis,
89  kIsPreSel && kIsSignal,
90  kNoShift,
91  weight);
92  Spectrum* NrmlWt_unf_TruepiKE =
93  new Spectrum(lNDDA, analysis_vars_unf[2].axis,
94  kIsPreSel && kIsSignal,
95  kNoShift,
96  weight);
97  Spectrum* NrmlWt_unf_TruepiCOS =
98  new Spectrum(lNDDA, analysis_vars_unf[3].axis,
99  kIsPreSel && kIsSignal,
100  kNoShift,
101  weight);
102 
103  //signal up by 10%
104  Spectrum* SigUp_unf_RecopiKE =
105  new Spectrum(lNDDA, analysis_vars_unf[0].axis,
106  kIsPreSel && kIsSignal,
107  kNoShift,
108  weightSigUp);
109  Spectrum* SigUp_unf_RecopiCOS =
110  new Spectrum(lNDDA, analysis_vars_unf[1].axis,
111  kIsPreSel && kIsSignal,
112  kNoShift,
113  weightSigUp);
114  Spectrum* SigUp_unf_TruepiKE =
115  new Spectrum(lNDDA, analysis_vars_unf[2].axis,
116  kIsPreSel && kIsSignal,
117  kNoShift,
118  weightSigUp);
119  Spectrum* SigUp_unf_TruepiCOS =
120  new Spectrum(lNDDA, analysis_vars_unf[3].axis,
121  kIsPreSel && kIsSignal,
122  kNoShift,
123  weightSigUp);
124 
125 
126 
127  //Save spectrum for fake data where fake data is weighted by trueKE
128  Spectrum* SigWtdTrueKE_unf_RecopiKE =
129  new Spectrum(lNDDA, analysis_vars_unf[0].axis,
130  kIsPreSel && kIsSignal,
131  kNoShift,
132  weightSigHigh);
133  Spectrum* SigWtdTrueKE_unf_RecopiCOS =
134  new Spectrum(lNDDA, analysis_vars_unf[1].axis,
135  kIsPreSel && kIsSignal,
136  kNoShift,
137  weightSigHigh);
138  Spectrum* SigWtdTrueKE_unf_TruepiKE =
139  new Spectrum(lNDDA, analysis_vars_unf[2].axis,
140  kIsPreSel && kIsSignal,
141  kNoShift,
142  weightSigHigh);
143  Spectrum* SigWtdTrueKE_unf_TruepiCOS =
144  new Spectrum(lNDDA, analysis_vars_unf[3].axis,
145  kIsPreSel && kIsSignal,
146  kNoShift,
147  weightSigHigh);
148 
149 
150 
151  // Save Spectrum for the FakeData,
152  std::vector<const ISyst*> systs = getAllXsecSysts_2017();
153  std::vector<SystShifts> shiftsups;
154  std::vector<SystShifts> shiftsdowns;
155 
156  std::vector<Spectrum*> spects_unf_RecopiKE;
157  std::vector<Spectrum*> spects_unf_TruepiKE;
158  std::vector<Spectrum*> spects_unf_RecopiCOS;
159  std::vector<Spectrum*> spects_unf_TruepiCOS;
160 
161  // for +1 up shifts
162  for(uint i = 0; i < systs.size(); i++){
163  shiftsups.push_back(SystShifts(systs[i],1));
164 
165  spects_unf_RecopiKE.push_back(new Spectrum(lNDDA, analysis_vars_unf[0].axis,
166  kIsPreSel && kIsSignal,
167  shiftsups[i],
168  kPPFXFluxCVWgt * kXSecCVWgt2017));
169  spects_unf_RecopiCOS.push_back(new Spectrum(lNDDA, analysis_vars_unf[1].axis,
170  kIsPreSel && kIsSignal,
171  shiftsups[i],
172  kPPFXFluxCVWgt * kXSecCVWgt2017));
173  spects_unf_TruepiKE.push_back(new Spectrum(lNDDA, analysis_vars_unf[2].axis,
174  kIsPreSel && kIsSignal,
175  shiftsups[i],
176  kPPFXFluxCVWgt * kXSecCVWgt2017));
177  spects_unf_TruepiCOS.push_back(new Spectrum(lNDDA, analysis_vars_unf[3].axis,
178  kIsPreSel && kIsSignal,
179  shiftsups[i],
180  kPPFXFluxCVWgt * kXSecCVWgt2017));
181 
182  }
183 
184 
185  // for -1 sigma
186  for(uint i = 0; i < systs.size(); i++){
187  shiftsdowns.push_back(SystShifts(systs[i],-1));
188 
189  spects_unf_RecopiKE.push_back(new Spectrum(lNDDA, analysis_vars_unf[0].axis,
190  kIsPreSel && kIsSignal,
191  shiftsdowns[i],
192  kPPFXFluxCVWgt * kXSecCVWgt2017));
193  spects_unf_RecopiCOS.push_back(new Spectrum(lNDDA, analysis_vars_unf[1].axis,
194  kIsPreSel && kIsSignal,
195  shiftsdowns[i],
196  kPPFXFluxCVWgt * kXSecCVWgt2017));
197  spects_unf_TruepiKE.push_back(new Spectrum(lNDDA, analysis_vars_unf[2].axis,
198  kIsPreSel && kIsSignal,
199  shiftsdowns[i],
200  kPPFXFluxCVWgt * kXSecCVWgt2017));
201  spects_unf_TruepiCOS.push_back(new Spectrum(lNDDA, analysis_vars_unf[3].axis,
202  kIsPreSel && kIsSignal,
203  shiftsdowns[i],
204  kPPFXFluxCVWgt * kXSecCVWgt2017));
205 
206  }
207 
208 
209 
210  lNDMC.Go();
211  lNDDA.Go();
212 
213  std::string outname = "UnfoldingSpectra_S17_11_16.root";
214  TFile *outf = new TFile(outname.c_str(), "recreate");
215 
216 
217  response_unf_RecopiKE->SaveTo(outf->mkdir("response_unf_RecopiKE"));
218  response_unf_RecopiCOS->SaveTo(outf->mkdir("response_unf_RecopiCOS"));
219 
220  NrmlWt_unf_RecopiKE->SaveTo(outf->mkdir("NrmlWt_unf_RecopiKE"));
221  NrmlWt_unf_RecopiCOS->SaveTo(outf->mkdir("NrmlWt_unf_RecopiCOS"));
222  NrmlWt_unf_TruepiKE->SaveTo(outf->mkdir("NrmlWt_unf_TruepiKE"));
223  NrmlWt_unf_TruepiCOS->SaveTo(outf->mkdir("NrmlWt_unf_TruepiCOS"));
224 
225 
226  SigUp_unf_RecopiKE->SaveTo(outf->mkdir("SigUp_unf_RecopiKE"));
227  SigUp_unf_RecopiCOS->SaveTo(outf->mkdir("SigUp_unf_RecopiCOS"));
228  SigUp_unf_TruepiKE->SaveTo(outf->mkdir("SigUp_unf_TruepiKE"));
229  SigUp_unf_TruepiCOS->SaveTo(outf->mkdir("SigUp_unf_TruepiCOS"));
230 
231  SigWtdTrueKE_unf_RecopiKE->SaveTo(outf->mkdir("SigWtdTrueKE_unf_RecopiKE"));
232  SigWtdTrueKE_unf_RecopiCOS->SaveTo(outf->mkdir("SigWtdTrueKE_unf_RecopiCOS"));
233  SigWtdTrueKE_unf_TruepiKE->SaveTo(outf->mkdir("SigWtdTrueKE_unf_TruepiKE"));
234  SigWtdTrueKE_unf_TruepiCOS->SaveTo(outf->mkdir("SigWtdTrueKE_unf_TruepiCOS"));
235 
236 
237 
238  uint numSysts = systs.size();
239  for(uint i = 0; i < systs.size(); i++){
240 
241  char name_up_unf_RecopiKE[50];
242  char name_up_unf_RecopiCOS[50];
243 
244  sprintf(name_up_unf_RecopiKE, "unf_RecopiKE_up_%i",i);
245  sprintf(name_up_unf_RecopiCOS, "unf_RecopiCOS_up_%i",i);
246 
247  spects_unf_RecopiKE[i]->SaveTo(outf->mkdir(name_up_unf_RecopiKE));
248  spects_unf_RecopiCOS[i]->SaveTo(outf->mkdir(name_up_unf_RecopiCOS));
249 
250  char name_down_unf_RecopiKE[50];
251  char name_down_unf_RecopiCOS[50];
252 
253  sprintf(name_down_unf_RecopiKE, "unf_RecopiKE_down_%i",i);
254  sprintf(name_down_unf_RecopiCOS, "unf_RecopiCOS_down_%i",i);
255 
256  spects_unf_RecopiKE[i+numSysts]->SaveTo(outf->mkdir(name_down_unf_RecopiKE));
257  spects_unf_RecopiCOS[i+numSysts]->SaveTo(outf->mkdir(name_down_unf_RecopiCOS));
258 
259  }
260 
261 
262  //Save truth spectra
263 
264  for(uint i = 0; i < systs.size(); i++){
265 
266  char name_up_unf_TruepiKE[50];
267  char name_up_unf_TruepiCOS[50];
268 
269  sprintf(name_up_unf_TruepiKE, "unf_TruepiKE_up_%i",i);
270  sprintf(name_up_unf_TruepiCOS, "unf_TruepiCOS_up_%i",i);
271 
272  spects_unf_TruepiKE[i]->SaveTo(outf->mkdir(name_up_unf_TruepiKE));
273  spects_unf_TruepiCOS[i]->SaveTo(outf->mkdir(name_up_unf_TruepiCOS));
274 
275  char name_down_unf_TruepiKE[50];
276  char name_down_unf_TruepiCOS[50];
277 
278  sprintf(name_down_unf_TruepiKE, "unf_TruepiKE_down_%i",i);
279  sprintf(name_down_unf_TruepiCOS, "unf_TruepiCOS_down_%i",i);
280 
281  spects_unf_TruepiKE[i+numSysts]->SaveTo(outf->mkdir(name_down_unf_TruepiKE));
282  spects_unf_TruepiCOS[i+numSysts]->SaveTo(outf->mkdir(name_down_unf_TruepiCOS));
283 
284  }
285 
286  outf->Close();
287 }
288 
289 
290 #endif
NDNCPi0Xsec.h
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
systs
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
cmf::Spectrum
std::vector< double > Spectrum
Definition: Constants.h:746
saveSpectraForUnf
void saveSpectraForUnf()
Definition: saveSpectraForUnf.C:35
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
ana::kTightBeamQualityCuts
const SpillCut kTightBeamQualityCuts([](const caf::SRSpillProxy *s){if(s->ismc) return true; if(s->trigger==2) return true;if(s->spilltimesec==0 &&s->deltaspilltimensec==0 &&s->widthx==0) return bool(s->isgoodspill);if(std::abs(s->deltaspilltimensec) > 0.5e9) return false;if(s->spillpot< 2e12) return false;if(s->hornI< -202|| s->hornI >-198) return false;return kBeamPositionCut(s)&&kBeamWidthCut(s);})
Definition: SpillCuts.h:10
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
ana::analysis_vars_unf
const mHistAxisDef analysis_vars_unf[kEnAngVars]
Definition: NCPi0SemiIncVars.h:466
ana::kSigUpST
const NuTruthVar kSigUpST([](const caf::SRNeutrinoProxy *nu){float MassOfPi0=0.135;float kinetic=0.0f;float en=0.0f;int nbOfPi=0;int nbofprim=nu->prim.size();for(int i=0;i< nbofprim;i++){if(nu->prim[i].pdg==111){nbOfPi++;en=nu->prim[i].p.E;kinetic=en-MassOfPi0;}}if(((nu->pdg==14 &&nu->pdgorig==14)||(nu->pdg==12 &&nu->pdgorig==12))&&(nu->iscc==0)&&(kinetic >0.1)) return 1.10;else return 1.;})
ana::kStandardDQCuts
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
outf
TFile * outf
Definition: testXsec.C:51
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
ana::Spectrum::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:534
ana::kIsPreSel
const Cut kIsPreSel
Definition: ncpi0Cuts.h:1099
ana::kSigHighWtST
const NuTruthVar kSigHighWtST([](const caf::SRNeutrinoProxy *nu){float MassOfPi0=0.135;float kinetic=0.0f;float en=0.0f;int nbOfPi=0;int nbofprim=nu->prim.size();for(int i=0;i< nbofprim;i++){if(nu->prim[i].pdg==111){nbOfPi++;en=nu->prim[i].p.E;kinetic=en-MassOfPi0;}}double wt=1.1+(kinetic *0.11);if(((nu->pdg==14 &&nu->pdgorig==14)||(nu->pdg==12 &&nu->pdgorig==12))&&(nu->iscc==0)&&(kinetic >0.1)) return wt;else return 1.;})
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
om::cout
OStream cout
Definition: OStream.cxx:6
ana::VarFromNuTruthVar
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:8
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::kXSecCVWgt2017ST
const NuTruthVar kXSecCVWgt2017ST
2017 analysis master weight
Definition: XsecTunes.h:35
ana::kXSecCVWgt2017
const Var kXSecCVWgt2017
Definition: XsecTunes.h:36
ana::_Var
Template for Vars applied to any type of object.
Definition: FwdDeclare.h:12
ana::kIsSignal
const Cut kIsSignal
Definition: ncpi0Cuts.h:1094
ana::kPPFXFluxCVWgtST
const NuTruthVar kPPFXFluxCVWgtST([](const caf::SRNeutrinoProxy *nu){ if(nu->rwgt.ppfx.cv!=nu->rwgt.ppfx.cv){return 1.f;}if(nu->rwgt.ppfx.cv >90){return 1.f;}return float(nu->rwgt.ppfx.cv);})
weight events with the flux PPFX Central value correction.
Definition: PPFXWeights.h:12
ana::getAllXsecSysts_2017
std::vector< const ISyst * > getAllXsecSysts_2017()
Get master XSec syst list for 2017 analyses.
Definition: XSecSystLists.cxx:617
string
enum BeamMode string
uint
unsigned int uint
Definition: DAQChannelMapBaseOLD.h:22
Generated on Thu May 6 2021 14:10:34 for NOvA by   doxygen 1.8.11