Make_NuMuCC_Inc_XS.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Vars/Vars.h"
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Core/Ratio.h"
8 
15 
17 
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TH2.h"
21 #include "TH3.h"
22 
23 #include <iostream>
24 #include <cmath>
25 #include <string>
26 #include <iomanip>
27 
28 #include "Utilities/func/MathUtil.h"
29 
30 // directory utilities
31 #include <boost/algorithm/string/classification.hpp> // Include boost::is_any_of
32 #include <boost/algorithm/string/split.hpp> // Include for boost::split
33 #include <boost/filesystem.hpp>
34 #include <boost/filesystem/convenience.hpp>
35 
36 using namespace ana;
37 using namespace boost::filesystem;
38 
39 //const double ntargs = 3.88998e+31; //pre-calculated with updated fiducial volume(1e6 interations)
40 const double ntargs = 3.89137e+31; //pre-calculated with updated fiducial volume(1e10 interations)
41 const double pot = 8.09e20; //data POT
42 
43 const std::string fake_xs = "ppfxwgt";
44 
45 const std::string xs_mode0 = "xs_nominal";
46 const std::string xs_mode1 = "xs_bkgcons";
47 
48 const std::string dirIn = "/pnfs/nova/data/XsecAna/numucc_inc/Feb2019_results";
49 
50 const int Nuniv_ppfx = 100;
51 const int Nuniv_genie = 1000;
52 const int Nsys = 8;
53 const int Nfoc = 11;
54 const std::string sys_name[Nsys] = {"lightdown","lightup","ckv","calibneg","calibpos","calibshape","MuES_onlyMuKEDw","MuES_onlyMuKEUp"};
55 const std::string foc_name[Nfoc] = {"2kA","02mmBeamSpotSize","1mmBeamShiftX","1mmBeamShiftY","3mmHorn1X","3mmHorn1Y","3mmHorn2X","3mmHorn2Y",
56  "7mmTargetZ","MagneticFieldinDecayPipe","1mmHornWater"};
57 
58 void make_xs(TDirectory* &dirInMC,TDirectory* &dirInDt, std::string sdOut, TH1* hInFlux, TFile* fileOut, Ratio inv_pur){
59 
60  std::cout<<"=> In make_xs() "<<sdOut<<std::endl;
61 
62  ///////////////////////////
63  ///Getting all:
64  //1. Efficiency:
65  Spectrum* fMCSig = ana::LoadFrom<Spectrum> (dirInMC->GetDirectory("MCSig")).release();
66  Spectrum* fMCSigNuTree = ana::LoadFrom<Spectrum> (dirInMC->GetDirectory("MCSigNuTre")).release();
67  TH3* h3Deff = ana::ToTH3(*fMCSig , pot, kPOT, angbinsCustom, mukebins, enubins);
68  TH3* h3Deff_den = ana::ToTH3(*fMCSigNuTree, pot, kPOT, angbinsCustom, mukebins, enubins);
69  h3Deff->Divide(h3Deff_den);
70 
71  //2. Signal:
72  Spectrum* fData = ana::LoadFrom<Spectrum> (dirInDt->GetDirectory(Form("%s",fake_xs.c_str()))).release();
73  IBkgdEstimator* fBkgEst = ana::LoadFrom<IBkgdEstimator>(dirInMC->GetDirectory("BkgdEst")).release();
74  Spectrum fBkg = fBkgEst->Background();
75  Spectrum fSignalEst = (*fData) - fBkg;
76 
77  //Background constraint
78  // Ratio fSignalEst_Bkg_Contraint = fData/inv_pur;
79  Ratio fSignalEst_Bkg_Contraint = *fData/Spectrum(inv_pur.ToTH1(), fData->POT(), fData->Livetime());
80 
81  //3. Unfolding:
82  ReweightableSpectrum* fRecoTrue = ana::LoadFrom<ReweightableSpectrum>(dirInMC->GetDirectory("RecoTrue")).release();
83 
84  UnfoldIterative unfold(*fRecoTrue,5);
85  Spectrum fSigUnfolded = unfold.Truth(fSignalEst);
86 
87  UnfoldIterative unfold_bkgc(*fRecoTrue,5);
88  Spectrum fSigUnfolded_bkgc = unfold_bkgc.Truth(Spectrum(fSignalEst_Bkg_Contraint.ToTH1(), fData->POT(), fData->Livetime()));
89 
90  //4. Calculation:
91  double integrated_flux = hInFlux->Integral();
92  TH3* h3SigUnfolded[2];
93  h3SigUnfolded[0] = ana::ToTH3(fSigUnfolded , pot, kPOT, angbinsCustom, mukebins, enubins);
94  h3SigUnfolded[1] = ana::ToTH3(fSigUnfolded_bkgc, pot, kPOT, angbinsCustom, mukebins, enubins);
95 
96  TH2* h2Dxs[2];
97  TH1* h1Dxs[2];
98  for(int ii=0;ii<2;ii++){
99  h3SigUnfolded[ii]->Divide(h3Deff);
100  h2Dxs[ii] = (TH2*)h3SigUnfolded[ii]->Project3D("yx");
101  h1Dxs[ii] = (TH1*)h3SigUnfolded[ii]->Project3D("z");
102 
103  h2Dxs[ii]->Scale(1./integrated_flux);
104  for(int ibin=1;ibin<=h1Dxs[ii]->GetNbinsX();ibin++){
105  double cont_num = h1Dxs[ii]->GetBinContent(ibin);
106  double cont_den = hInFlux->GetBinContent(ibin);
107  if(cont_den>0)h1Dxs[ii]->SetBinContent(ibin,cont_num/cont_den);
108  else h1Dxs[ii]->SetBinContent(ibin,0);
109  h1Dxs[ii]->SetBinError(ibin,0);
110  }
111  h2Dxs[ii]->Scale(1./ntargs);
112  h1Dxs[ii]->Scale(1./ntargs);
113  }
114 
115  //Storing
116  //xs_nominal: calcualting the signal estomation by S-B
117  fileOut->cd();
118  fileOut->cd(xs_mode0.c_str());
119  h2Dxs[0]->Write(Form("hXS2Dsec_%s",sdOut.c_str()));
120  h1Dxs[0]->Write(Form("hXS1Dsec_%s",sdOut.c_str()));
121  //xs_bkgcons: calcualting the signal estomation by P*S
122  fileOut->cd();
123  fileOut->cd(xs_mode1.c_str());
124  h2Dxs[1]->Write(Form("hXS2Dsec_%s",sdOut.c_str()));
125  h1Dxs[1]->Write(Form("hXS1Dsec_%s",sdOut.c_str()));
126 
127  //Releasing memory:
128  delete fMCSig;
129  delete fMCSigNuTree;
130  delete fData;
131  delete fBkgEst;
132  delete fRecoTrue;
133 
134 }
135 
136 void make_nominal_xs(TFile* fpur, TFile* fsys, TH1* hInFlux)
137 {
138  std::cout <<"=> In make_nominal_xs() "<< std::endl;
139 
140  // Efficiency:
141  Spectrum* fMCSig = ana::LoadFrom<Spectrum> (fsys->GetDirectory("cv/MCSig")).release();
142  Spectrum* fMCSigNuTree = ana::LoadFrom<Spectrum> (fsys->GetDirectory("cv/MCSigNuTre")).release();
143  TH3* h3Deff = ana::ToTH3(*fMCSig , pot, kPOT, angbinsCustom, mukebins, enubins);
144  TH3* h3Deff_den = ana::ToTH3(*fMCSigNuTree, pot, kPOT, angbinsCustom, mukebins, enubins);
145  h3Deff->Divide(h3Deff_den);
146 
147  // Signal
148  Spectrum* fData = Spectrum::LoadFrom(fpur->GetDirectory("purity/MCSelected")).release();
149  IBkgdEstimator* fBkgEst = ana::LoadFrom<IBkgdEstimator>(fpur->GetDirectory("purity/BkgdEst")).release();
150  Spectrum fBkg = fBkgEst->Background();
151  Spectrum fSignalEst = (*fData) - fBkg;
152 
153  // Get unfolding matrix
154  ReweightableSpectrum* fRecoTrue = ana::LoadFrom<ReweightableSpectrum>(fsys->GetDirectory("cv/RecoTrue")).release();
155 
156  // unfold
157  UnfoldIterative unfold(*fRecoTrue,5);
158  Spectrum fSigUnfolded = unfold.Truth(fSignalEst);
159 
160  // project
161  double integrated_flux = hInFlux->Integral();
162  TH3* h3SigUnfolded = ana::ToTH3(fSigUnfolded , pot, kPOT, angbinsCustom, mukebins, enubins);
163  TH2* h2Dxs;
164  TH1* h1Dxs;
165  h3SigUnfolded->Divide(h3Deff);
166  h2Dxs = (TH2*)h3SigUnfolded->Project3D("yx");
167  h1Dxs = (TH1*)h3SigUnfolded->Project3D("z");
168 
169  h2Dxs->Scale(1./integrated_flux);
170  for(int ibin=1;ibin<=h1Dxs->GetNbinsX();ibin++){
171  double cont_num = h1Dxs->GetBinContent(ibin);
172  double cont_den = hInFlux->GetBinContent(ibin);
173  if(cont_den>0)h1Dxs->SetBinContent(ibin,cont_num/cont_den);
174  else h1Dxs->SetBinContent(ibin,0);
175  h1Dxs->SetBinError(ibin,0);
176  }
177  h2Dxs->Scale(1./ntargs);
178  h1Dxs->Scale(1./ntargs);
179 
180  // outupt
181  TFile* fileOut = new TFile("nominal_xsec.root","recreate");
182  fileOut->cd();
183  h2Dxs->Write("hXS2Dsec_nominal");
184  h1Dxs->Write("hXS1Dsec_nominal");
185 }
186 
187 void Make_NuMuCC_Inc_XS(const char* fOutname, bool nominalXsec = false){
188 
189  std::cout<<"=> start "<<std::endl;
190  //this script does the cv, det, muon energy scale and focusing systematics
191 
192  TFile* fpur = new TFile(Form("%s/purity.root" ,dirIn.c_str()),"read");
193  TFile* fsys = new TFile(Form("%s/systs_all.root" ,dirIn.c_str()),"read");
194  TFile* fmc;
195  //flux:
196  TH1* hFlux_cv = Spectrum::LoadFrom(fsys->GetDirectory("fluxes/CV"))->ToTH1(pot);
197  hFlux_cv->Scale(1e-4); // nu/m^2 to nu/cm^2
198 
199  // Read fake data
200  TFile* fflux = new TFile(Form("%s/systs_all.root" ,dirIn.c_str()),"read");
201  TH1* hFlux_ppfx[Nuniv_ppfx];
202  for(int i=0;i<Nuniv_ppfx;i++){
203  hFlux_ppfx[i] = Spectrum::LoadFrom( fflux->GetDirectory(Form("fluxes/ppfx%d",i)))->ToTH1(pot);
204  hFlux_ppfx[i]->Scale(1e-4);
205  }
206 
207  TH1* hFlux_foc[Nfoc][2]; //0:up, 1: dw
208  for(int i=0;i<Nfoc;i++){
209  hFlux_foc[i][0] = Spectrum::LoadFrom( fflux->GetDirectory(Form("fluxes/%s_Up",foc_name[i].c_str())))->ToTH1(pot);
210  hFlux_foc[i][1] = Spectrum::LoadFrom( fflux->GetDirectory(Form("fluxes/%s_Dw",foc_name[i].c_str())))->ToTH1(pot);
211  hFlux_foc[i][0]->Scale(1e-4);
212  hFlux_foc[i][1]->Scale(1e-4);
213  }
214 
215  ///Background constraint
216  Spectrum* MCSelected = Spectrum::LoadFrom(fpur->GetDirectory("purity/MCSelected")).release();
217  IBkgdEstimator* MCBkgEst = ana::LoadFrom<IBkgdEstimator>(fpur->GetDirectory("purity/BkgdEst")).release();
218  Spectrum MCBkg = MCBkgEst->Background();
219  Spectrum MCSelected_minus_MCBkd = *MCSelected - MCBkg;
220  Ratio MC_Inv_Pur = *MCSelected / MCSelected_minus_MCBkd;
221 
222  //void make_xs(TDirectory* &dirInMC,TDirectory* &dirInDt, std::string sdOut, TH1* hInFlux, TFile* fOut, Ratio inv_pur)
223  TDirectory* dir_data = fsys->GetDirectory("fakedata");
224  std::string namesyst = "";
225  TDirectory* dir_mc;
226 
227  // Nominal XSec from MC
228  if(nominalXsec){
229  make_nominal_xs(fpur, fsys, hFlux_cv);
230  return;
231  }
232 
233  //file Out
234  TFile* fOut = new TFile(fOutname, "RECREATE");
235  fOut->mkdir(xs_mode0.c_str());
236  fOut->mkdir(xs_mode1.c_str());
237 
238  //CV:
239  {
240  dir_mc = fsys->GetDirectory("cv");
241  make_xs(dir_mc,dir_data,"cv",hFlux_cv,fOut,MC_Inv_Pur);
242  }
243  //DET SYST
244  for(int i=0;i<Nsys;i++){
245  dir_mc = fsys->GetDirectory(sys_name[i].c_str());
246  make_xs(dir_mc,dir_data,sys_name[i],hFlux_cv,fOut,MC_Inv_Pur);
247  }
248  //FOC
249  for(int i=0;i<Nfoc;i++){
250  //Up
251  dir_mc = fsys->GetDirectory((foc_name[i]+"_Up").c_str());
252  make_xs(dir_mc,dir_data,foc_name[i]+"_Up",hFlux_foc[i][0],fOut,MC_Inv_Pur);
253  //Down
254  dir_mc = fsys->GetDirectory((foc_name[i]+"_Dw").c_str());
255  make_xs(dir_mc,dir_data,foc_name[i]+"_Dw",hFlux_foc[i][1],fOut,MC_Inv_Pur);
256  }
257  //PPFX
258  for(int i=0;i<Nuniv_ppfx;i++){
259  int idx_file = i%4;
260  if(idx_file==0)
261  {
262  if(boost::filesystem::exists(Form("%s/ppfx_univ_%d_to_%d.root",dirIn.c_str(),i,i+3)))
263  fmc = new TFile(Form("%s/ppfx_univ_%d_to_%d.root",dirIn.c_str(),i,i+3),"read");
264  else fmc = NULL;
265  }
266  if(fmc)
267  {
268  dir_mc = fmc->GetDirectory(Form("ppfx%d",i));
269  make_xs(dir_mc,dir_data,Form("ppfx%d",i),hFlux_ppfx[i],fOut,MC_Inv_Pur);
270  }
271  }
272  //GENIE
273  for(int i=0;i<Nuniv_genie;i++){
274  int idx_file = i%4;
275  if(idx_file==0)
276  {
277  if(boost::filesystem::exists(Form("%s/genie_univ_%d_to_%d.root",dirIn.c_str(),i,i+3)))
278  fmc = new TFile(Form("%s/genie_univ_%d_to_%d.root",dirIn.c_str(),i,i+3),"read");
279  else fmc = NULL;
280  }
281  if(fmc)
282  {
283  dir_mc = fmc->GetDirectory(Form("genie%d",i));
284  make_xs(dir_mc,dir_data,Form("genie%d",i),hFlux_cv,fOut,MC_Inv_Pur);
285  }
286  }
287 
288  fOut->Close();
289 
290 }
291 
292 
const int Nuniv_ppfx
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const double pot
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
void make_nominal_xs(TFile *fpur, TFile *fsys, TH1 *hInFlux)
const std::string dirIn
Spectrum with the value of a second variable, allowing for reweighting
const std::string foc_name[Nfoc]
TFile * fmc
Definition: bdt_com.C:14
const int Nsys
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
virtual Spectrum Background() const =0
const int Nuniv_genie
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const std::string fake_xs
void Make_NuMuCC_Inc_XS(const char *fOutname, bool nominalXsec=false)
const Binning enubins
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void make_xs(TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio inv_pur)
const double ntargs
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
virtual Spectrum Truth(const Spectrum &reco) override
const std::string xs_mode1
const Binning mukebins
Definition: NumuCCIncBins.h:90
TFile * fflux
Definition: Xdiff_gwt.C:118
const std::string xs_mode0
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
Float_t e
Definition: plot.C:35
const int Nfoc
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
Definition: UtilsExt.cxx:162
const std::string sys_name[Nsys]
enum BeamMode string