Functions | Variables
Make_NuMuCC_Inc_XS_v2.C File Reference
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Ratio.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "NDAna/numucc_inc/NumuCCIncAnalysis.h"
#include "NDAna/numucc_inc/NumuCCIncBins.h"
#include "NDAna/numucc_inc/NumuCCIncVars.h"
#include "NDAna/numucc_inc/NumuCCIncCuts.h"
#include "NDAna/numucc_inc/NDXSecMuonPID.h"
#include "CAFAna/Unfold/UnfoldIterative.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include <iostream>
#include <cmath>
#include <string>
#include <iomanip>
#include <map>
#include "Utilities/func/MathUtil.h"
#include <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/filesystem.hpp>
#include <boost/filesystem/convenience.hpp>

Go to the source code of this file.

Functions

void make_xs (TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio *inv_pur)
 
void make_xs_1D (TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio *inv_pur)
 
void Make_NuMuCC_Inc_XS_v2 (const std::string inputDir, const std::string fOutname)
 

Variables

const double ntargs = 3.89137e+31
 
const double pot = 8.09e20
 
const std::string fake_xs = "ppfxwgt"
 
const std::string xs_mode0 = "xs_nominal"
 
const std::string xs_mode1 = "xs_bkgcons"
 
std::string dirIn = "input_files"
 
const int Nuniv_ppfx = 0
 
const int Nuniv_genie = 0
 
std::map< std::string, std::stringcoreFiles
 
const std::vector< std::stringunfoldVars3D
 
const std::vector< std::stringunfoldVars1D
 
std::map< std::string, std::stringvar2FD
 
const std::map< std::string, std::stringsystFiles
 
const std::map< std::string, std::stringbeamSysts
 
const std::vector< std::stringupdownlabels = {"_Up", "_Dw"}
 

Function Documentation

void Make_NuMuCC_Inc_XS_v2 ( const std::string  inputDir,
const std::string  fOutname 
)

Definition at line 315 of file Make_NuMuCC_Inc_XS_v2.C.

References ana::assert(), beamSysts, coreFiles, om::cout, dirIn, e, allTimeWatchdog::endl, fmc, makeBrightnessMap::fOut, MECModelEnuComparisons::i, ana::Spectrum::LoadFrom(), make_xs(), make_xs_1D(), Nuniv_genie, Nuniv_ppfx, pot, string, systFiles, art::to_string(), unfoldVars1D, unfoldVars3D, updownlabels, PandAna.Demos.tute_pid_validation::var, var2FD, xs_mode0, and xs_mode1.

