nue_decomp_scales.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void nue_decomp_scales(bool makeFile = false)
3 {
4  std::cout << "Sorry, you must run in compiled mode" << std::endl;
5 }
6 #else
7 
8 #include "CAFAna/Core/HistAxis.h"
9 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Cuts/TimingCuts.h"
12 #include "CAFAna/Cuts/Cuts.h"
13 #include "CAFAna/Cuts/NueCutsSecondAna.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Decomp/IDecomp.h"
26 #include "CAFAna/Vars/HistAxes.h"
30 #include "CAFAna/Vars/Vars.h"
31 #include "CAFAna/Vars/XsecTunes.h"
32 
34 
35 #include "TCanvas.h"
36 #include "TFile.h"
37 #include "TH1.h"
38 #include "TLegend.h"
39 #include "TPad.h"
40 
41 #include <iostream>
42 #include <memory>
43 
44 using namespace ana;
45 
46 
47 void nue_decomp_scales(bool makeFile = false)
48 {
49 
50  const HistAxis axisNue = {"Nue Energy / Selection Bin", kNue2020Binning, kNue2020AnaBin};
51  std::string filename = "nue2020_decomp_spectra_calibdown.root";
52 
53  if(makeFile){
54 
55 
56  SpectrumLoader loader_data("def_snapshot prod_caf_R19-11-18-prod5reco.d.f.h_nd_numi_fhc_period235910_v1 with limit 200");
57  SpectrumLoader loader_mc("def_snapshot prod_caf_R19-11-18-prod5reco.d.h_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1 with limit 200");
58 
60  loader_data.SetSpillCut(kStandardSpillCuts);
61 
62  BENDecomp* decompBEN = new BENDecomp(loader_mc, loader_data, axisNue, kNue2020ND, kNoShift, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
63  MichelDecomp* decompMichel = new MichelDecomp(loader_mc, loader_data, axisNue, kNue2020ND, decompBEN, kNoShift, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
64  ProportionalDecomp* decompProportional = new ProportionalDecomp(loader_mc, loader_data, axisNue, kNue2020ND, kNoShift, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2020);
65 
66  Spectrum* data = new Spectrum(loader_data, axisNue, kNue2020ND, kNoShift);
67  PredictionNoExtrap* mc_nom = new PredictionNoExtrap(loader_mc, kNullLoader, kNullLoader, axisNue.GetLabels()[0], axisNue.GetBinnings()[0], axisNue.GetVars()[0],
69 
70  loader_mc.Go();
71  loader_data.Go();
72  TFile* fout = new TFile(filename.c_str(), "recreate");
73 
74  data->SaveTo(fout, "data");
75  mc_nom->SaveTo(fout, "mc_nom");
76  decompBEN->SaveTo(fout, "decompBEN");
77  decompMichel->SaveTo(fout, "decompMichel");
78  decompProportional->SaveTo(fout, "decompProportional");
79 
80  delete fout;
81  }
82  else {
83 
84  TFile* fin = new TFile(filename.c_str(), "read");
85  std::cout << "Extracting scales...";
86 
87  Spectrum* spec_data = Spectrum::LoadFrom(fin, "data").release();
88  double pot = spec_data->POT();
89  std::cout << pot << std::endl;
90  PredictionNoExtrap* mc_nom = PredictionNoExtrap::LoadFrom(fin, "mc_nom").release();
91  MichelDecomp* michel = MichelDecomp::LoadFrom(fin, "decompMichel").release();
92 
93  TH1* data = spec_data->ToTH1(pot);
94  TH1* nue_nom = mc_nom->ComponentCC(12, 12).Unoscillated().ToTH1(pot);
95  TH1* numu_nom = mc_nom->ComponentCC(14, 14).Unoscillated().ToTH1(pot);
96  TH1* antinue_nom = mc_nom->ComponentCC(-12, -12).Unoscillated().ToTH1(pot);
97  TH1* antinumu_nom = mc_nom->ComponentCC(-14, -14).Unoscillated().ToTH1(pot);
98  TH1* nc_nom = mc_nom->ComponentNC().ToTH1(pot);
99 
100  // Michel
101 
102  TH1* nue_ratio = (TH1*)((michel->NueComponent().ToTH1(pot))->Clone("combo_beam_nue"));
103  nue_ratio->Divide(nue_nom);
104 
105  TH1* antinue_ratio = (TH1*)((michel->AntiNueComponent().ToTH1(pot))->Clone("combo_beam_antinue"));
106  antinue_ratio->Divide(antinue_nom);
107 
108  TH1* antinumu_ratio = (TH1*)((michel->AntiNumuComponent().ToTH1(pot))->Clone("combo_beam_antinumu"));
109  antinumu_ratio->Divide(antinumu_nom);
110 
111  TH1* numu_ratio = (TH1*)((michel->NumuComponent().ToTH1(pot))->Clone("combo_beam_numu"));
112  numu_ratio->Divide(numu_nom);
113 
114  TH1* nc_ratio = (TH1*)((michel->NCComponent().ToTH1(pot))->Clone("combo_beam_nc"));
115  nc_ratio->Divide(nc_nom);
116 
117  TFile* outputFile = new TFile("nue2017_decomp_ratios.root", "update");
118  outputFile->cd();
119  nue_ratio->Write();
120  numu_ratio->Write();
121  nc_ratio->Write();
122  antinue_ratio->Write();
123  antinumu_ratio->Write();
124  outputFile->Close();
125  delete outputFile;
126 
127  // TCanvas *c0 = new TCanvas("c0", "c0", 200, 10, 800, 600);
128  // c0->cd();
129  // nue_ratio->GetYaxis()->SetRangeUser(0., 2.);
130  // nue_ratio->SetLineColor(kPink);
131  // nue_ratio->Draw("hist");
132  // numu_ratio->SetLineColor(kGreen+2);
133  // numu_ratio->Draw("hist same");
134  // nc_ratio->SetLineColor(kBlue);
135  // nc_ratio->Draw("hist same");
136  // c0->Update();
137  // c0->Print("michel_ratios.pdf");
138 
139 
140 
141 
142  }
143 
144 
145 }
146 #endif
OscillatableSpectrum ComponentCC(int from, int to) const override
TString fin
Definition: Style.C:24
Spectrum ComponentNC() const override
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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
virtual Spectrum AntiNumuComponent() const override
virtual Spectrum NCComponent() const override
string filename
Definition: shutoffs.py:106
void SetSpillCut(const SpillCut &cut)
virtual Spectrum NumuComponent() const override
const Cut kNue2020ND
Definition: NueCuts2020.h:178
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
virtual Spectrum AntiNueComponent() const override
static std::unique_ptr< MichelDecomp > LoadFrom(TDirectory *dir, const std::string &name)
const XML_Char const XML_Char * data
Definition: expat.h:268
const Binning kNue2020Binning
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
virtual Spectrum NueComponent() const override
#define pot
void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: BENDecomp.cxx:387
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
std::vector< float > Spectrum
Definition: Constants.h:610
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
double POT() const
Definition: Spectrum.h:227
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
void SaveTo(TDirectory *dir, const std::string &name) const override
const HistAxis axisNue
Definition: syst_header.h:378
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
Splits Data proportionally according to MC.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kNue2020AnaBin([](const caf::SRProxy *sr){int selBin=kNue2020SelectionBin(sr);float nuE=kNueEnergy2020(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2020, official Binning.
Definition: NueCuts2020.h:191
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Prediction that just uses FD MC, with no extrapolation.
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
void nue_decomp_scales(bool makeFile=false)
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106
enum BeamMode string