CalcCutValsPlot.C
Go to the documentation of this file.
1 //......................................................
5 
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/Spectrum.h"
10 
13 
15 
16 #include "Utilities/rootlogon.C"
17 
18 #include "OscLib/OscCalcPMNSOpt.h"
19 
20 #include "TCanvas.h"
21 #include "TColor.h"
22 #include "TFile.h"
23 #include "TGraph.h"
24 #include "TGraphAsymmErrors.h"
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "TLatex.h"
28 #include "TLegend.h"
29 #include "TStyle.h"
30 #include "TSystem.h"
31 
32 #include <fstream>
33 #include <iostream>
34 #include <iomanip>
35 #include <string>
36 
37 bool isFHC = false; // Global variable for whether I am running over FHC or RHC?
38 std::string sFHC = "RHC"; // Global variable for string related to FHC or RHC.
39 
40 //......................................................
41 using namespace ana;
42 
43 //......................................................
44 TH1D* GetHistVectors( TFile* InF, std::string InDir );
45 //......................................................
46 void CutThreshold( TH1D* MyHist, std::string Key );
47 //......................................................
48 void MakeCanvas( TH1D* Period1, TH1D* Period2, TH1D* FHCHigh, TH1D* RHCHigh, bool isBPF, std::string Ident );
49 //......................................................
50 
51 // ---- Define my output directory
53 
54 double POTNorm = 0;
55 double LivNorm = 0;
56 
57 double CutAim = 0;
58 
59 //......................................................
60 void CalcCutValsPlot( std::string sIdentifier, bool fhc ) {
61  // Set my global variables.
62  isFHC = fhc;
63  sFHC = (isFHC == true ? "fhc":"rhc" );
64 
65  if ( fhc && sIdentifier == "Period1" ) {
68  CutAim = 7.59;
69  }
70  else if ( fhc && sIdentifier == "Period2" ) {
73  CutAim = 37.88;
74  }
75  else if ( fhc && sIdentifier == "HighGain" ) {
78  CutAim = 80.41;
79  }
80  else if (!fhc && sIdentifier == "HighGain" ) {
83  CutAim = 53.54;
84  }
85  else {
86  std::cout << "\n\n I don't recognise sIdentifier = " << sIdentifier << " quitting.\n\n" << std::endl;
87  return;
88  }
89 
90  gStyle->SetHistMinimumZero();
91 
92  // Where are my input files?
93  //std::string InBase = "/nova/app/users/karlwarb/Workspace/Post-Ana18/CosRej/CalcCutVals/";
94  //std::string MyFile = "CalcCutVals-"+sFHC+"-"+sIdentifier+".root";
95  std::string InBase = "/nova/ana/users/karlwarb/PostAna18/CosRej/CalcCutVals-190404/";
96  std::string MyFile = "CalcVals-"+sFHC+"-"+sIdentifier+".root";
97 
98  std::string FileLoc = InBase + MyFile;
99  std::cerr << "Loading MyFile - " << FileLoc << std::endl;
100  TFile *ThisFile = TFile::Open( FileLoc.c_str() );
101 
102  // --- Get my output file sorted out...
103  OutDir = "Plots-"+sFHC+"-"+sIdentifier+"/";
104  gSystem -> MakeDirectory( OutDir.c_str() );
105  std::string OutName = OutDir+"CalcCutValsPlots_"+sFHC+"-"+sIdentifier+".root";
106  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
107 
108  // -----------------------------------------------------------
109  // --- Define my histograms.
110  // -----------------------------------------------------------
111  // ================
112  // = Data Spectra =
113  // ================
114  // --- BPF spectra.
115  /*
116  TH1D* hData_FHC_BPF_MLP_Per1 = GetHistVectors( ThisFile, "Data_FHC_BPF_MLP_Per1" );
117  TH1D* hData_FHC_BPF_BDT_Per1 = GetHistVectors( ThisFile, "Data_FHC_BPF_BDT_Per1" );
118  TH1D* hData_FHC_BPF_MLP_Per2 = GetHistVectors( ThisFile, "Data_FHC_BPF_MLP_Per2" );
119  TH1D* hData_FHC_BPF_BDT_Per2 = GetHistVectors( ThisFile, "Data_FHC_BPF_BDT_Per2" );
120  TH1D* hData_FHC_BPF_MLP_High = GetHistVectors( ThisFile, "Data_FHC_BPF_MLP_High" );
121  TH1D* hData_FHC_BPF_BDT_High = GetHistVectors( ThisFile, "Data_FHC_BPF_BDT_High" );
122  TH1D* hData_RHC_BPF_MLP_High = GetHistVectors( ThisFile, "Data_RHC_BPF_MLP_High" );
123  TH1D* hData_RHC_BPF_BDT_High = GetHistVectors( ThisFile, "Data_RHC_BPF_BDT_High" );
124 
125  // --- Kalman spectra.
126  TH1D* hData_FHC_Kal_MLP_Per1 = GetHistVectors( ThisFile, "Data_FHC_Kal_MLP_Per1" );
127  TH1D* hData_FHC_Kal_BDT_Per1 = GetHistVectors( ThisFile, "Data_FHC_Kal_BDT_Per1" );
128  TH1D* hData_FHC_Kal_MLP_Per2 = GetHistVectors( ThisFile, "Data_FHC_Kal_MLP_Per2" );
129  TH1D* hData_FHC_Kal_BDT_Per2 = GetHistVectors( ThisFile, "Data_FHC_Kal_BDT_Per2" );
130  TH1D* hData_FHC_Kal_MLP_High = GetHistVectors( ThisFile, "Data_FHC_Kal_MLP_High" );
131  TH1D* hData_FHC_Kal_BDT_High = GetHistVectors( ThisFile, "Data_FHC_Kal_BDT_High" );
132  TH1D* hData_RHC_Kal_MLP_High = GetHistVectors( ThisFile, "Data_RHC_Kal_MLP_High" );
133  TH1D* hData_RHC_Kal_BDT_High = GetHistVectors( ThisFile, "Data_RHC_Kal_BDT_High" );
134  */
135  // ================
136  // = Mont Spectra =
137  // ================
138  // --- BPF spectra.
139  TH1D* hMont_FHC_BPF_MLP_Per1 = GetHistVectors( ThisFile, "Mont_FHC_BPF_MLP_Per1" );
140  TH1D* hMont_FHC_BPF_BDT_Per1 = GetHistVectors( ThisFile, "Mont_FHC_BPF_BDT_Per1" );
141  TH1D* hMont_FHC_BPF_MLP_Per2 = GetHistVectors( ThisFile, "Mont_FHC_BPF_MLP_Per2" );
142  TH1D* hMont_FHC_BPF_BDT_Per2 = GetHistVectors( ThisFile, "Mont_FHC_BPF_BDT_Per2" );
143  TH1D* hMont_FHC_BPF_MLP_High = GetHistVectors( ThisFile, "Mont_FHC_BPF_MLP_High" );
144  TH1D* hMont_FHC_BPF_BDT_High = GetHistVectors( ThisFile, "Mont_FHC_BPF_BDT_High" );
145  TH1D* hMont_RHC_BPF_MLP_High = GetHistVectors( ThisFile, "Mont_RHC_BPF_MLP_High" );
146  TH1D* hMont_RHC_BPF_BDT_High = GetHistVectors( ThisFile, "Mont_RHC_BPF_BDT_High" );
147 
148  // --- Kalman spectra.
149  TH1D* hMont_FHC_Kal_MLP_Per1 = GetHistVectors( ThisFile, "Mont_FHC_Kal_MLP_Per1" );
150  TH1D* hMont_FHC_Kal_BDT_Per1 = GetHistVectors( ThisFile, "Mont_FHC_Kal_BDT_Per1" );
151  TH1D* hMont_FHC_Kal_MLP_Per2 = GetHistVectors( ThisFile, "Mont_FHC_Kal_MLP_Per2" );
152  TH1D* hMont_FHC_Kal_BDT_Per2 = GetHistVectors( ThisFile, "Mont_FHC_Kal_BDT_Per2" );
153  TH1D* hMont_FHC_Kal_MLP_High = GetHistVectors( ThisFile, "Mont_FHC_Kal_MLP_High" );
154  TH1D* hMont_FHC_Kal_BDT_High = GetHistVectors( ThisFile, "Mont_FHC_Kal_BDT_High" );
155  TH1D* hMont_RHC_Kal_MLP_High = GetHistVectors( ThisFile, "Mont_RHC_Kal_MLP_High" );
156  TH1D* hMont_RHC_Kal_BDT_High = GetHistVectors( ThisFile, "Mont_RHC_Kal_BDT_High" );
157 
158 
159  std::cout << "\n\nOut of there..."
160  << "\n\t hMont_FHC_BPF_MLP_Per1 -> " << hMont_FHC_BPF_MLP_Per1->Integral()
161  << "\n\t hMont_FHC_Kal_MLP_Per1 -> " << hMont_FHC_Kal_MLP_Per1->Integral()
162  << "\n\t hMont_FHC_BPF_BDT_High -> " << hMont_FHC_BPF_BDT_High->Integral()
163  << "\n\t hMont_RHC_Kal_MLP_High -> " << hMont_RHC_Kal_MLP_High->Integral()
164  << std::endl;
165 
166  // --- Now lets start drawing things!!!
167  std::cerr << "\n\nNow lets start plotting things..." << std::endl;
168  OutFile -> cd();
169 
170  // --- BPF --- MLP
171  MakeCanvas( hMont_FHC_BPF_MLP_Per1, hMont_FHC_BPF_MLP_Per2, hMont_FHC_BPF_MLP_High, hMont_RHC_BPF_MLP_High, true, sFHC+"-"+sIdentifier+"-BPF-MLP" );
172 
173  // --- BPF --- BDT
174  MakeCanvas( hMont_FHC_BPF_BDT_Per1, hMont_FHC_BPF_BDT_Per2, hMont_FHC_BPF_BDT_High, hMont_RHC_BPF_BDT_High, true, sFHC+"-"+sIdentifier+"-BPF-BDT" );
175 
176  // --- Kal --- MLP
177  MakeCanvas( hMont_FHC_Kal_MLP_Per1, hMont_FHC_Kal_MLP_Per2, hMont_FHC_Kal_MLP_High, hMont_RHC_Kal_MLP_High, true, sFHC+"-"+sIdentifier+"-Kal-MLP" );
178 
179  // --- Kal --- BDT
180  MakeCanvas( hMont_FHC_Kal_BDT_Per1, hMont_FHC_Kal_BDT_Per2, hMont_FHC_Kal_BDT_High, hMont_RHC_Kal_BDT_High, true, sFHC+"-"+sIdentifier+"-Kal-BDT" );
181 
182 
183  return;
184 }
185 
186 //......................................................
187 TH1D* GetHistVectors( TFile* InF, std::string InDir ) {
188  InF -> cd();
189  TH1D* TempHist;
190 
191  if ( InDir.find("Data") != std::string::npos ) {
192 
193  // --- If Cosmics then get a Spectrum and normalise by livetime
194  std::unique_ptr<Spectrum> Temp = Spectrum::LoadFrom( InF -> GetDirectory( InDir.c_str() ) );
195  TempHist = Temp -> ToTH1( LivNorm, kLivetime );
196  }
197  else if ( InDir.find("Mont") != std::string::npos ) {
198 
199  // --- If Monte Carlo then get a PredictionNoExtrap and normalise by POT.
200  std::unique_ptr<PredictionNoExtrap> TempSpec = PredictionNoExtrap::LoadFrom( InF -> GetDirectory( InDir.c_str() ) );
201  // --- Need to set my oscillation calculator...
204  calc.SetdCP ( M_PI*0.915869 ); // preliminary Ana2017
205  calc.SetDmsq32( 0.00248809 );
206  calc.SetTh23 ( asin(sqrt(0.587685)) );
207  // --- Now I can convert to a histogram
208  TempHist = TempSpec -> Predict(&calc).ToTH1( POTNorm );
209  }
210  else {
211 
212  // --- If not, then I passed the wrong argument...
213  std::cout << "\nYou passed " << InDir << " which doesn't match anything!" << std::endl;
214  return TempHist;
215  }
216 
217  TempHist->GetXaxis()->SetTitle( "Cosmic Rej Alg Score" );
218  // --- And return!
219  return TempHist;
220 }
221 
222 //......................................................
223 void CutThreshold( TH1D* MyHist, std::string Key ) {
224 
225  //for (int bin=0; bin < MyHist->GetNbinsX()+1; ++bin) {
226  //for (int bin=MyHist->GetNbinsX()+1; bin>=0; --bin) {
227  for (int bin = MyHist->GetNbinsX(); bin >= 1; --bin) {
228 
229  double NewBinCent = MyHist->GetBinLowEdge(bin);
230  double NewInt = MyHist->Integral( bin, MyHist->GetNbinsX() );
231  double NewIntFrac = MyHist->Integral( bin, MyHist->GetNbinsX() ) / MyHist->Integral();
232 
233  if ( NewInt > CutAim ) {
234 
235  double LooBinCent = MyHist->GetBinLowEdge(bin-1);
236  double LooIntFrac = MyHist->Integral( bin-1, MyHist->GetNbinsX() ) / MyHist->Integral();
237  double LooInt = MyHist->Integral( bin-1, MyHist->GetNbinsX() );
238 
239  double StrBinCent = MyHist->GetBinLowEdge(bin+1);
240  double StrIntFrac = MyHist->Integral( bin+1, MyHist->GetNbinsX() ) / MyHist->Integral();
241  double StrInt = MyHist->Integral( bin+1, MyHist->GetNbinsX() );
242 
243  std::cout << std::setprecision(3) << std::fixed
244  << "\t Found my threshold for when looking at the " << Key << " Var..."
245  << "\n\t\t Loose Cut bin " << LooBinCent << " --> Int " << LooInt << ", IntFrac " << LooIntFrac
246  << "\n\t\t Match Cut bin " << NewBinCent << " --> Int " << NewInt << ", IntFrac " << NewIntFrac
247  << "\n\t\t Strict Cut bin " << StrBinCent << " --> Int " << StrInt << ", IntFrac " << StrIntFrac
248  << std::endl;
249 
250  break;
251  }
252  }
253 }
254 
255 //......................................................
256 void MakeCanvas( TH1D* Period1, TH1D* Period2, TH1D* FHCHigh, TH1D* RHCHigh, bool isBPF, std::string Ident ) {
257 
258  std::cout << "\n\nNow looking at " << Ident
259  /*
260  << "\n\t Period 1 Integral is: " << Period1 -> Integral()
261  << "\n\t Period 2 Integral is: " << Period2 -> Integral()
262  << "\n\t FHC High Integral is: " << FHCHigh -> Integral()
263  << "\n\t RHC High Integral is: " << RHCHigh -> Integral()
264  */
265  << "\n" << std::endl;
266 
267  CutThreshold( Period1, "Period 1" );
268  CutThreshold( Period2, "Period 2" );
269  CutThreshold( FHCHigh, "FHC High" );
270  CutThreshold( RHCHigh, "RHC High" );
271 
272  // ---- Let make a dummy canvas...
273  TCanvas *HistCan = new TCanvas( Ident.c_str(), Ident.c_str() );
274  if (Ident.find("MLP") != std::string::npos) {
275  HistCan -> SetLogy();
276  Period1 -> GetXaxis() -> SetRangeUser( 0.95, 1.0 ); Period1 -> GetYaxis() -> SetRangeUser( 1e-3, 1e2 );
277  } else {
278  Period1 -> GetXaxis() -> SetRangeUser( 0.4, 0.8 );
279 
280  Period1 -> Rebin( 4 ); Period1 -> GetXaxis() -> SetRangeUser( 0.4, 0.8 );
281  Period2 -> Rebin( 4 ); Period2 -> GetXaxis() -> SetRangeUser( 0.4, 0.8 );
282  FHCHigh -> Rebin( 4 ); FHCHigh -> GetXaxis() -> SetRangeUser( 0.4, 0.8 );
283  RHCHigh -> Rebin( 4 ); RHCHigh -> GetXaxis() -> SetRangeUser( 0.4, 0.8 );
284  }
285 
286  Period1 -> SetLineColor( 1 );
287  Period2 -> SetLineColor( 2 );
288  FHCHigh -> SetLineColor( 4 );
289  RHCHigh -> SetLineColor( 6 );
290 
291  Period1 -> Draw();
292  Period2 -> Draw("same");
293  FHCHigh -> Draw("same");
294  RHCHigh -> Draw("same");
295 
296  TLegend *l = new TLegend(0.15,0.65,0.4,0.88);
297  l -> AddEntry( Period1, "Period 1", "l" );
298  l -> AddEntry( Period2, "Period 2", "l" );
299  l -> AddEntry( FHCHigh, "FHC High", "l" );
300  l -> AddEntry( RHCHigh, "RHC High", "l" );
301  l -> Draw();
302 
303  HistCan -> Write ( Ident.c_str() );
304  HistCan -> SaveAs( TString(OutDir)+TString(Ident)+TString(".pdf" ) );
305  HistCan -> SaveAs( TString(OutDir)+TString(Ident)+TString(".png" ) );
306 }
307 //......................................................
tree Draw("slc.nhit")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
bool isFHC
std::string OutDir
void CutThreshold(TH1D *MyHist, std::string Key)
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
correl_xv GetYaxis() -> SetDecimals()
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kAna2017Period5POT
Definition: Exposures.h:171
void CalcCutValsPlot(std::string sIdentifier, bool fhc)
double POTNorm
TH2 * Rebin(TH2 *hist, std::vector< double > edges)
OStream cerr
Definition: OStream.cxx:7
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
osc::OscCalcDumb calc
const double kAna2017Period3Livetime
Definition: Exposures.h:199
TH1D * GetHistVectors(TFile *InF, std::string InDir)
#define M_PI
Definition: SbMath.h:34
const double kAna2018RHCPOT
Definition: Exposures.h:208
double LivNorm
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
correl_xv GetXaxis() -> SetDecimals()
leg AddEntry(GRdata,"data","p")
hmean SetLineColor(4)
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void MakeCanvas(TH1D *Period1, TH1D *Period2, TH1D *FHCHigh, TH1D *RHCHigh, bool isBPF, std::string Ident)
void SetTh23(const T &th23) override
const double kAna2017Period5Livetime
Definition: Exposures.h:198
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
std::string sFHC
const double kSecondAnaPeriod1Livetime
Definition: Exposures.h:99
void SetdCP(const T &dCP) override
bool isBPF
const double kSecondAnaPeriod2Livetime
Definition: Exposures.h:100
TFile * OutFile
cosmicTree SaveAs("cosmicTree.root")
Float_t e
Definition: plot.C:35
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
const double kAna2017Period3POT
Definition: Exposures.h:172
c cd(1)
double CutAim
gPad SetLogy()
const double kAna2018RHCLivetime
Definition: Exposures.h:214
T asin(T number)
Definition: d0nt_math.hpp:60
gm Write()
enum BeamMode string