315  {
316  ::dirIn = inputDir;
317  std::cout<<"=> start "<<std::endl;
318  //this script does the cv, det, muon energy scale and focusing systematics
319 
320  // Load core files
321  TFile* fsys = new TFile((dirIn + "/" + coreFiles["nominal"]).c_str());
322  TFile* fluxFile = new TFile((dirIn + "/" + coreFiles["fluxes"]).c_str());
323  TFile* fakedataFile = new TFile((dirIn + "/" + coreFiles["fakedata"]).c_str());
324  TFile* fmc;
325  assert(fsys && fsys->IsOpen());
326  assert(fluxFile && fluxFile->IsOpen());
327  assert(fakedataFile && fakedataFile->IsOpen());
328 
329  //flux:
330  TH1* hFlux_cv = Spectrum::LoadFrom(fluxFile->GetDirectory("fluxes/CV"))->ToTH1(pot);
331  hFlux_cv->Scale(1e-4); // nu/m^2 to nu/cm^2
332 
333  // Read ppfx
334  TH1* hFlux_ppfx[Nuniv_ppfx];
335  for(int i=0;i<Nuniv_ppfx;i++){
336  hFlux_ppfx[i] = Spectrum::LoadFrom(fluxFile->GetDirectory(Form("fluxes/ppfx%d",i)))->ToTH1(pot);
337  hFlux_ppfx[i]->Scale(1e-4);
338  }
339 
340  //file Out
341  TFile* fOut = new TFile(fOutname.c_str(), "RECREATE");
342  fOut->mkdir(xs_mode0.c_str());
343  fOut->mkdir(xs_mode1.c_str());
344 
345  // Load all the unfolding spectrums we'll need. Load 3D, 1D separately to keep track.
346  std::map<std::string, TDirectory*> dir_data;
347  std::map<std::string, TDirectory*> nuTrueDirs;
348  std::map<std::string, Spectrum*> MCSelected;
349  std::map<std::string, Spectrum*> MCSignal;
350  std::map<std::string, Spectrum*> MCSignal_Reco;
351  std::map<std::string, Spectrum*> MCBkg;
352  std::map<std::string, Ratio*> MC_Inv_Pur;
353  for (std::string var : unfoldVars3D){
354  dir_data[var] = fakedataFile->GetDirectory(var2FD[var].c_str(), true);
355  assert(dir_data[var] && dir_data[var] != NULL);
356  MCSelected[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSel_Reco").c_str(), true)).release();
357  MCSignal_Reco[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSig_Reco").c_str(), true)).release();
358  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTrue").c_str(), true);
359  if (!nuTrueDirs[var] || nuTrueDirs[var] == NULL || nuTrueDirs[var] == 0)
360  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTre").c_str(), true); // Fix typo in dir name
361  assert(nuTrueDirs[var] && nuTrueDirs[var] != NULL);
362  MCSignal[var] = Spectrum::LoadFrom(nuTrueDirs[var]).release();
363  MC_Inv_Pur[var] = new Ratio(*MCSelected[var], *MCSignal_Reco[var]);
364  }
365  for (std::string var : unfoldVars1D){
366  dir_data[var] = fakedataFile->GetDirectory(var2FD[var].c_str(), true);
367  assert(dir_data[var] && dir_data[var] != NULL);
368  MCSelected[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSel_Reco").c_str(), true)).release();
369  MCSignal_Reco[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSig_Reco").c_str(), true)).release();
370  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTrue").c_str(), true);
371  if (!nuTrueDirs[var] || nuTrueDirs[var] == NULL || nuTrueDirs[var] == 0)
372  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTre").c_str(), true); // Fix typo in dir name
373  assert(nuTrueDirs[var] && nuTrueDirs[var] != NULL);
374  MCSignal[var] = Spectrum::LoadFrom(nuTrueDirs[var]).release();
375  MC_Inv_Pur[var] = new Ratio(*MCSelected[var], *MCSignal_Reco[var]);
376  }
377 
378  // Holder for the current var directory
379  TDirectory* dir_mc;
380  //CV:
381  for(std::string var : unfoldVars3D)
382  {
383  dir_mc = fsys->GetDirectory(("cv/" + var).c_str(), true);
384  assert(dir_mc && dir_mc != NULL);
385  make_xs(dir_mc, dir_data[var], "cv_" + var, hFlux_cv, fOut, MC_Inv_Pur[var]);
386  }
387 
388  for(std::string var : unfoldVars1D)
389  {
390  dir_mc = fsys->GetDirectory(("cv/" + var).c_str(), true);
391  assert(dir_mc && dir_mc != NULL);
392  make_xs_1D(dir_mc, dir_data[var], "cv_" + var, hFlux_cv, fOut, MC_Inv_Pur[var]);
393  }
394 
395  // Systematic files and Muon Energy
396  for (std::pair<std::string, std::string> syst : systFiles){
397  std::string systName = syst.first;
398  std::string systFilename = syst.second;
399  TFile * systFile = new TFile((dirIn + "/" + systFilename).c_str());
400  for(std::string var : unfoldVars3D)
401  {
402  dir_mc = systFile->GetDirectory((systName + "/" + var).c_str(), true);
403  assert(dir_mc && dir_mc != NULL);
404  make_xs(dir_mc, dir_data[var], (systName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
405  }
406 
407  for(std::string var : unfoldVars1D)
408  {
409  dir_mc = systFile->GetDirectory((systName + "/" + var).c_str(), true);
410  assert(dir_mc && dir_mc != NULL);
411  make_xs_1D(dir_mc, dir_data[var], (systName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
412  }
413  systFile->Close();
414  }
415 
416  // Beam Systematics
417  for (std::pair<std::string, std::string> beamsyst : beamSysts){
418  std::string beamsystName = beamsyst.first;
419  std::string beamFilename = beamsyst.second;
420  TFile * beamsystFile = new TFile((dirIn + "/" + beamFilename).c_str());
421 
422  std::string beamsystupName = beamsystName + updownlabels[0];
423  std::string beamsystdwName = beamsystName + updownlabels[1];
424  for(std::string var : unfoldVars3D)
425  {
426  dir_mc = beamsystFile->GetDirectory((beamsystupName + "/" + var).c_str(), true);
427  assert(dir_mc && dir_mc != NULL);
428  make_xs(dir_mc, dir_data[var], (beamsystupName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
429 
430  dir_mc = beamsystFile->GetDirectory((beamsystdwName + "/" + var).c_str(), true);
431  assert(dir_mc && dir_mc != NULL);
432  make_xs(dir_mc, dir_data[var], (beamsystdwName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
433  }
434 
435  for(std::string var : unfoldVars1D)
436  {
437  dir_mc = beamsystFile->GetDirectory((beamsystupName + "/" + var).c_str(), true);
438  assert(dir_mc && dir_mc != NULL);
439  make_xs_1D(dir_mc, dir_data[var], (beamsystupName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
440 
441  dir_mc = beamsystFile->GetDirectory((beamsystdwName + "/" + var).c_str(), true);
442  assert(dir_mc && dir_mc != NULL);
443  make_xs_1D(dir_mc, dir_data[var], (beamsystdwName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
444  }
445 
446  beamsystFile->Close();
447  }
448 
449  //PPFX
450  TFile * ppfxFile;
451  for(int i=0;i<Nuniv_ppfx;i++){
452  // Only update file every 4 univs (Please can we do this in a smarter way?)
453  int idx_file = i%4;
454  if(idx_file==0)
455  {
456  // Close and delete file if it exists
457  if (ppfxFile && ppfxFile != NULL){
458  ppfxFile->Close();
459  delete ppfxFile;
460  ppfxFile = NULL;
461  }
462 
463  // Replace the file
464  std::string ppfxFilename = dirIn + "/ppfx_univ_" + to_string(i) + "_to_" + to_string(i+3) + ".root";
465  if(boost::filesystem::exists(ppfxFilename.c_str()))
466  ppfxFile = new TFile(ppfxFilename.c_str(), "read");
467 
468  }
469  if (!ppfxFile || ppfxFile == NULL || ppfxFile->IsZombie())
470  continue;
471 
472  // Perform unfolding for all 3D and 1D vars included for this PPFX
473  std::string ppfxLabel = "ppfx" + to_string(i);
474  for(std::string var : unfoldVars3D)
475  {
476  dir_mc = ppfxFile->GetDirectory((ppfxLabel + "/" + var).c_str(), true);
477  assert(dir_mc && dir_mc != NULL);
478  make_xs(dir_mc, dir_data[var], (ppfxLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
479  }
480 
481  for(std::string var : unfoldVars1D)
482  {
483  dir_mc = ppfxFile->GetDirectory((ppfxLabel + "/" + var).c_str(), true);
484  assert(dir_mc && dir_mc != NULL);
485  make_xs_1D(dir_mc, dir_data[var], (ppfxLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
486  }
487  }
488  if (ppfxFile && ppfxFile != NULL){
489  ppfxFile->Close();
490  delete ppfxFile;
491  ppfxFile = NULL;
492  }
493 
494  //GENIE
495  TFile * genieFile;
496  for(int i=0;i<Nuniv_genie;i++){
497  int idx_file = i%4;
498  if(idx_file==0)
499  {
500  // Close and delete file if it exists
501  if (genieFile && genieFile != NULL){
502  genieFile->Close();
503  delete genieFile;
504  genieFile = NULL;
505  }
506 
507  // Replace the file
508  std::string genieFilename = dirIn + "/genie_univ_" + to_string(i) + "_to_" + to_string(i+3) + ".root";
509  if(boost::filesystem::exists(genieFilename.c_str()))
510  genieFile = new TFile(genieFilename.c_str(), "read");
511 
512  }
513  if (!genieFile || genieFile == NULL || genieFile->IsZombie())
514  continue;
515 
516  std::string genieLabel = "genie" + to_string(i);
517  for(std::string var : unfoldVars3D)
518  {
519  dir_mc = genieFile->GetDirectory((genieLabel + "/" + var).c_str(), true);
520  assert(dir_mc && dir_mc != NULL);
521  make_xs(dir_mc, dir_data[var], (genieLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
522  }
523 
524  for(std::string var : unfoldVars1D)
525  {
526  dir_mc = genieFile->GetDirectory((genieLabel + "/" + var).c_str(), true);
527  assert(dir_mc && dir_mc != NULL);
528  make_xs_1D(dir_mc, dir_data[var], (genieLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
529  }
530  }
531  if (genieFile && genieFile != NULL){
532  genieFile->Close();
533  delete genieFile;
534  genieFile = NULL;
535  }
536 
537  fOut->Close();
538 }
const int Nuniv_genie
const int Nuniv_ppfx
TFile * fmc
Definition: bdt_com.C:14
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
const std::map< std::string, std::string > systFiles
const std::string xs_mode0
const std::vector< std::string > unfoldVars1D
const std::map< std::string, std::string > beamSysts
const std::string xs_mode1
std::string dirIn
const std::vector< std::string > updownlabels
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
std::map< std::string, std::string > var2FD
const double pot
void make_xs_1D(TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio *inv_pur)
assert(nhit_max >=nhit_nbins)
std::map< std::string, std::string > coreFiles
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35
const std::vector< std::string > unfoldVars3D
void make_xs(TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio *inv_pur)
enum BeamMode string
void make_xs ( TDirectory *&  dirInMC,
TDirectory *&  dirInDt,
std::string  sdOut,
TH1 *  hInFlux,
TFile *  fileOut,
Ratio inv_pur 
)

Getting all:

Definition at line 108 of file Make_NuMuCC_Inc_XS_v2.C.

References ana::angbinsCustom, ana::IBkgdEstimator::Background(), om::cout, ana::eavailbins, allTimeWatchdog::endl, fake_xs, make_syst_table_plots::ibin, ana::kPOT, ana::Spectrum::Livetime(), ana::LoadFrom< IBkgdEstimator >(), ana::mukebins, ntargs, pot, ana::Spectrum::POT(), runNovaSAM::release, ana::Spectrum::Scale(), ana::Ratio::ToTH1(), ana::ToTH3(), ana::UnfoldIterative::Truth(), xs_mode0, and xs_mode1.

Referenced by Make_NuMuCC_Inc_XS_v2().

108  {
109 
110  std::cout<<"=> In make_xs() "<<sdOut<<std::endl;
111 
112  ///////////////////////////
113  ///Getting all:
114  //1. Efficiency:
115  TDirectory * nuTrueDir = dirInMC->GetDirectory("MCSigNuTrue");
116  if (!nuTrueDir || nuTrueDir == NULL || nuTrueDir == 0)
117  nuTrueDir = dirInMC->GetDirectory("MCSigNuTre"); // Try both titles to fix typo
118  Spectrum* fMCSig = ana::LoadFrom<Spectrum> (dirInMC->GetDirectory("MCSig")).release();
119  Spectrum* fMCSigNuTrue = ana::LoadFrom<Spectrum> (nuTrueDir).release();
120  TH3* h3Deff = ana::ToTH3(*fMCSig , pot, kPOT, angbinsCustom, mukebins, eavailbins);
121  TH3* h3Deff_den = ana::ToTH3(*fMCSigNuTrue, pot, kPOT, angbinsCustom, mukebins, eavailbins);
122  TH3* h3DSel = (TH3*) h3Deff->Clone("h3DSel");
123  TH3* h3DSig = (TH3*) h3Deff_den->Clone("h3DSig");
124  TH2* h2Deff = (TH2*) h3Deff->Project3D("yx");
125  TH2* h2Deff_den = (TH2*) h3Deff_den->Project3D("yx");
126  h3Deff->Divide(h3Deff_den);
127  h2Deff->Divide(h2Deff_den);
128  h3Deff_den->Delete();
129  h2Deff_den->Delete();
130  delete fMCSig;
131  delete fMCSigNuTrue;
132 
133  //2. Signal:
134  Spectrum* fData = ana::LoadFrom<Spectrum> (dirInDt->GetDirectory(fake_xs.c_str())).release();
135  IBkgdEstimator* fBkgEst = ana::LoadFrom<IBkgdEstimator>(dirInMC->GetDirectory("BkgdEst")).release();
136  Spectrum *fBkg = new Spectrum(fBkgEst->Background());
137  Spectrum *fSignalEst = new Spectrum(*fData - *fBkg);
138  delete fBkgEst;
139  delete fBkg;
140 
141  //Background constraint
142  // Ratio fSignalEst_Bkg_Constraint = fData/inv_pur;
143  Ratio * fSignalEst_Bkg_Constraint = new Ratio(*fData, Spectrum(inv_pur->ToTH1(), fData->POT(), fData->Livetime()));
144 
145  //3. Unfolding:
146  ReweightableSpectrum* fRecoTrue = ana::LoadFrom<ReweightableSpectrum>(dirInMC->GetDirectory("RecoTrue")).release();
147 
148  UnfoldIterative * unfold = new UnfoldIterative(*fRecoTrue,5);
149  Spectrum * fSigUnfolded = new Spectrum(unfold->Truth(*fSignalEst));
150  delete unfold;
151 
152  UnfoldIterative * unfold_bkgc = new UnfoldIterative(*fRecoTrue,5);
153  Spectrum * fSigUnfolded_bkgc = new Spectrum(unfold_bkgc->Truth(Spectrum(fSignalEst_Bkg_Constraint->ToTH1(), fData->POT(), fData->Livetime())));
154  delete unfold_bkgc;
155  delete fRecoTrue;
156 
157  //4. Calculation:
158  double integrated_flux = hInFlux->Integral();
159  TH3* h3SigEst[2];
160  h3SigEst[0] = ana::ToTH3(*fSignalEst, pot, kPOT, angbinsCustom, mukebins, eavailbins);
161  h3SigEst[1] = ana::ToTH3(*fSignalEst_Bkg_Constraint, angbinsCustom, mukebins, eavailbins);
162  TH3* h3SigUnfolded[2];
163  h3SigUnfolded[0] = ana::ToTH3(*fSigUnfolded , pot, kPOT, angbinsCustom, mukebins, eavailbins);
164  h3SigUnfolded[1] = ana::ToTH3(*fSigUnfolded_bkgc, pot, kPOT, angbinsCustom, mukebins, eavailbins);
165  delete fSignalEst;
166  delete fSignalEst_Bkg_Constraint;
167  delete fSigUnfolded;
168  delete fSigUnfolded_bkgc;
169 
170  TH2* h2Dsigest[2];
171  TH1* h1Dsigest[2];
172  TH2* h2Dxs[2];
173  TH1* h1Dxs[2];
174  for(int ii=0;ii<2;ii++){
175  h3SigUnfolded[ii]->Divide(h3Deff);
176  h2Dxs[ii] = (TH2*)h3SigUnfolded[ii]->Project3D("yx");
177  h1Dxs[ii] = (TH1*)h3SigUnfolded[ii]->Project3D("z");
178 
179  h2Dxs[ii]->Scale(1./integrated_flux);
180  for(int ibin=1;ibin<=h1Dxs[ii]->GetNbinsX();ibin++){
181  double cont_num = h1Dxs[ii]->GetBinContent(ibin);
182  double cont_den = hInFlux->GetBinContent(ibin);
183  if(cont_den>0)h1Dxs[ii]->SetBinContent(ibin,cont_num/cont_den);
184  else h1Dxs[ii]->SetBinContent(ibin,0);
185  h1Dxs[ii]->SetBinError(ibin,0);
186  }
187  h2Dxs[ii]->Scale(1./ntargs);
188  h1Dxs[ii]->Scale(1./ntargs);
189  }
190 
191  //Storing
192  //xs_nominal: calcualting the signal estomation by S-B
193  fileOut->cd();
194  h3DSel->Write(("h3DSel_" + sdOut).c_str());
195  h3DSig->Write(("h3DSig_" + sdOut).c_str());
196  h3Deff->Write(("h3DEff_" + sdOut).c_str());
197  h2Deff->Write(("h2DEff_" + sdOut).c_str());
198  fileOut->cd(xs_mode0.c_str());
199  h2Dxs[0]->Write(("hXS2Dsec_" + sdOut).c_str());
200  h1Dxs[0]->Write(("hXS1Dsec_" + sdOut).c_str());
201  h3SigUnfolded[0]->Write(("hXS3Dsec_" + sdOut).c_str());
202  h3SigEst[0]->Write(("h3SigEst_" + sdOut).c_str());
203  //xs_bkgcons: calcualting the signal estomation by P*S
204  fileOut->cd();
205  fileOut->cd(xs_mode1.c_str());
206  h2Dxs[1]->Write(("hXS2Dsec_" + sdOut).c_str());
207  h1Dxs[1]->Write(("hXS1Dsec_" + sdOut).c_str());
208  h3SigUnfolded[1]->Write(("hXS3Dsec_" + sdOut).c_str());
209  h3SigEst[1]->Write(("h3SigEst_" + sdOut).c_str());
210 
211  //Leftover cleanup
212  h3DSel->Delete();
213  h3DSig->Delete();
214  h3Deff->Delete();
215  h2Deff->Delete();
216 }
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
Spectrum with the value of a second variable, allowing for reweighting
const std::string fake_xs
const Binning eavailbins
const double ntargs
virtual Spectrum Background() const =0
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const std::string xs_mode0
const std::string xs_mode1
std::vector< float > Spectrum
Definition: Constants.h:759
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 Binning mukebins
Definition: NumuCCIncBins.h:90
const double pot
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
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
void make_xs_1D ( TDirectory *&  dirInMC,
TDirectory *&  dirInDt,
std::string  sdOut,
TH1 *  hInFlux,
TFile *  fileOut,
Ratio inv_pur 
)

Getting all:

Definition at line 218 of file Make_NuMuCC_Inc_XS_v2.C.

References ana::IBkgdEstimator::Background(), om::cout, allTimeWatchdog::endl, fake_xs, ana::Spectrum::Livetime(), ana::LoadFrom< IBkgdEstimator >(), ntargs, pot, ana::Spectrum::POT(), runNovaSAM::release, ana::Spectrum::Scale(), ana::Ratio::ToTH1(), ana::Spectrum::ToTH1(), ana::UnfoldIterative::Truth(), xs_mode0, and xs_mode1.

Referenced by Make_NuMuCC_Inc_XS_v2().

218  {
219  std::cout<<"=> In make_xs_1D() "<< sdOut << std::endl;
220 
221  ///////////////////////////
222  ///Getting all:
223  //1. Efficiency:
224  TDirectory * nuTrueDir = dirInMC->GetDirectory("MCSigNuTrue");
225  if (!nuTrueDir || nuTrueDir == NULL || nuTrueDir == 0)
226  nuTrueDir = dirInMC->GetDirectory("MCSigNuTre"); // Try both titles to fix typo
227  Spectrum* fMCSig = ana::LoadFrom<Spectrum> (dirInMC->GetDirectory("MCSig")).release();
228  Spectrum* fMCSigNuTrue = ana::LoadFrom<Spectrum> (nuTrueDir).release();
229  TH1* h1Deff = fMCSig->ToTH1(pot);
230  TH1* h1Deff_den = fMCSigNuTrue->ToTH1(pot);
231  TH1* h1DSel = (TH1*) h1Deff->Clone("h1DSel");
232  TH1* h1DSig = (TH1*) h1Deff_den->Clone("h1DSig");
233  h1Deff->Divide(h1Deff_den);
234  h1Deff_den->Delete();
235  delete fMCSig;
236  delete fMCSigNuTrue;
237 
238  //2. Signal:
239  Spectrum* fData = ana::LoadFrom<Spectrum> (dirInDt->GetDirectory(fake_xs.c_str())).release();
240  IBkgdEstimator* fBkgEst = ana::LoadFrom<IBkgdEstimator>(dirInMC->GetDirectory("BkgdEst")).release();
241  Spectrum *fBkg = new Spectrum(fBkgEst->Background());
242  Spectrum *fSignalEst = new Spectrum(*fData - *fBkg);
243  delete fBkgEst;
244  delete fBkg;
245 
246  //Background constraint
247  // Ratio fSignalEst_Bkg_Constraint = fData/inv_pur;
248  Ratio * fSignalEst_Bkg_Constraint = new Ratio(*fData, Spectrum(inv_pur->ToTH1(), fData->POT(), fData->Livetime()));
249 
250  //3. Unfolding:
251  ReweightableSpectrum* fRecoTrue = ana::LoadFrom<ReweightableSpectrum>(dirInMC->GetDirectory("RecoTrue")).release();
252 
253  UnfoldIterative * unfold = new UnfoldIterative(*fRecoTrue,5);
254  Spectrum * fSigUnfolded = new Spectrum(unfold->Truth(*fSignalEst));
255  delete unfold;
256 
257  UnfoldIterative * unfold_bkgc = new UnfoldIterative(*fRecoTrue,5);
258  Spectrum * fSigUnfolded_bkgc = new Spectrum(unfold_bkgc->Truth(Spectrum(fSignalEst_Bkg_Constraint->ToTH1(), fData->POT(), fData->Livetime())));
259  delete unfold_bkgc;
260  delete fRecoTrue;
261 
262  //4. Calculation:
263  double integrated_flux = hInFlux->Integral();
264  TH1* h1SigEst[2];
265  h1SigEst[0] = fSignalEst->ToTH1(pot);
266  h1SigEst[1] = fSignalEst_Bkg_Constraint->ToTH1();
267  TH1* h1SigUnfolded[2];
268  h1SigUnfolded[0] = fSigUnfolded->ToTH1(pot);
269  h1SigUnfolded[1] = fSigUnfolded_bkgc->ToTH1(pot);
270  delete fSignalEst;
271  delete fSignalEst_Bkg_Constraint;
272  delete fSigUnfolded;
273  delete fSigUnfolded_bkgc;
274 
275  TH1* h1Dsigest[2];
276  TH1* h1Dxs[2];
277  for(int ii=0;ii<2;ii++){
278  h1SigUnfolded[ii]->Divide(h1Deff);
279  h1Dxs[ii] = h1SigUnfolded[ii];
280 
281  h1Dxs[ii]->Scale(1./integrated_flux);
282  // for(int ibin=1;ibin<=h1Dxs[ii]->GetNbinsX();ibin++){
283  // double cont_num = h1Dxs[ii]->GetBinContent(ibin);
284  // double cont_den = hInFlux->GetBinContent(ibin);
285  // if(cont_den>0)h1Dxs[ii]->SetBinContent(ibin,cont_num/cont_den);
286  // else h1Dxs[ii]->SetBinContent(ibin,0);
287  // h1Dxs[ii]->SetBinError(ibin,0);
288  // }
289  h1Dxs[ii]->Scale(1./ntargs);
290  }
291 
292  //Storing
293  //xs_nominal: calcualting the signal estomation by S-B
294  fileOut->cd();
295  h1DSel->Write(("h1DSel_" + sdOut).c_str());
296  h1DSig->Write(("h1DSig_" + sdOut).c_str());
297  h1Deff->Write(("h1DEff_" + sdOut).c_str());
298  fileOut->cd(xs_mode0.c_str());
299  h1Dxs[0]->Write(("hXS1Dsec_" + sdOut).c_str());
300  h1SigUnfolded[0]->Write(("hXS1Dsec_" + sdOut).c_str());
301  h1SigEst[0]->Write(("h1SigEst_" + sdOut).c_str());
302  //xs_bkgcons: calcualting the signal estomation by P*S
303  fileOut->cd();
304  fileOut->cd(xs_mode1.c_str());
305  h1Dxs[1]->Write(("hXS1Dsec_" + sdOut).c_str());
306  h1SigUnfolded[1]->Write(("hXS1Dsec_" + sdOut).c_str());
307  h1SigEst[1]->Write(("h1SigEst_" + sdOut).c_str());
308 
309  //Leftover cleanup
310  h1DSel->Delete();
311  h1DSig->Delete();
312  h1Deff->Delete();
313 }
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
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
Spectrum with the value of a second variable, allowing for reweighting
const std::string fake_xs
const double ntargs
virtual Spectrum Background() const =0
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const std::string xs_mode0
const std::string xs_mode1
std::vector< float > Spectrum
Definition: Constants.h:759
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 double pot
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230

Variable Documentation

const std::map<std::string, std::string> beamSysts
Initial value:
= {
}

Definition at line 91 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

std::map<std::string, std::string> coreFiles
Initial value:
{
{"nominal", "syst_specs_cv.root"},
{"fluxes", "fake_data_fluxes.root"},
{"fakedata", "fake_data_fluxes.root"},
{"genie", ""},
{"ppfx", ""}
}

Definition at line 53 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

std::string dirIn = "input_files"

Definition at line 48 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const std::string fake_xs = "ppfxwgt"

Definition at line 44 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by make_xs(), and make_xs_1D().

const double ntargs = 3.89137e+31

Definition at line 41 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by make_xs(), and make_xs_1D().

const int Nuniv_genie = 0

Definition at line 51 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const int Nuniv_ppfx = 0

Definition at line 50 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const double pot = 8.09e20

Definition at line 42 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2(), make_xs(), and make_xs_1D().

const std::map<std::string, std::string> systFiles
Initial value:
= {
{"ckv", "syst_specs_ckv.root"},
{"angleup", "syst_specs_angleup.root"},
{"angledw", "syst_specs_angledw.root"},
{"neutron_Up", "syst_specs_neutron.root"},
{"neutron_Dw", "syst_specs_neutron.root"},
{"neutron_nom", "syst_specs_neutron_nom.root"}
}

Definition at line 76 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const std::vector<std::string> unfoldVars1D
Initial value:
{
"ENu",
"Q2"
}

Definition at line 65 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const std::vector<std::string> unfoldVars3D
Initial value:
{
"EAvail"
}

Definition at line 61 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const std::vector<std::string> updownlabels = {"_Up", "_Dw"}

Definition at line 105 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

std::map<std::string, std::string> var2FD
Initial value:
{
{"EAvail", "fakedata"},
{"ENu", "fakedata_enu"},
{"Q2", "fakedata_q2"}
}

Definition at line 70 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2().

const std::string xs_mode0 = "xs_nominal"

Definition at line 45 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2(), make_xs(), and make_xs_1D().

const std::string xs_mode1 = "xs_bkgcons"

Definition at line 46 of file Make_NuMuCC_Inc_XS_v2.C.

Referenced by Make_NuMuCC_Inc_XS_v2(), make_xs(), and make_xs_1D().