getCrossSectionAnalysis_Spectra.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/ISyst.h"
7 #include "CAFAna/XSec/Flux.h"
9 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
16 
20 
21 using namespace ana;
22 
23 const Var kcvn ([](const caf::SRProxy* sr){
24  return (sr->sel.cvn.nueid);
25  });
26 
27 const Var kcvn2017([] (const caf::SRProxy* sr){
28  return (sr->sel.cvn2017.nueid);
29  });
30 
31 const Var kcvnProd3([] (const caf::SRProxy* sr){
32  return (sr->sel.cvnProd3Train.nueid);
33  });
34 
35 const Var klid([](const caf::SRProxy* sr){
36  return sr->sel.lid.ann;
37  });
38 
39 //sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.photonid
40 const Var kprongCVN([](const caf::SRProxy* sr){
41  if (sr->vtx.nelastic < 1) return -1000.f;
42  if (sr->vtx.elastic[0].fuzzyk.npng < 1) return -1000.f;
43 
44  int nProngs = sr->vtx.elastic[0].fuzzyk.npng;
45  float maxProngCVNe = -5.;
46  for(int iProng = 0; iProng < nProngs; iProng++){
47  if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >
48  maxProngCVNe) {
49  maxProngCVNe =
50  sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;
51  }
52  }
53  return maxProngCVNe;
54  });
55 
56 const int kNumPIDVars = 1;
58  //{"cvne", {"CVNe", Binning::Simple(40,0,1), kcvn}},
59  //{"cvne2017", {"CVNe 2017", pidbins, kcvn2017}},
60  //{"cvneProd3", {"CVNe Prod 3", Binning::Simple(40,0,1), kcvnProd3}},
61  //{"lid", {"LID", Binning::Simple(40,-0.1,1.1), klid}},
62  {"prongCVN", {"Electron Prong Selector", pidbins, kprongCVN}},
63  //{"twoviewCVN", {"Electron Prong Selector 2 View",
64  // Binning::Simple(40,0,1), ktwoviewCVN}},
65  //{"flatfluxCVN", {"Electron Prong Selector Flat Flux",
66  // Binning::Simple(40,0,1), kflatfluxCVN}},
67  //{"genielikeCVN", {"Electron Prong Selector GENIE Like",
68  // Binning::Simple(40,0,1), kgenielikeCVN}},
69 };
70 
71 
72 const Var kMySignalWgt([](const caf::SRProxy* sr)->float
73  {
74  if(sr->mc.nnu == 0) return 0.;
75  assert(sr->mc.nnu == 1);
76  if( (sr->mc.nu[0].iscc)
77  && (abs(sr->mc.nu[0].pdg) == 12)) return 1.06;
78  else return 1.;
79  });
80 
81 const Var kQ2Wgt([](const caf::SRProxy* sr)->float
82  {
83  if(sr->mc.nnu == 0) return 1;
84 
85  double q2 = sr->mc.nu[0].q2;
86  if(q2 < 0) return 1.0;
87  double wt = 0.8 + (q2*0.08);
88  if (wt <= 0) wt = 1.0;
89 
90  return wt;
91  });
92 
93 const NuTruthVar kMySignalWgtNT([](const caf::SRNeutrinoProxy* sr)->float
94  {
95  if( (sr->iscc)
96  && (abs(sr->pdg) == 12)) return 1.06;
97  else return 1.;
98  });
99 
100 const NuTruthVar kQ2WgtNT([](const caf::SRNeutrinoProxy* sr)->float
101  {
102  double q2 = sr->q2;
103  if(q2 < 0) return 1.0;
104  double wt = 0.8 + (q2*0.08);
105  if (wt <= 0) wt = 1.0;
106 
107  return wt;
108  });
109 
111 {
112  std::vector<std::string> dataset_labels = {"nominal","fakedata",
113  "fakedata_q2", "fakedata_sig"};
114 
115  std::vector<std::string> datasets =
116  {
117  NominalMC_entire_RHC_prod4, NominalMC_entire_RHC_prod4 //, NominalFakeData_2_full_RHC
118  };
119 
120  std::vector<Var> weights = { (kPPFXFluxCVWgt * kXSecCVWgt2018),
122  (kPPFXFluxCVWgt * kXSecCVWgt2018 * kQ2Wgt),
123  (kPPFXFluxCVWgt * kXSecCVWgt2018 * kMySignalWgt)
124  };
125 
126  std::vector<NuTruthVar> weights_nt = { (kPPFXFluxCVWgtST *
128  (kPPFXFluxCVWgtST *
130  (kPPFXFluxCVWgtST *
131  kXSecCVWgt2018_smallerDISScale_NT *
132  kQ2WgtNT),
133  (kPPFXFluxCVWgtST *
134  kXSecCVWgt2018_smallerDISScale_NT *
136  };
137 
138 
139  std::vector<ana::SpectrumLoaderBase*> loaders;
140  for(uint idataset = 0; idataset < 4; ++idataset){
142  uint holder = idataset;
143  if(idataset > 0) idataset = 1;
144  loader = new SpectrumLoader(datasets[idataset]);
146  loaders.push_back(loader);
147  idataset = holder;
148  }//end loop over loaders
149 
150  std::vector<std::vector<ana::Spectrum*>> spec_analysis_3d;
151  std::vector<ana::Spectrum*> spec_template;
152  std::vector<ana::Spectrum*> spec_unfolding_1d;
153  std::vector<ana::Spectrum*> spec_unfolding_2d;
154  std::vector<ana::Spectrum*> spec_efficiency_num_1d;
155  std::vector<ana::Spectrum*> spec_efficiency_num_2d;
156  std::vector<ana::Spectrum*> spec_efficiency_denom_1d;
157  std::vector<ana::Spectrum*> spec_efficiency_denom_2d;
158 
159  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
160  std::vector<ana::NuTruthVar> ppfx_univ_st = ana::GetkPPFXFluxUnivWgtST();
161 
162  for(uint iload = 0; iload < loaders.size(); ++iload){
163 
164  std::vector<ana::Spectrum*> vspec_analysis_3d;
165  for(int ichn = 0; ichn < kcNumChns; ichn++){
166  vspec_analysis_3d.push_back(new Spectrum(*loaders[iload],
167  analysis_vars[1].axis, // costheta
168  analysis_vars[2].axis, // electron E
169  analysis_vars[0].axis, // recoE
170  kcPresel && chns[ichn].cut,
171  kNoShift,
172  weights[iload]));
173  }//loop over interaction channels
174  spec_analysis_3d.push_back(vspec_analysis_3d);
175 
176 
177  spec_template.push_back(new Spectrum(*loaders[iload],
178  analysis_vars[1].axis, // costheta
179  analysis_vars[2].axis, // electron E
180  pid_vars[0].axis, // prongCVN
181  kcPresel,
182  kNoShift,
183  weights[iload]));
184 
185  spec_unfolding_1d.push_back(new Spectrum(*loaders[iload],
186  analysis_vars[0].axis, // recoE
187  analysis_vars[5].axis, // true E
188  kcPresel &&
190  kNoShift,
191  weights[iload]));
192  spec_unfolding_2d.push_back(new Spectrum(*loaders[iload],
195  kcPresel &&
197  kNoShift,
198  weights[iload]));
199  spec_efficiency_num_1d.push_back(new Spectrum(*loaders[iload],
201  kcPresel &&
203  kNoShift,
204  weights[iload]));
205  spec_efficiency_num_2d.push_back(new Spectrum(*loaders[iload],
207  kcPresel &&
209  kNoShift,
210  weights[iload]));
211  spec_efficiency_denom_1d.push_back(new Spectrum(*loaders[iload],
212  analysis_vars_truth[0].axis,
213  kIsTrueSigST,
214  kNoShift,
215  weights_nt[iload]));
216  spec_efficiency_denom_2d.push_back(new Spectrum(*loaders[iload],
218  kIsTrueSigST,
219  kNoShift,
220  weights_nt[iload]));
221  }
222 
223  for(uint iload = 0; iload < loaders.size(); ++iload){
224  loaders[iload]->Go();
225  }
226 
227 
228  TFile* outf =
229  new TFile("fCrossSection_Nominal_FakeData_prongcvn_new.root", "recreate");
230 
231  for(uint iload = 0; iload < loaders.size(); ++iload){
232  for(int ichn = 0; ichn < kcNumChns; ichn++){
233  char name[50];
234  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
235  "analysis", chns[ichn].name.c_str());
236  spec_analysis_3d[iload][ichn]->SaveTo(outf, name);
237  }//loop over interaction channels
238 
239  char name[50];
240  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
241  "template", "tot");
242  spec_template[iload]->SaveTo(outf, name);
243  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
244  "unfolding", "1d");
245  spec_unfolding_1d[iload]->SaveTo(outf, name);
246  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
247  "unfolding", "2d");
248  spec_unfolding_2d[iload]->SaveTo(outf, name);
249  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
250  "efficiency", "num_1d");
251  spec_efficiency_num_1d[iload]->SaveTo(outf, name);
252  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
253  "efficiency", "num_2d");
254  spec_efficiency_num_2d[iload]->SaveTo(outf, name);
255  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
256  "efficiency", "denom_1d");
257  spec_efficiency_denom_1d[iload]->SaveTo(outf, name);
258  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
259  "efficiency", "denom_2d");
260  spec_efficiency_denom_2d[iload]->SaveTo(outf, name);
261  }//loader
262 
263  outf->Close();
264 
265 
266 
267 
268 }
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
const XML_Char * name
Definition: expat.h:151
const std::string NominalMC_entire_RHC_prod4
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
const NuTruthCut kIsTrueSigST
const Var klid([](const caf::SRProxy *sr){return sr->sel.lid.ann;})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kcvn([](const caf::SRProxy *sr){return(sr->sel.cvn.nueid);})
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const SelDef chns[knumchns]
const Cut kcPresel
Definition: NueCCIncCuts.h:139
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const Var kcvnProd3([](const caf::SRProxy *sr){return(sr->sel.cvnProd3Train.nueid);})
std::vector< Dataset > datasets
Definition: Periods.h:3
const mHistAxisSTDef analysis_vars_truth[kvAnalysisVarsST]
Definition: NueCCIncVars.h:489
void abs(TH1 *hist)
const mHistAxisDef pid_vars[kNumPIDVars]
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void SetSpillCut(const SpillCut &cut)
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
const ana::Binning pidbins
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
Double_t q2[12][num]
Definition: f2_nu.C:137
void getCrossSectionAnalysis_Spectra()
const Var kcvn2017([](const caf::SRProxy *sr){return(sr->sel.cvn2017.nueid);})
const mHistAxisDef analysis_vars[6]
if(dump)
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
caf::Proxy< float > nueid
Definition: SRProxy.h:906
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
Var weights
caf::StandardRecord * sr
HistAxis HistAxisFromNuTruthHistAxis(NuTruthHistAxis ntha, double _default)
Definition: HistAxis.cxx:9
TFile * outf
Definition: testXsec.C:51
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:47
loader
Definition: demo0.py:10
const HistAxis kTrueElecEVsCosStandardAxis
Definition: NueCCIncVars.h:533
std::vector< float > Spectrum
Definition: Constants.h:759
const int kNumPIDVars
const Var kprongCVN([](const caf::SRProxy *sr){if(sr->vtx.nelastic< 1) return-1000.f;if(sr->vtx.elastic[0].fuzzyk.npng< 1) return-1000.f;int nProngs=sr->vtx.elastic[0].fuzzyk.npng;float maxProngCVNe=-5.;for(int iProng=0;iProng< nProngs;iProng++){if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >maxProngCVNe){maxProngCVNe=sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;}}return maxProngCVNe;})
const Var kQ2Wgt([](const caf::SRProxy *sr) ->float{if(sr->mc.nnu==0) return 1;double q2=sr->mc.nu[0].q2;if(q2< 0) return 1.0;double wt=0.8+(q2 *0.08);if(wt<=0) wt=1.0;return wt;})
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< caf::SRELid > lid
Definition: SRProxy.h:1261
Base class for the various types of spectrum loader.
const HistAxis kRecoElecEVsCosStandardAxis("Reconstructed E_{e} vs cos #theta; Elec_{e} (GeV); cos #{theta}", double_diff_bins, kRecoElecEVsCos)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< float > ann
Definition: SRProxy.h:971
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const int kcNumChns
Definition: NueCCIncCuts.h:293
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< Loaders * > loaders
Definition: syst_header.h:386
assert(nhit_max >=nhit_nbins)
const NuTruthVar kQ2WgtNT([](const caf::SRNeutrinoProxy *sr) ->float{double q2=sr->q2;if(q2< 0) return 1.0;double wt=0.8+(q2 *0.08);if(wt<=0) wt=1.0;return wt;})
caf::Proxy< float > q2
Definition: SRProxy.h:557
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const NuTruthHistAxis kTrueElecEVsCosStandardAxisST("True E_{e} vs cos #theta;cs Elec_{e} (GeV); cos #{theta}", double_diff_bins, kTrueElecEVsCosST)
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var kMySignalWgt([](const caf::SRProxy *sr) ->float{if(sr->mc.nnu==0) return 0.;assert(sr->mc.nnu==1);if((sr->mc.nu[0].iscc) &&(abs(sr->mc.nu[0].pdg)==12)) return 1.06;else return 1.;})
const NuTruthVar kMySignalWgtNT([](const caf::SRNeutrinoProxy *sr) ->float{if((sr->iscc) &&(abs(sr->pdg)==12)) return 1.06;else return 1.;})
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
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
unsigned int uint