Functions | Variables
GeniePredictionToRoot.C File Reference
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/XSec/TargetCount.h"
#include "NDAna/numucc_inc/NumuCCIncCuts.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include <string>
#include <vector>

Go to the source code of this file.

Functions

void GeniePredictionToRoot (string filein, string fileout)
 

Variables

vector< stringgenieWeights
 
const Binning mukebins_fine = Binning::Simple(200, 0.5, 2.5)
 
const Binning q2bins_fine = Binning::Simple(500, 0.0, 5.0)
 
const Binning enubins_fine = Binning::Simple(400, 0.0, 10.0)
 
double pot = 8.09E20
 
double N_THROWS = 1E6
 
const string FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/real_data_fluxes.root"
 
const string FLUX_DIR = "fluxes/CV"
 
const string ENU_FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/enuFineBinnedFluxes.root"
 

Function Documentation

void GeniePredictionToRoot ( string  filein,
string  fileout 
)

Definition at line 45 of file GeniePredictionToRoot.C.

References ana::angbinsCustom, om::cout, e, ana::eavailbins, allTimeWatchdog::endl, ENU_FLUX_FILE, fin, FLUX_DIR, FLUX_FILE, submit_syst::fout, genieWeights, ana::kBinContent, ana::Spectrum::LoadFrom(), N_THROWS, ana::TargetCount::NNucleons(), pot, runNovaSAM::release, ana::Spectrum::ToTH1(), ana::ToTH2Helper(), ana::ToTH3Helper(), ana::vtxmax(), ana::vtxmin(), and ana::weight.

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
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
const Binning eavailbins
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
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
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 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
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10

Variable Documentation

const string ENU_FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/enuFineBinnedFluxes.root"

Definition at line 43 of file GeniePredictionToRoot.C.

Referenced by GeniePredictionToRoot().

const Binning enubins_fine = Binning::Simple(400, 0.0, 10.0)

Definition at line 37 of file GeniePredictionToRoot.C.

const string FLUX_DIR = "fluxes/CV"

Definition at line 42 of file GeniePredictionToRoot.C.

Referenced by GeniePredictionToRoot().

const string FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/real_data_fluxes.root"

Definition at line 41 of file GeniePredictionToRoot.C.

Referenced by GeniePredictionToRoot().

vector<string> genieWeights
Initial value:
{
"unweight",
"ppfx",
"xsec",
"both"
}

Definition at line 27 of file GeniePredictionToRoot.C.

Referenced by GeniePredictionToRoot().

const Binning mukebins_fine = Binning::Simple(200, 0.5, 2.5)

Definition at line 35 of file GeniePredictionToRoot.C.

double N_THROWS = 1E6

Definition at line 40 of file GeniePredictionToRoot.C.

Referenced by GeniePredictionToRoot().

double pot = 8.09E20

Definition at line 39 of file GeniePredictionToRoot.C.

Referenced by GeniePredictionToRoot().

const Binning q2bins_fine = Binning::Simple(500, 0.0, 5.0)

Definition at line 36 of file GeniePredictionToRoot.C.