Flux.cxx
Go to the documentation of this file.
1 #include "CAFAna/XSec/Flux.h"
2 
4 #include "CAFAna/Vars/Vars.h"
8 
9 #include "TFile.h"
10 
11 #include <iostream>
12 
13 namespace ana
14 {
15  //----------------------------------------------------------------------
17  {
18  // It doesn't matter exactly what process we choose, so long as it
19  // corresponds to one of the modes GENIE considers as a top-level
20  // cross-section. Here we got for CC QE on Carbon.
21  return nu->pdg == pdg && nu->iscc &&
22  nu->mode == caf::kQE && nu->inttype == caf::kCCQE &&
23  nu->tgtZ == 6 && nu->tgtA == 12 && !nu->ischarm &&
24  nu->isvtxcont && nu->vtx.Z() < 1280;
25  // TODO nu->isvtxcont isn't a great cut here. It lets in M/C, which we
26  // don't want, and it seems like it might potentially cut out events
27  // right at the edge of the detector. Cut out M/C by hand with Z cut,
28  // that's not a great idea either.
29 
30  //(Leo Aliaga): this process does not work for very low neutrino energy (< 200 MeV)(NOvA-doc-32230-v1, slides 18-23)
31  }
32 
33  //----------------------------------------------------------------------
35  {
36  //This process works for very low neutrino energy where the muon mass becomes relevant.
37  //The total change for the integrated flux after appliying PPFX correction is about 1%
38  //(NOvA-doc-32230-v1, slides 18-23).
39  return nu->pdg == pdg && !nu->iscc &&
40  nu->mode == caf::kQE && nu->inttype == caf::kNCQE &&
41  nu->tgtZ == 6 && nu->tgtA == 12 && !nu->ischarm &&
42  nu->isvtxcont && nu->vtx.Z() < 1280 && nu->hitnuc == 2112;
43  }
44 
45  //----------------------------------------------------------------------
46  bool IsFiducial(const TVector3& nuVtx,
47  const TVector3& min,
48  const TVector3& max)
49  {
50 
51  return ( nuVtx.X() > min.X() && nuVtx.X() < max.X() &&
52  nuVtx.Y() > min.Y() && nuVtx.Y() < max.Y() &&
53  nuVtx.Z() > min.Z() && nuVtx.Z() < max.Z());
54  }
55 
56 
57  //----------------------------------------------------------------------
59  const Binning& bins,
60  int pdg,
61  const TVector3* min,
62  const TVector3* max,
63  const SystShifts& shift,
64  const NuTruthVar weight)
65  {
66  const NuTruthVar energyVar([](const caf::SRNeutrinoProxy* nu)
67  {
68  return nu->E;
69  });
70 
71  double nuclei;
72  if(min == NULL || max == NULL ){
73  static TargetCount tc(6); // Count all the carbons
74  nuclei = tc.NNucleons()/12; // #nucleons -> #nuclei
75  }
76  else{
77  static TargetCount tc(*min, *max, 1000000, 6); // Count all the carbons
78  nuclei = tc.NNucleons()/12; // #nucleons -> #nuclei
79  }
80 
81  const NuTruthVar oneOverXSecWei(
82  [pdg, nuclei](const caf::SRNeutrinoProxy* nu)
83  {
84  // GENIE uses GeV internally. We ultimately want a flux in m^-2
85  const double GeV2perm2 = 2.56819e31;
86  const double xsec = nu->xsec / GeV2perm2;
87  return 1/(xsec*nuclei);
88  }
89  );
90 
91  const NuTruthCut cut([pdg, min, max](const caf::SRNeutrinoProxy* nu)
92  {
93  if(min && max && !IsFiducial(nu->vtx, *min, *max))
94  return false;
95 
96  //We use NC QEL on Carbon as default
97  //return IsCCQEOnCarbon(nu, pdg);
98  return IsNCQEOnCarbon(nu, pdg);
99  });
100 
101  SystShifts filteredSysts;
102  for (const auto & syst : shift.ActiveSysts())
103  {
104  auto beamSyst = dynamic_cast<const BeamSyst*>(syst);
105 
106  if (!beamSyst) continue;
107  else filteredSysts.SetShift(beamSyst, shift.GetShift(syst));
108  }
109 
110  return new Spectrum(loader,
111  NuTruthHistAxis("True neutrino energy (GeV)", bins, energyVar),
112  cut, filteredSysts, oneOverXSecWei*weight);
113 
114  }
115 
116  //----------------------------------------------------------------------
117 
118 
120  const Binning& bins,
121  int pdg,
122  const TVector3* min,
123  const TVector3* max,
124  const SystShifts& shift,
125  const NuTruthVar cv_weight,
126  const std::vector< NuTruthVar>& vuniv_weight){
127 
128  std::vector<Spectrum*> vtmp_flux;
129 
130  ///Vector of multi-universe spectra:
131  unsigned int Nuniv = vuniv_weight.size();
132  for(unsigned int iuniv=0;iuniv<Nuniv;iuniv++){
133  vtmp_flux.push_back( DeriveFlux(loader,bins,pdg,min,max,shift,vuniv_weight[iuniv]) );
134  }
135 
136  ///Spectrum for the flux central value:
137  vtmp_flux.push_back( DeriveFlux(loader,bins,pdg,min,max,shift,cv_weight));
138 
139  return vtmp_flux;
140  }
141 
142  //----------------------------------------------------------------------
143  TH1D* GetFractionalError(TH1D* hcvIn, std::vector<TH1D*> vhIn){
144 
145  if( vhIn.size() == 0 ){
146  std::cerr<<"There is not universe to calculate the uncertaintity"<<std::endl;
147  exit(1);
148  }
149 
150  TH1D* hfe = (TH1D*) vhIn[0]->Clone("hfe");
151 
152  for(int ibin=1;ibin<=hfe->GetXaxis()->GetNbins();ibin++){
153 
154  //mean:
155  double mean = hcvIn->GetBinContent(ibin);
156  if(mean<=0)continue;
157  //error:
158  double err = 0.0;
159  for(unsigned int iuniv=0;iuniv<vhIn.size();iuniv++){
160  err += pow(vhIn[iuniv]->GetBinContent(ibin) - mean,2);
161  }
162  err /= double(vhIn.size());
163  err = sqrt(err);
164  err /= mean;
165  hfe->SetBinContent(ibin,err);
166  }
167 
168  return hfe;
169  }
170 
171  //----------------------------------------------------------------------
172 
174  int pdg,
175  const int Nbins,
176  double newbins[]){
177 
178  //Looking for the precalculated flux for this analysis
179  bool not_found;
180  TFile* fluxfileIn;
181  TH1D* htmp_cv;
182  std::string nupdg;
183  std::vector<TH1D*> vtmp_univ,vuniv;
184  std::vector<TH1D*> voutput; //index 0: cv with new bins, index 1: error.
185 
186  //Checking neutrino type:
187  if(pdg == 14)nupdg = "numu";
188  else if(pdg ==-14)nupdg = "numubar";
189  else if(pdg == 12)nupdg = "nue";
190  else if(pdg ==-12)nupdg = "nuebar";
191  else{
192  std::cerr<<"Check your requested pdg... only 14, -14, 12, -12 !"<<std::endl;
193  exit (1);
194  }
195 
196  //Checking input file:
197  const char* CAFAnaPrivateDir = std::getenv("SRT_PRIVATE_CONTEXT");
198  std::string fileInName = std::string(CAFAnaPrivateDir)+"/CAFAna/data/beam/flux_TA_ND_fid_"+name_ana+".root";
199  not_found = gSystem->AccessPathName(fileInName.c_str()); // true if NOT found
200  if(not_found){
201  const char* CAFAnaPublicDir = std::getenv("SRT_PUBLIC_CONTEXT");
202  fileInName = std::string(CAFAnaPublicDir) + "/CAFAna/data/beam/flux_TA_ND_fid_"+name_ana+".root";
203  not_found = gSystem->AccessPathName(fileInName.c_str());
204  }
205 
206  if(not_found){
207  std::cerr<<"No precalculated flux file for this analysis found"<<std::endl;
208  exit (1);
209  }
210 
211  //Getting histograms:
212  fluxfileIn = new TFile(fileInName.c_str(),"read");
213 
214  //Reading CV:
215  htmp_cv = (TH1D*)fluxfileIn->Get(Form("%s/cv_%s",nupdg.c_str(),nupdg.c_str()));
216 
217  //Reading Universes:
218  TDirectory *udir = fluxfileIn->GetDirectory(nupdg.c_str());
219  TIter next(udir->GetListOfKeys());
220  TKey* key;
221  while((key= (TKey*)next())){
222  TClass* cl = gROOT->GetClass(key->GetClassName());
223  if(!cl->InheritsFrom("TH1D"))continue;
224  std::string hName = key->GetName();
225  if(hName.find("univ") > 100)continue;
226  vtmp_univ.push_back((TH1D*)key->ReadObj());
227  }
228 
229  //Rebinning:
230  voutput.push_back( (TH1D*)htmp_cv->Rebin(Nbins,"hnewcv",newbins));
231  for(unsigned int i=0;i<vtmp_univ.size();i++){
232  vuniv.push_back( (TH1D*)vtmp_univ[i]->Rebin(Nbins,Form("hnewuniv%d",i),newbins));
233  }
234 
235  voutput.push_back( GetFractionalError(voutput[0],vuniv) );
236 
237  return voutput;
238 
239  }
240 
241  //----------------------------------------------------------------------
242 
243  void CalculateIntegratedFlux(std::vector<Spectrum*> flux_spectra,
244  double& integral_cv,
245  double& integral_unc,
246  int pdg,
247  const double emin,
248  const double emax){
249  //to be implemented
250  }
251 
252  //----------------------------------------------------------------------
254  double& integral_cv,
255  double& integral_unc,
256  int pdg,
257  const double emin,
258  const double emax){
259 
260  //Checking min and max energies:
261  if(emin>=emax){
262  std::cerr<<"check your energy inputs... emax should be greater emin"<<std::endl;
263  exit (1);
264  }
265 
266  double array_for_integral[2] = {emin,emax};
267  std::vector<TH1D*> vtmp = CalculateFluxAndNoCorrelatedErrorBand(name_ana,pdg,1,array_for_integral);
268 
269  integral_cv = vtmp[0]->GetBinContent(1);
270  integral_unc = vtmp[1]->GetBinContent(1);
271 
272  }
273 }
caf::Proxy< int > tgtZ
Definition: SRProxy.h:563
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
bool IsNCQEOnCarbon(const caf::SRNeutrinoProxy *nu, int pdg)
Definition: Flux.cxx:34
Beam systematics:
Definition: BeamSysts.h:37
const Var weight
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
caf::Proxy< int > tgtA
Definition: SRProxy.h:562
TH2 * Rebin(TH2 *hist, std::vector< double > edges)
OStream cerr
Definition: OStream.cxx:7
bool IsCCQEOnCarbon(const caf::SRNeutrinoProxy *nu, int pdg)
Definition: Flux.cxx:16
caf::Proxy< int > mode
Definition: SRProxy.h:543
caf::Proxy< float > E
Definition: SRProxy.h:520
caf::Proxy< int > hitnuc
Definition: SRProxy.h:532
caf::Proxy< int > inttype
Definition: SRProxy.h:534
TH1D * GetFractionalError(TH1D *hcvIn, std::vector< TH1D * > vhIn)
Definition: Flux.cxx:143
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
neutral current quasi-elastic
Definition: SREnums.h:70
const double emin
double NNucleons() const
Number of nucleons (mass * avogadro&#39;s number)
Definition: TargetCount.h:31
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:573
caf::Proxy< bool > isvtxcont
Definition: SRProxy.h:541
const double emax
bool IsFiducial(const TVector3 &nuVtx, const TVector3 &min, const TVector3 &max)
Definition: Flux.cxx:46
std::string getenv(std::string const &name)
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
T GetShift(const ISyst *syst) const
void CalculateIntegratedFlux(std::vector< Spectrum * > flux_spectra, double &integral_cv, double &integral_unc, int pdg, const double emin, const double emax)
Definition: Flux.cxx:243
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
Double_t xsec[nknots]
Definition: testXsec.C:47
charged current quasi-elastic
Definition: SREnums.h:69
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Base class for the various types of spectrum loader.
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Cut cut
Definition: exporter_fd.C:30
const int Nuniv
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:58
exit(0)
std::vector< Spectrum * > DeriveFluxCVAndUniverses(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar cv_weight, const std::vector< NuTruthVar > &vuniv_weight)
Definition: Flux.cxx:119
std::vector< const ISyst * > ActiveSysts() const
Definition: SystShifts.cxx:220
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
caf::Proxy< float > xsec
Definition: SRProxy.h:576
std::vector< TH1D * > CalculateFluxAndNoCorrelatedErrorBand(std::string name_ana, int pdg, const int Nbins, double newbins[])
Definition: Flux.cxx:173
Template for Cut and SpillCut.
Definition: Cut.h:15
_HistAxis< NuTruthVar > NuTruthHistAxis
Definition: HistAxis.h:106
Template for Var and SpillVar.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10
void next()
Definition: show_event.C:84
caf::Proxy< bool > ischarm
Definition: SRProxy.h:539
enum BeamMode string