calculateComponentsNumu.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////
2 // //
3 // WRONG SIGN NORMALIZATION STUDY SELECTION //
4 // NUMU ONLY! //
5 // //
6 /////////////////////////////////////////////////////////////
7 
10 #include "CAFAna/Fit/Fit.h"
12 #include "CAFAna/Analysis/Style.h"
13 #include "CAFAna/Vars/FitVars.h"
14 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/Utilities.h"
17 #include "CAFAna/Analysis/Calcs.h"
24 #include "CAFAna/Core/Progress.h"
31 #include "OscLib/OscCalcPMNSOpt.h"
32 #include "OscLib/OscCalcDumb.h"
33 #include "Utilities/func/MathUtil.h"
35 
36 using namespace ana;
37 
38 #include "TFile.h"
39 #include "TH2.h"
40 #include "TMath.h"
41 #include "TCanvas.h"
42 #include "TGraph.h"
43 #include "TLatex.h"
44 #include "TLegend.h"
45 #include "TBox.h"
46 #include "TLine.h"
47 #include "TSystem.h"
48 #include "TColor.h"
49 #include "TPaletteAxis.h"
50 #include "THStack.h"
51 #include "TString.h"
52 #include "TEllipse.h"
53 
54 #include <iostream>
55 #include <iomanip>
56 #include <string>
57 #include <vector>
58 #include "Utilities/rootlogon.C"
59 
61 {
62  TFile *inFileMC = new TFile("/nova/ana/users/howard/Ana2019-WrongSign/WrongSign-XCheck2019-NumuMC.root","read");
63 
65 
66  // ----- Load files -----
67  std::vector<std::string> typeString = {"RS","WS"};
68 
69  PredictionNoOsc *pRHCNumu[2]; // FSProton RS Enhanced, FSProton WS Enhanced
70  Spectrum *sNCNumu[2]; // FSProton RS Enhanced, FSProton WS Enhanced
71  Spectrum *sNCAntiNumu[2]; // FSProton RS Enhanced, FSProton WS Enhanced
72 
73  PredictionNoOsc *pRHCNumuTot;
74  Spectrum *sNCNumuTot;
75  Spectrum *sNCAntiNumuTot;
76 
77  for( int iMd=0; iMd<2; ++iMd ){
78  pRHCNumu[iMd] = PredictionNoOsc::LoadFrom(inFileMC, TString::Format("pNumuFSProton%s",
79  typeString[iMd].c_str()).Data() ).release();
80  sNCNumu[iMd] = Spectrum::LoadFrom(inFileMC, TString::Format("sNumuFSProton%sNC",
81  typeString[iMd].c_str()).Data() ).release();
82  sNCAntiNumu[iMd] = Spectrum::LoadFrom(inFileMC, TString::Format("sNumuFSProton%sNCAnti",
83  typeString[iMd].c_str()).Data() ).release();
84  }
85 
86  pRHCNumuTot = PredictionNoOsc::LoadFrom(inFileMC, "pNumuTot").release();
87  sNCNumuTot = Spectrum::LoadFrom(inFileMC, "sNumuTotNC").release();
88  sNCAntiNumuTot = Spectrum::LoadFrom(inFileMC, "sNumuTotNCAnti").release();
89 
90  delete inFileMC;
91 
92  // --- USE POT FROM PERIOD DATA ---
93  //TFile *inFileData = new TFile( "/nova/ana/users/howard/Ana2019-WrongSign/WrongSign-XCheck2019-NumuDataPeriods4-6.root", "read" );
94  TFile *inFileData = new TFile( "/nova/ana/users/howard/Ana2019-WrongSign/WrongSign-XCheck2019-NumuFullDataPeriods.root", "read" );
95  Spectrum *dataSpecRS = Spectrum::LoadFrom(inFileData, "sQ1NumuData_FSProtonRS").release();
96  Spectrum *dataSpecWS = Spectrum::LoadFrom(inFileData, "sQ1NumuData_FSProtonWS").release();
97  Spectrum *dataSpecTot = Spectrum::LoadFrom(inFileData, "sQ1NumuDataTot").release();
98 
99  auto pot = dataSpecRS->POT();
100  delete inFileData;
101 
102  std::cout << "----------------------------" << std::endl;
103  std::cout << " POT: " << pot << std::endl;
104  std::cout << "----------------------------" << std::endl;
105  ///////////////////////////////////////////////////////////////////////////////////
106 
107  // ----- Numu -----
108  // The 2 is b/c RS Enhanced, WS Enhanced
109  TH1D *numuS[3][2]; // make this also have 3 categories, so I can handle NC as part of signal (and NCAnti alone as part of BG)
110  TH1D *numuB[3][2];
111 
112  TH1D *numuTotS[3];
113  TH1D *numuTotB[3];
114 
115  TH1D *numuDataH[2];
116  TH1D *numuDataTotH;
117 
118  for( int iMd=0; iMd<2; ++iMd ){
119  // Get predictions
120  numuS[0][iMd] = pRHCNumu[iMd]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
121  numuS[0][iMd]->SetName( TString::Format("NueFSProton_%s",typeString[iMd].c_str()) );
122  numuS[1][iMd] = pRHCNumu[iMd]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
123  numuS[1][iMd]->SetName( TString::Format("NumuFSProton_%s",typeString[iMd].c_str()) );
124  numuB[0][iMd] = pRHCNumu[iMd]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
125  numuB[0][iMd]->SetName( TString::Format("NuebarFSProton_%s",typeString[iMd].c_str()) );
126  numuB[1][iMd] = pRHCNumu[iMd]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
127  numuB[1][iMd]->SetName( TString::Format("NumubarFSProton_%s",typeString[iMd].c_str()) );
128  // since split NC, handle differently
129  numuB[2][iMd] = sNCAntiNumu[iMd]->ToTH1(pot); // NCAnti
130  numuB[2][iMd]->SetName( TString::Format("NcAntiFSProton_%s",typeString[iMd].c_str()) );
131  numuS[2][iMd] = sNCNumu[iMd]->ToTH1(pot); // NC
132  numuS[2][iMd]->SetName( TString::Format("NcFSProton_%s",typeString[iMd].c_str()) );
133  }
134 
135  // Handle Data
136  numuDataH[0] = dataSpecRS->ToTH1(pot);
137  numuDataH[0]->SetTitle("numuData_RS");
138  numuDataH[1] = dataSpecWS->ToTH1(pot);
139  numuDataH[1]->SetTitle("numuData_WS");
140 
141  numuDataTotH = dataSpecTot->ToTH1(pot);
142  numuDataTotH->SetTitle("numuDataTot");
143 
144  // Total ones
145  numuTotS[0] = pRHCNumuTot->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
146  numuTotS[0]->SetName( "NueFSProtonTot" );
147  numuTotS[1] = pRHCNumuTot->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
148  numuTotS[1]->SetName( "NumuFSProtonTot" );
149  numuTotB[0] = pRHCNumuTot->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
150  numuTotB[0]->SetName( "NuebarFSProtonTot" );
151  numuTotB[1] = pRHCNumuTot->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
152  numuTotB[1]->SetName( "NumubarFSProtonTot" );
153  numuTotB[2] = sNCAntiNumuTot->ToTH1(pot); // NCAnti
154  numuTotB[2]->SetName( "NcAntiFSProtonTot" );
155  numuTotS[2] = sNCNumuTot->ToTH1(pot); // NC
156  numuTotS[2]->SetName( "NcFSProtonTot" );
157 
158  // PRINT OUT SUMS
159  std::cout << "___________MC SELECTIONS___________" << std::endl;
160  std::cout << std::endl;
161  std::cout << "Using FSProton splits:" << std::endl;
162  std::cout << " numu(bar) = " << numuB[1][0]->Integral() + numuB[1][1]->Integral() + numuS[1][0]->Integral() + numuS[1][1]->Integral() << std::endl;
163  std::cout << " nue(bar) = " << numuB[0][0]->Integral() + numuB[0][1]->Integral() + numuS[0][0]->Integral() + numuS[0][1]->Integral() << std::endl;
164  std::cout << " NC = " << numuB[2][0]->Integral() + numuB[2][1]->Integral() + numuS[2][0]->Integral() + numuS[2][1]->Integral() << std::endl;
165  std::cout << "------------------" << std::endl;
166  std::cout << " WSCC = " << numuS[1][0]->Integral() + numuS[1][1]->Integral() + numuS[0][0]->Integral() + numuS[0][1]->Integral() << std::endl;
167  std::cout << " WSNC = " << numuS[2][0]->Integral() + numuS[2][1]->Integral() << std::endl;
168  std::cout << " RSCC = " << numuB[1][0]->Integral() + numuB[1][1]->Integral() + numuB[0][0]->Integral() + numuB[0][1]->Integral() << std::endl;
169  std::cout << " RSNC = " << numuB[2][0]->Integral() + numuB[2][1]->Integral() << std::endl;
170  std::cout << " Tot Data = " << numuDataH[0]->Integral() + numuDataH[1]->Integral() << std::endl;
171 
172  std::cout << std::endl;
173  std::cout << std::endl;
174 
175  std::cout << "Using totals:"<< std::endl;
176  std::cout << " numu(bar) = " << numuTotB[1]->Integral() + numuTotS[1]->Integral() << std::endl;
177  std::cout << " nue(bar) = " << numuTotB[0]->Integral() + numuTotS[0]->Integral() << std::endl;
178  std::cout << " NC = " << numuTotB[2]->Integral() + numuTotS[2]->Integral() << std::endl;
179  std::cout << "------------------"<< std::endl;
180  std::cout << " WSCC = " << numuTotS[0]->Integral() + numuTotS[1]->Integral() << std::endl;
181  std::cout << " WSNC = " << numuTotS[2]->Integral() << std::endl;
182  std::cout << " RSCC = " << numuTotB[0]->Integral() + numuTotB[1]->Integral() << std::endl;
183  std::cout << " RSNC = " << numuTotB[2]->Integral() << std::endl;
184  std::cout << " Tot Data = " << numuDataTotH->Integral() << std::endl;
185 
186  std::cout << std::endl;
187  std::cout << std::endl;
188  std::cout << "Done." << std::endl;
189 }
void calculateComponentsNumu()
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
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
static std::unique_ptr< PredictionNoOsc > LoadFrom(TDirectory *dir, const std::string &name)
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Charged-current interactions.
Definition: IPrediction.h:39
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
#define pot
double POT() const
Definition: Spectrum.h:227
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Prediction that wraps a simple Spectrum.