CalculateXSec.C
Go to the documentation of this file.
1 ////////////////////////////////////////
2 /// \brief Logic to go from unfolded data (in histogram form) to a cross-section
3 /// \author Connor Johnson
4 /// \date July 2019
5 ////////////////////////////////////////
6 
7 /// CAFAna includes
8 #include "CAFAna/Vars/Vars.h"
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Core/Ratio.h"
18 
20 
21 #include "TFile.h"
22 #include "TDirectory.h"
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TH3.h"
26 
27 #include <vector>
28 #include <string>
29 
30 /// Map the variable names to their folder in the fake data (could preferably switch the labelling of the fake data instead)
31 map<string, vector<ana::Binning > > vars3D{
33 };
34 
35 map<string, string> vars1D{
36  {"ENu", "E_{#nu} (GeV)"},
37  {"Q2", "q^{2} (GeV^{2})"}
38 };
39 
40 /// Standard systmatics to check
41 const vector<string> standardSysts =
42 {
43  "cv",
44  "lightup",
45  "lightdown",
46  "calibpos",
47  "calibneg",
48  "calibshift",
49  "muoneup",
50  "muonedw",
51  "ckv",
52  "angleup",
53  "angledw",
54  "neutron_Up",
55  "neutron_Dw",
56  "neutron_nom",
57  "nom_good_seventh"
58 };
59 
60 /// Beam systematics (each has _Up and _Dw variations)
61 const vector<string> beamSysts = {
62  "2kA",
63  "02mmBeamSpotSize",
64  "1mmBeamShiftX",
65  "1mmBeamShiftY",
66  "3mmHorn1X",
67  "3mmHorn1Y",
68  "3mmHorn2X",
69  "3mmHorn2Y",
70  "7mmTargetZ",
71  "MagneticFieldinDecayPipe",
72  "1mmHornWater"
73 };
74 
75 // GENIE and PPFX
76 const unsigned int N_GENIE = 1000;
77 const unsigned int N_PPFX = 100;
78 
79 double pot = 8.09E20;
80 double N_THROWS = 1E6;
81 const string FLUX_FILE = "/pnfs/nova/persistent/users/connorj/NumuCC/flux_universes_complete.root";
82 const string FLUX_DIR = "nominalFlux";
83 
84 // Cross section functions
85 void ComputeXSec3D(string systVarName, TDirectory * varDir, vector<ana::Binning> varBins, double intFlux, double nTargs)
86 {
87  TH1 * unfoldedData = (TH1*) varDir->Get("UnfoldedData");
88  TH1 * efficiency = (TH1*) varDir->Get("Efficiency");
89  unique_ptr<TH1> effCopy(std::move((TH1*) efficiency->Clone("EffCopy")));
90  TH3 * eff3D = ana::ToTH3Helper(std::move(effCopy), varBins[0], varBins[1], varBins[2], ana::kBinContent);
91  eff3D->SetName("Efficiency3D");
92  TH1 * purity = (TH1*) varDir->Get("Purity");
93  TMatrixD * covariance = (TMatrixD*) varDir->Get("Covariance");
94 
95  TH1 * crossSection = (TH1*) unfoldedData->Clone(("xsec_" + systVarName).c_str());
96  crossSection->SetTitle((systVarName + ";cos(#theta_{#mu});T_{#mu} (GeV);EAvail (GeV);#frac{d^2#sigma}{d#cos{#theta_{#mu}}dT_{#mu}} (cm^{2} / GeV^2 / Nucleon) (Flux Avg)").c_str());
97  crossSection->Divide(eff3D);
98  crossSection->Scale(1 / intFlux);
99  crossSection->Scale(1 / nTargs);
100 
101  // Save Results
102  unfoldedData ->Write();
103  efficiency ->Write();
104  eff3D ->Write();
105  purity ->Write();
106  covariance ->Write();
107  crossSection ->Write();
108 
109  // Clear Memory
110  unfoldedData ->Delete();
111  efficiency ->Delete();
112  eff3D ->Delete();
113  purity ->Delete();
114  covariance ->Delete();
115  crossSection ->Delete();
116 }
117 
118 void ComputeXSec1D(string systVarName, TDirectory * varDir, double intFlux, double nTargs, string xAxisTitle = "")
119 {
120  TH1 * unfoldedData = (TH1*) varDir->Get("UnfoldedData");
121  TH1 * efficiency = (TH1*) varDir->Get("Efficiency");
122  TH1 * purity = (TH1*) varDir->Get("Purity");
123 
124  TH1 * crossSection = (TH1*) unfoldedData->Clone(("xsec_" + systVarName).c_str());
125  crossSection->SetTitle(systVarName.c_str());
126  crossSection->GetYaxis()->SetTitle("#sigma(cm^{2} / GeV / Nucleon) (Flux Avg)");
127  if (!xAxisTitle.empty())
128  crossSection->GetXaxis()->SetTitle(xAxisTitle.c_str());
129  crossSection->Divide(efficiency);
130  crossSection->Scale(1 / intFlux);
131  crossSection->Scale(1 / nTargs);
132 
133  // Save Results
134  unfoldedData ->Write();
135  efficiency ->Write();
136  purity ->Write();
137  crossSection ->Write();
138 
139  // Clear Memory
140  unfoldedData ->Delete();
141  efficiency ->Delete();
142  purity ->Delete();
143  crossSection ->Delete();
144 }
145 
146 void ComputeXSecENu(string systVarName, TDirectory * varDir, TH1 * fluxHist, double nTargs, string xAxisTitle = "")
147 {
148  TH1 * unfoldedData = (TH1*) varDir->Get("UnfoldedData");
149  TH1 * efficiency = (TH1*) varDir->Get("Efficiency");
150  TH1 * purity = (TH1*) varDir->Get("Purity");
151 
152  TH1 * crossSection = (TH1*) unfoldedData->Clone(("xsec_" + systVarName).c_str());
153  crossSection->SetTitle(systVarName.c_str());
154  crossSection->GetYaxis()->SetTitle("#sigma(cm^{2} / GeV / Nucleon)");
155  if (!xAxisTitle.empty())
156  crossSection->GetXaxis()->SetTitle(xAxisTitle.c_str());
157  crossSection->Divide(efficiency);
158  crossSection->Divide(fluxHist);
159  crossSection->Scale(1 / nTargs);
160 
161  // Save Results
162  unfoldedData ->Write();
163  efficiency ->Write();
164  purity ->Write();
165  crossSection ->Write();
166 
167  unfoldedData ->Delete();
168  efficiency ->Delete();
169  purity ->Delete();
170  crossSection ->Delete();
171 }
172 
173 void CalculateSyst(string syst, TFile * inputFile, TFile * outFile, TH1 * enuBinnedFlux, double intFlux, double nNucleons)
174 {
175  TDirectory * systDir = inputFile->GetDirectory(syst.c_str());
176  if (systDir && systDir != NULL){
177  TDirectory * outSystDir = outFile->mkdir(syst.c_str());
178  cout << "\n~~~~~~~~~~~~~~~~~~~\n~~~~~~~~~~~~~~~~~~~\nCalculating XSec of systematic: " << syst << endl;
179  for(pair<string, vector<ana::Binning> > varWBins : vars3D)
180  {
181  string var = varWBins.first;
182  vector<ana::Binning> varBins = varWBins.second;
183  TDirectory * varDir = systDir->GetDirectory(var.c_str());
184  if (!varDir || varDir == NULL)
185  continue;
186  TDirectory * outVarDir = outSystDir->mkdir(var.c_str());
187  outVarDir->cd();
188  ComputeXSec3D((syst + "_" + var), varDir, varBins, intFlux, nNucleons);
189  }
190 
191  for(pair<string, string> varAndTitle : vars1D)
192  {
193  string var = varAndTitle.first;
194  string varTitle = varAndTitle.second;
195  TDirectory * varDir = systDir->GetDirectory(var.c_str());
196  if (!varDir || varDir == NULL)
197  continue;
198  TDirectory * outVarDir = outSystDir->mkdir(var.c_str());
199  outVarDir->cd();
200  if (var == "ENu")
201  ComputeXSecENu((syst + "_" + var), varDir, enuBinnedFlux, nNucleons, varTitle);
202  else
203  ComputeXSec1D((syst + "_" + var), varDir, intFlux, nNucleons, varTitle);
204  }
205  }
206 }
207 
208 void CalculateXSec(string unfoldedData, string outputFile = "xsec_results.root")
209 {
210  TFile * inputFile = new TFile(unfoldedData.c_str(), "READ");
211  TFile * fluxFile = new TFile(FLUX_FILE.c_str(), "READ");
212  TFile * outFile = new TFile(outputFile.c_str(), "RECREATE");
213  outFile->cd();
214 
215  cout << "Starting target count." << endl;
217  double nNucleons = targetCount->NNucleons();
218  cout << "Number of nucleons: " << nNucleons << endl;
219 
220  // Nominal flux
221  TH1* nominalFlux = ana::Spectrum::LoadFrom(fluxFile->GetDirectory((FLUX_DIR).c_str(), true))->ToTH1(pot);
222  nominalFlux->Scale(1e-4); // From Nu/m^2 to Nu/cm^2. For result in terms of cm^2, rather than m^2.
223  nominalFlux->SetName("nominalFlux");
224  nominalFlux->Write();
225  double intFlux = nominalFlux->Integral();
226  double nuEnergyRange = nominalFlux->GetXaxis()->GetXmax() - nominalFlux->GetXaxis()->GetXmin();
227  cout << "Integrated flux over ENu (" << nuEnergyRange << " GeV) = " << intFlux << endl;
228 
229  TH1 * enuBinnedFlux = nominalFlux->Rebin(ana::enuedges.size() - 1, "enuBinnedFlux", ana::enuedges.data());
230  enuBinnedFlux->Write();
231 
232  vector<string> upDownLabel {"_Up", "_Dw"};
233 
234  for (string syst : standardSysts)
235  CalculateSyst(syst, inputFile, outFile, enuBinnedFlux, intFlux, nNucleons);
236 
237  for (string systBase : beamSysts)
238  for (string label : upDownLabel){
239  TH1* beamUnivFlux = ana::Spectrum::LoadFrom(fluxFile->GetDirectory(("beam/" + systBase + label).c_str(), true))->ToTH1(pot);
240  beamUnivFlux->Scale(1e-4);
241  double intBeamFlux = beamUnivFlux->Integral();
242  TH1 * enuBinnedBeamFlux = beamUnivFlux->Rebin(ana::enuedges.size() - 1, (systBase + label + "FineBinnedFlux").c_str(), ana::enuedges.data());
243  CalculateSyst(systBase + label, inputFile, outFile, enuBinnedBeamFlux, intBeamFlux, nNucleons);
244  }
245 
246  for(unsigned int igenie = 0; igenie < N_GENIE; ++igenie)
247  CalculateSyst("genie" + to_string(igenie), inputFile, outFile, enuBinnedFlux, intFlux, nNucleons);
248 
249  for(unsigned int ippfx = 0; ippfx < N_PPFX; ++ippfx){
250  TH1* ppfxUnivFlux = ana::Spectrum::LoadFrom(fluxFile->GetDirectory(("ppfx/ppfx" + to_string(ippfx)).c_str(), true))->ToTH1(pot);
251  ppfxUnivFlux->Scale(1e-4);
252  double intPPFXFlux = ppfxUnivFlux->Integral();
253  TH1 * enuBinnedPPFXFlux = ppfxUnivFlux->Rebin(ana::enuedges.size() - 1, ("ppfx" + to_string(ippfx) + "FineBinnedFlux").c_str(), ana::enuedges.data());
254  CalculateSyst("ppfx" + to_string(ippfx), inputFile, outFile, enuBinnedPPFXFlux, intPPFXFlux, nNucleons);
255  }
256 }
TFile * inputFile
Definition: PlotXSec.C:134
void ComputeXSecENu(string systVarName, TDirectory *varDir, TH1 *fluxHist, double nTargs, string xAxisTitle="")
map< string, string > vars1D
Definition: CalculateXSec.C:35
void CalculateSyst(string syst, TFile *inputFile, TFile *outFile, TH1 *enuBinnedFlux, double intFlux, double nNucleons)
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
map< string, vector< ana::Binning > > vars3D
Logic to go from unfolded data (in histogram form) to a cross-section.
Definition: CalculateXSec.C:31
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
void ComputeXSec1D(string systVarName, TDirectory *varDir, double intFlux, double nTargs, string xAxisTitle="")
const Binning eavailbins
void ComputeXSec3D(string systVarName, TDirectory *varDir, vector< ana::Binning > varBins, double intFlux, double nTargs)
Definition: CalculateXSec.C:85
const string FLUX_DIR
Definition: CalculateXSec.C:82
const vector< string > standardSysts
Standard systmatics to check.
Definition: CalculateXSec.C:41
const unsigned int N_PPFX
Definition: CalculateXSec.C:77
const char * label
double NNucleons() const
Number of nucleons (mass * avogadro&#39;s number)
Definition: TargetCount.h:31
const XML_Char const XML_Char * data
Definition: expat.h:268
TFile * outFile
Definition: PlotXSec.C:135
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
double pot
Definition: CalculateXSec.C:79
void CalculateXSec(string unfoldedData, string outputFile="xsec_results.root")
const unsigned int N_GENIE
Definition: CalculateXSec.C:76
Regular histogram.
Definition: UtilsExt.h:30
OStream cout
Definition: OStream.cxx:6
const Binning mukebins
Definition: NumuCCIncBins.h:90
double N_THROWS
Definition: CalculateXSec.C:80
const std::vector< double > enuedges
Definition: NumuCCIncBins.h:94
const string FLUX_FILE
Definition: CalculateXSec.C:81
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const TVector3 vtxmax(160, 160, 1000)
Float_t e
Definition: plot.C:35
const vector< string > beamSysts
Beam systematics (each has _Up and _Dw variations)
Definition: CalculateXSec.C:61
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
void efficiency()
Definition: efficiency.C:58
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10