getFitTemplates.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/ISyst.h"
8 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
15 
19 
20 using namespace ana;
21 
22 const Var kcvn ([](const caf::SRProxy* sr){
23  return (sr->sel.cvn.nueid);
24  });
25 
26 const Var kcvn2017([] (const caf::SRProxy* sr){
27  return (sr->sel.cvn2017.nueid);
28  });
29 
30 const Var kcvnProd3([] (const caf::SRProxy* sr){
31  return (sr->sel.cvnProd3Train.nueid);
32  });
33 
34 const Var klid([](const caf::SRProxy* sr){
35  return sr->sel.lid.ann;
36  });
37 
38 //sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.photonid
39 const Var kprongCVN([](const caf::SRProxy* sr){
40  if (sr->vtx.nelastic < 1) return -1000.f;
41  if (sr->vtx.elastic[0].fuzzyk.npng < 1) return -1000.f;
42 
43  int nProngs = sr->vtx.elastic[0].fuzzyk.npng;
44  float maxProngCVNe = -5.;
45  for(int iProng = 0; iProng < nProngs; iProng++){
46  if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >
47  maxProngCVNe) {
48  maxProngCVNe =
49  sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;
50  }
51  }
52  return maxProngCVNe;
53  });
54 
55 const int kNumPIDVars = 1;
57  //{"cvne", {"CVNe", Binning::Simple(40,0,1), kcvn}},
58  //{"cvne2017", {"CVNe 2017", pidbins, kcvn2017}},
59  //{"cvneProd3", {"CVNe Prod 3", pidbins, kcvnProd3}},
60  //{"lid", {"LID", pidbins, klid}},
61  {"prongCVN", {"Electron Prong Selector", pidbins, kprongCVN}},
62  //{"twoviewCVN", {"Electron Prong Selector 2 View",
63  // pidbins, ktwoviewCVN}},
64  //{"flatfluxCVN", {"Electron Prong Selector Flat Flux",
65  // pidbins, kflatfluxCVN}},
66  //{"genielikeCVN", {"Electron Prong Selector GENIE Like",
67  // pidbins, kgenielikeCVN}},
68 };
69 
70 
71 
73 {
74  std::vector<std::string> dataset_labels = {"nominal","lightdown",
75  "lightup", "ckv", "calibneg",
76  "calibpos","calibshape"};
77 
78  std::vector<std::string> datasets =
79  {
87  };
88 
89  //Loaders, Chns
90  std::vector<std::vector<ana::Spectrum*>> spec_double_templates;
91  std::vector<std::vector<ana::Spectrum*>> spec_double_xsec_templates;
92  std::vector<std::vector<ana::Spectrum*>> spec_double_ppfx_templates;
93 
94  std::vector<ana::SpectrumLoaderBase*> loaders;
95  for(uint idataset = 0; idataset < datasets.size(); ++idataset){
97  loader = new SpectrumLoader(datasets[idataset]);
99  loaders.push_back(loader);
100  }//end loop over loaders
101 
102  //Get All the Alternative flux and genie universe weights and shifts
103  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
104  GenieMultiverseParameters verse(100,
106  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
107 
108  for(uint iload = 0; iload < loaders.size(); ++iload){
109  std::vector<ana::Spectrum*> vspec;
110 
111  for(int ichn = 0; ichn < kcNumChns; ichn++){
112  vspec.push_back(new Spectrum(*loaders[iload],
113  analysis_vars[1].axis, // costheta
114  analysis_vars[2].axis, // electronE
115  pid_vars[0].axis, // prongcvn
116  kcPresel && chns[ichn].cut,
117  kNoShift,
119  }//loop over interaction channels
120 
121  //loop over all the flux and genie universes
122  //there are 100 of them so we can do it together
123  if(iload == 0){
124  std::vector<ana::Spectrum*> vspec_xsec;
125  std::vector<ana::Spectrum*> vspec_flux;
126 
127  for(int iuniv = 0; iuniv < 100; iuniv++){
128  vspec_xsec.push_back(new Spectrum(*loaders[iload],
129  analysis_vars[1].axis, //costheta
130  analysis_vars[2].axis, //electronE
131  pid_vars[0].axis, // prongcvn
132  kcPresel,
133  genie_shifts[iuniv],
135  vspec_flux.push_back(new Spectrum(*loaders[iload],
136  analysis_vars[1].axis, //costheta
137  analysis_vars[2].axis, //electronE
138  pid_vars[0].axis, //prongcnv
139  kcPresel,
140  kNoShift,
141  ppfx_univ[iuniv]* kXSecCVWgt2018));
142  }//loop over universes
143 
144  spec_double_xsec_templates.push_back(vspec_xsec);
145  spec_double_ppfx_templates.push_back(vspec_flux);
146  }//if statement
147 
148  spec_double_templates.push_back(vspec);
149  }//loader
150 
151 
152  for(uint iload = 0; iload < loaders.size(); ++iload){
153  loaders[iload]->Go();
154  }
155 
156  TFile* outf =
157  new TFile("fFitTemplates_CovarianceMatrix_newPIDBins_prongcvn.root",
158  "recreate");
159 
160  for(uint iload = 0; iload < loaders.size(); ++iload){
161  for(int ichn = 0; ichn < kcNumChns; ichn++){
162  char name[50];
163  sprintf(name, "%s_%s_%s_%s", dataset_labels[iload].c_str(),
164  pid_vars[0].name.c_str(), "double", chns[ichn].name.c_str());
165 
166  spec_double_templates[iload][ichn]->SaveTo(outf, name);
167  }//loop over interaction channels
168 
169  //loop over all the flux and genie universes
170  //there are 100 of them so we can do it together
171  if(iload == 0){
172  for(int iuniv = 0; iuniv < 100; iuniv++){
173  char name[50];
174  sprintf(name, "%s_%i_%s_%s", "genie", iuniv,
175  pid_vars[0].name.c_str(), "double");
176  spec_double_xsec_templates[iload][iuniv]->SaveTo(outf, name);
177 
178  sprintf(name, "%s_%i_%s_%s", "ppfx", iuniv,
179  pid_vars[0].name.c_str(), "double");
180  spec_double_ppfx_templates[iload][iuniv]->SaveTo(outf, name);
181  }//loop over universes
182  }//if statement
183  }//loader
184  outf->Close();
185 }
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
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const SelDef chns[knumchns]
const Cut kcPresel
Definition: NueCCIncCuts.h:139
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
std::vector< Dataset > datasets
Definition: Periods.h:3
void SetSpillCut(const SpillCut &cut)
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 std::string DWNLightFiles_full_RHC_prod4
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
const ana::Binning pidbins
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
const std::string ShapeCalibFiles_full_RHC_prod4
const Var kcvn2017([](const caf::SRProxy *sr){return(sr->sel.cvn2017.nueid);})
const mHistAxisDef analysis_vars[6]
if(dump)
const std::string UPCalibFiles_full_RHC_prod4
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
caf::Proxy< float > nueid
Definition: SRProxy.h:906
caf::StandardRecord * sr
TFile * outf
Definition: testXsec.C:51
const std::string UPLightFiles_full_RHC_prod4
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:759
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.
caf::Proxy< float > ann
Definition: SRProxy.h:971
const Var kcvn([](const caf::SRProxy *sr){return(sr->sel.cvn.nueid);})
const Cut cut
Definition: exporter_fd.C:30
const int kNumPIDVars
void getFitTemplates()
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const mHistAxisDef pid_vars[kNumPIDVars]
const int kcNumChns
Definition: NueCCIncCuts.h:293
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
const Var kcvnProd3([](const caf::SRProxy *sr){return(sr->sel.cvnProd3Train.nueid);})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const std::string DWNCalibFiles_full_RHC_prod4
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const std::string CherenkovFiles_full_RHC_prod4
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
std::string name
const Var klid([](const caf::SRProxy *sr){return sr->sel.lid.ann;})
unsigned int uint