GeniePredictionToRoot.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////
2 /// \brief Macro to conver the Genie predictions to xsec results
3 /// \author Connor Johnson
4 /// \date Jan 2020
5 //////////////////////////////////////////////////////////////////
6 
7 #include "CAFAna/Core/Spectrum.h"
10 
11 // Binning
13 
14 // ROOT includes
15 #include "TFile.h"
16 #include "TDirectory.h"
17 #include "TH1.h"
18 #include "TH2.h"
19 #include "TH3.h"
20 
21 #include <string>
22 #include <vector>
23 
24 using namespace ana;
25 
26 vector<string> genieWeights
27 {
28  "unweight",
29  "ppfx",
30  "xsec",
31  "both"
32 };
33 
34 // Bins
35 const Binning mukebins_fine = Binning::Simple(200, 0.5, 2.5);
36 const Binning q2bins_fine = Binning::Simple(500, 0.0, 5.0);
37 const Binning enubins_fine = Binning::Simple(400, 0.0, 10.0);
38 
39 double pot = 8.09E20;
40 double N_THROWS = 1E6;
41 const string FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/real_data_fluxes.root";
42 const string FLUX_DIR = "fluxes/CV";
43 const string ENU_FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/enuFineBinnedFluxes.root";
44 
45 void GeniePredictionToRoot(string filein, string fileout)
46 {
47  TFile * fin = TFile::Open(filein.c_str(), "READ");
48  TFile * fluxFile =TFile::Open(FLUX_FILE.c_str(), "READ");
49  TFile * enuFluxFile = TFile::Open(ENU_FLUX_FILE.c_str(), "READ");
50  TFile * fout = TFile::Open(fileout.c_str(), "RECREATE");
51  fout->cd();
52 
53  // Load the target count and flux info
55  double nNucleons = targetCount->NNucleons();
56  cout << "Number of nucleons: " << nNucleons << endl;
57 
58  TH1* derivedFlux = ana::Spectrum::LoadFrom(fluxFile->GetDirectory((FLUX_DIR).c_str(), true))->ToTH1(pot);
59  derivedFlux->Scale(1e-4); // From Nu/m^2 to Nu/cm^2. For result in terms of cm^2, rather than m^2.
60  derivedFlux->SetName("derivedFlux");
61  double intFlux = derivedFlux->Integral();
62  double nuEnergyRange = derivedFlux->GetXaxis()->GetXmax() - derivedFlux->GetXaxis()->GetXmin();
63  cout << "Integrated flux over ENu (" << nuEnergyRange << " GeV) = " << intFlux << endl;
64 
65  // Finer bin flux for ENu
66  TH1* enuBinnedFlux = ana::Spectrum::LoadFrom(enuFluxFile->GetDirectory("enubinnedFlux", true))->ToTH1(pot);
67  enuBinnedFlux->Scale(1e-4); // From Nu/m^2 to Nu/cm^2.
68  enuBinnedFlux->SetName("enuBinnedFlux");
69 
70  TDirectory * geniePredOld = fin->GetDirectory("geniePrediction");
71  TDirectory * geniePredNew = fout->mkdir("geniePrediction");
72  for (string weight : genieWeights)
73  {
74  TDirectory * weightDirOld = geniePredOld->GetDirectory(weight.c_str());
75  TDirectory * weightDirNew = geniePredNew->mkdir(weight.c_str());
76  weightDirNew->cd();
77 
78  // Load spectrums
79  Spectrum * eavailspec = Spectrum::LoadFrom(weightDirOld->GetDirectory("EAvail")).release();
80  Spectrum * mukinspec = Spectrum::LoadFrom(weightDirOld->GetDirectory("MuKin")).release();
81  Spectrum * enuspec = Spectrum::LoadFrom(weightDirOld->GetDirectory("ENu")).release();
82  Spectrum * q2spec = Spectrum::LoadFrom(weightDirOld->GetDirectory("Q2")).release();
83 
84  // Convert to ROOT histograms
85  std::unique_ptr<TH1> eavailhist1D(std::move(eavailspec->ToTH1(pot)));
86  std::unique_ptr<TH1> mukinhist1D(std::move(mukinspec ->ToTH1(pot)));
87  TH3 * eavailhist = ana::ToTH3Helper(std::move(eavailhist1D), ana::angbinsCustom, mukebins_fine, ana::eavailbins, ana::kBinContent);
88  TH2 * mukinhist = ana::ToTH2Helper(std::move(mukinhist1D), ana::angbinsCustom, mukebins_fine, ana::kBinContent);
89  TH1 * enuhist = enuspec->ToTH1(pot);
90  TH1 * q2hist = q2spec ->ToTH1(pot);
91 
92  // Calculate the full xsec
93  eavailhist->Scale(1.0 / (intFlux * nNucleons));
94  mukinhist->Scale(1.0 / (intFlux * nNucleons));
95  enuhist->Divide(enuBinnedFlux);
96  enuhist->Scale(1.0 / nNucleons); // Do not flux average ENu
97  q2hist->Scale(1.0 / (intFlux * nNucleons));
98 
99  // Bin-area normalize the differentials (all but width, and not the 3D until it's in 2D form)
100  mukinhist ->Scale(1, "width");
101  q2hist ->Scale(1, "width");
102 
103  // Switch ENu from sigma to sigma / ENu
104  // for (int ibinx = 1; ibinx <= enuhist->GetNbinsX(); ++ibinx)
105  // enuhist->SetBinContent(ibinx, enuhist->GetBinContent(ibinx) / enuhist->GetXaxis()->GetBinCenter(ibinx));
106 
107  // Save the results
108  eavailhist->Write("EAvail");
109  mukinhist ->Write("MuKin");
110  enuhist ->Write("ENu");
111  q2hist ->Write("Q2");
112  }
113  fout->Close();
114  fluxFile->Close();
115  fin->Close();
116 }
TString fin
Definition: Style.C:24
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
const string ENU_FLUX_FILE
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const Var weight
const TVector3 vtxmin(-130,-176, 225)
TH3 * ToTH3Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Helper for ana::ToTH3.
Definition: UtilsExt.cxx:186
void GeniePredictionToRoot(string filein, string fileout)
const Binning eavailbins
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
double NNucleons() const
Number of nucleons (mass * avogadro&#39;s number)
Definition: TargetCount.h:31
double pot
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const Binning enubins_fine
Regular histogram.
Definition: UtilsExt.h:30
OStream cout
Definition: OStream.cxx:6
vector< string > genieWeights
double N_THROWS
const string FLUX_DIR
TH2 * ToTH2Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Helper for ana::ToTH2.
Definition: UtilsExt.cxx:135
const Binning q2bins_fine
const string FLUX_FILE
const TVector3 vtxmax(160, 160, 1000)
const Binning mukebins_fine
Float_t e
Definition: plot.C:35
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10