ND_pred_MC_nue_signal.C
Go to the documentation of this file.
1 //Program writen by Prabhjot Singh on August 11, 2016 at 12:15 pm.
2 //
3 //Purpose: Make ND prediction over MC histograms and their ratios.
4 // ND prediction over MC histograms helps to make FD prediction for the nue signal (numu to nue signal)
5 //
6 //
7 
8 #include<Header.h>
9 #include<HistogramAttr.h>
10 
11 //root libraries
12 #include "vector.h"
13 
15 
16  //disable stat box.
17  gStyle->SetOptStat(" ");
18 
19  TString det_type[] = {"Near"};
20  TString interaction_type[] = {"NuMuCC"};
21  TString selection_type[] = {"NuMuSel"};
22  TString file_type[] = {"Beam"};
23 
24  bool isMC;
25 
26  TString epoch;
27  TString dirPOT;
28  TString filenamePOT;
29  TString POTfile;
30  TString dir;
31  TString filename;
32  TString rootfile;
33  TString histdir;
34  TString histname_reco;
35  TString histname_true;
36  TString histname_recotrue;
37  TString MC;
38  TString Data;
39 
40  Double_t nd_mc_pot;
41  Double_t nd_data_pot;
42 
43  Float_t margin = 0.15;
44 
45  //directory path where the root file exists with the spill and POT summary
46  dirPOT = "/nova/ana/users/prabhjot/FNEX/nue_analysis/";
47 
48  //name of the root file that has spill and POT summary
49  filenamePOT = "eventlisttrees_prod_pid_R16-03-03-prod2reco.a_fhc.root";
50 
51  POTfile = dirPOT+filenamePOT;
52 
53  //ND MC and Data POT
54  nd_mc_pot = POTinformationMC();
55  nd_data_pot = POTinformationData();
56 
57  //directoy path where the root file with histograms exists
58  dir = "/nova/ana/users/prabhjot/FNEX/nue_analysis/";
59 
60  //name of the root file that has all the histograms
61  filename = "fnex_plotpoint_hist_nueanalysis.root";
62 
63  rootfile = dir+filename;
64 
65  //name of the histogram directory
66  histdir = "plot/SkimmedListsConcatenatedByEpoch/";
67 
68  //initial name of the reco energy histograms
69  histname_reco = "neutrinoE";
70 
71  //initial name of the true energy histograms
72  histname_true = "neutrinoTrueE";
73 
74  //initial name of the reco vs true histogram
75  histname_recotrue = "recoToTrue";
76 
77  //combined epoch name in the histogram names
78  epoch = "Epoch0z";
79 
80  //MC histogram
81  MC = "MC";
82 
83  //Data histogram
84  Data = "Data";
85 
86  //Open root file
87  TFile* f = new TFile(rootfile, "read");
88  f->cd();
89 
90  //ND data histogram that we want to use.
91  TH1F* histdata = (TH1F*)f->Get(histdir+histname_reco+Data+det_type[0]+selection_type[0]+epoch);
92 
93  //ND MC for the reco energy for Beam, numu CC, numu selection and near detector
94  TH1F* histMC_reco = (TH1F*)f->Get(histdir+histname_reco+MC+file_type[0]+interaction_type[0]+det_type[0]+selection_type[0]+epoch);
95 
96  //ND MC for the true energy for Beam, numu CC, numu selection and near detector
97  TH1F* histMC_true = (TH1F*)f->Get(histdir+histname_true+MC+file_type[0]+interaction_type[0]+det_type[0]+selection_type[0]+epoch);
98 
99  //ND MC reco. vs. true energy plot for Beam, numu CC, numu selection and near detector
100  TH2F* hist_recotrue = (TH2F*)f->Get(histdir+histname_recotrue+MC+file_type[0]+interaction_type[0]+det_type[0]+selection_type[0]+epoch);
101 
102  //Normalize data and MC histograms by POT
103  Double_t ratio_pot = histdata->Integral()/histMC_reco->Integral();
104  std::cout << "Normalization ratio = " << ratio_pot << std::endl;
105 
106  histMC_reco->Scale(ratio_pot);
107  histMC_true->Scale(ratio_pot);
108 
109  std::cout << " " << std::endl;
110  std::cout << "Histogram directory = " << histdir << std::endl;
111  std::cout << "Reco. Histogram name = " << histname_reco << std::endl;
112  std::cout << "True. Histogram name = " << histname_true << std::endl;
113  std::cout << "Data = " << Data << std::endl;
114  std::cout << "MC = " << MC << std::endl;
115  std::cout << "File Type = " << file_type[0] << std::endl;
116  std::cout << "Interaction type = " << interaction_type[0] << std::endl;
117  std::cout << "Detector type = " << det_type[0] << std::endl;
118  std::cout << "Selection type = " << selection_type[0] << std::endl;
119  std::cout << " " << std::endl;
120 
121  //Draw data histogram
122  new TCanvas;
123  gPad->SetLeftMargin(margin);
124  gPad->SetBottomMargin(margin);
125  histdata->Draw();
126  HistogramAttr1D(histdata, "E_{#nu}^{reco} (GeV)", "N^{Data}_{#nu_{#mu}, S_{#mu}} Events");
127  // histdata->GetXaxis()->SetRangeUser(0, 5);
128  histdata->SetLineWidth(2);
129  histdata->SetLineColor(kBlack);
130  //Set legends for the data plot
131  TLegend *leg_data = new TLegend(.60,.73,.78,.84);
132  leg_data->SetBorderSize(0);
133  leg_data->SetFillStyle(0);
134  leg_data->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
135  leg_data->AddEntry(histdata, Data);
136  leg_data->SetTextSize(.05);
137  leg_data->Draw();
138  gPad->Print("Images/NDdata_recoE.pdf");
139 
140  //Draw reco energy MC histogram
141  new TCanvas;
142  gPad->SetLeftMargin(margin);
143  gPad->SetBottomMargin(margin);
144  histMC_reco->Draw();
145  HistogramAttr1D(histMC_reco, "E_{#nu}^{reco} (GeV)", "N^{MC}_{#nu_{#mu}, S_{#mu}} Events");
146  // histMC_reco->GetXaxis()->SetRangeUser(0, 5);
147  histMC_reco->SetLineWidth(2);
148  histMC_reco->SetLineColor(kRed);
149  //Set legends for the MC reco energy plot
150  TLegend *leg_MC_reco = new TLegend(.60,.73,.78,.84);
151  leg_MC_reco->SetBorderSize(0);
152  leg_MC_reco->SetFillStyle(0);
153  leg_MC_reco->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
154  leg_MC_reco->AddEntry(histMC_reco, MC);
155  leg_MC_reco->SetTextSize(.05);
156  leg_MC_reco->Draw();
157  gPad->Print("Images/NDMC_recoE.pdf");
158 
159  //Draw true energy MC histogram
160  new TCanvas;
161  gPad->SetLeftMargin(margin);
162  gPad->SetBottomMargin(margin);
163  histMC_true->Draw();
164  HistogramAttr1D(histMC_true, "E_{#nu}^{true} (GeV)", "N^{MC}_{#nu_{#mu}, S_{#mu}} Events");
165  // histMC_true->GetXaxis()->SetRangeUser(0, 5);
166  histMC_true->SetLineWidth(2);
167  histMC_true->SetLineColor(kBlue);
168  //Set legends for the MC true energy plot
169  TLegend *leg_MC_true = new TLegend(.60,.73,.78,.84);
170  leg_MC_true->SetBorderSize(0);
171  leg_MC_true->SetFillStyle(0);
172  leg_MC_true->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
173  leg_MC_true->AddEntry(histMC_true, MC);
174  leg_MC_true->SetTextSize(.05);
175  leg_MC_true->Draw();
176  gPad->Print("Images/NDMC_trueE.pdf");
177 
178  //Draw 2D reco. vs. true energy histogram.
179  new TCanvas;
180  gPad->SetLeftMargin(margin);
181  gPad->SetBottomMargin(margin);
182  gPad->SetRightMargin(1.2*margin);
183  gPad->SetLogz();
184  // hist_recotrue->Draw("colz text");
185  HistogramAttr2D(hist_recotrue, "E_{#nu}^{true} (GeV)", "E_{#nu}^{reco} (GeV)", "N^{MC}_{#nu_{#mu}, S_{#mu}} Events");
186  hist_recotrue->Scale(ratio_pot);
187  gPad->Print("Images/NDMC_recoEtrueE.pdf");
188  // gPad->Print("Images/ND_truevsreco_old.C");
189 
190  //Draw data and MC histograms on the same canvas as well.
191  new TCanvas;
192  gPad->SetLeftMargin(margin);
193  gPad->SetBottomMargin(margin);
194  TH1F* histdata_clone = (TH1F*)histdata->Clone("histdata_clone");
195  TH1F* histMC_reco_clone = (TH1F*)histMC_reco->Clone("histMC_reco_clone");
196  TH1F* histMC_true_clone = (TH1F*)histMC_true->Clone("histMC_true_clone");
197  histdata_clone->Draw();
198  histMC_reco_clone->Draw("same");
199  histMC_true_clone->Draw("same");
200  HistogramAttr1D(histdata_clone, "E_{#nu} (GeV)", "Events");
201  // histdata_clone->GetYaxis()->SetRangeUser(0, 100e3);
202  //legends for the data and MC plots
203  TLegend *leg_data_MC = new TLegend(.53,.58,.73,.85);
204  leg_data_MC->SetBorderSize(0);
205  leg_data_MC->SetFillStyle(0);
206  leg_data_MC->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
207  leg_data_MC->AddEntry(histdata_clone, Data);
208  leg_data_MC->AddEntry(histMC_reco_clone, MC+" "+(histMC_reco_clone->GetXaxis()->GetTitle()));
209  leg_data_MC->AddEntry(histMC_true_clone, MC+" "+(histMC_true_clone->GetXaxis()->GetTitle()));
210  leg_data_MC->SetTextSize(.05);
211  leg_data_MC->Draw();
212  gPad->Print("Images/NDSpec.pdf");
213 
214  //Take ratio of N(data)/N(MC) as a function of reco energy.
215  TH1F* ratio_dataByMC = ratio(histMC_reco, histdata);
216  new TCanvas;
217  gPad->SetLeftMargin(margin);
218  gPad->SetBottomMargin(margin);
219  ratio_dataByMC->Draw("E");
220  ratio_dataByMC->SetMarkerStyle(20);
221  HistogramAttr1D(ratio_dataByMC, (histdata->GetXaxis()->GetTitle()), "N^{Data}_{#nu_{#mu}, S_{#mu}} / N^{MC}_{#nu_{#mu}, s_{#mu}}");
222  gPad->Print("Images/NDdataMCratio_recoE.pdf");
223 
224  //make a modify 2D N(MC) plot by multiplying the original 2D ND(MC) plot by the ratio of the N(data)/N(MC) as a function of reco energy.
225  TH2F* hist_recotrue_modified = modify2D(hist_recotrue, ratio_dataByMC);
226  new TCanvas;
227  gPad->SetLogz();
228  gPad->SetLeftMargin(margin);
229  gPad->SetBottomMargin(margin);
230  gPad->SetRightMargin(1.2*margin);
231  hist_recotrue_modified->Draw("colz text");
232  // hist_recotrue_modified->GetXaxis()->SetRangeUser(0, 5);
233  // hist_recotrue_modified->GetYaxis()->SetRangeUser(0, 5);
234  HistogramAttr2D(hist_recotrue_modified, "E_{#nu}^{true} (GeV)", "E_{#nu}^{reco} (GeV)", "N^{MC}_{#nu_{#mu}, S_{#mu}} Events");
235  gPad->Print("Images/NDSpec_weighted.pdf");
236  // gPad->Print("Images/ND_truevsreco_modified_old.C");
237 
238  //make a ND prediction as a function of true energy.
239  TH1F* ND_prediction = NDPrediction(hist_recotrue_modified);
240  new TCanvas;
241  gPad->SetLogz();
242  gPad->SetLeftMargin(margin);
243  gPad->SetBottomMargin(margin);
244  ND_prediction->Draw("HIST");
245  // ND_prediction->GetXaxis()->SetRangeUser(0, 5);
246  HistogramAttr1D(ND_prediction, hist_recotrue_modified->GetXaxis()->GetTitle(), "N^{Prediction}_{#nu_{#mu}, S_{#mu}} Events");
247  ND_prediction->SetLineColor(28);
248  //Set legends for the ND prediction true energy plot
249  TLegend* leg_NDpred_true = new TLegend(.60,.73,.78,.84);
250  leg_NDpred_true->SetBorderSize(0);
251  leg_NDpred_true->SetFillStyle(0);
252  leg_NDpred_true->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
253  leg_NDpred_true->AddEntry(ND_prediction, "Prediction");
254  leg_NDpred_true->SetTextSize(.05);
255  leg_NDpred_true->Draw();
256  gPad->Print("Images/NDPred_trueE.pdf");
257 
258  //take ratio of ND prediction and ND MC as a function of true energy.
259  TH1F* ratio_ND_predictionToMC = ratio(histMC_true, ND_prediction);
260  new TCanvas;
261  gPad->SetLeftMargin(margin);
262  gPad->SetBottomMargin(margin);
263  ratio_ND_predictionToMC->Draw();
264  HistogramAttr1D(ratio_ND_predictionToMC, "E_{#nu}^{true} (GeV)", "N^{Prediction}_{#nu_{#mu}, S_{#mu}} / N^{MC}_{#nu_{#mu}, s_{#mu}}");
265  // ratio_ND_predictionToMC->GetXaxis()->SetRangeUser(0, 5);
266  ratio_ND_predictionToMC->SetLineWidth(2);
267  ratio_ND_predictionToMC->SetLineColor(28);
268  ratio_ND_predictionToMC->SetMarkerStyle(20);
269  ratio_ND_predictionToMC->SetMarkerColor(28);
270  gPad->Print("Images/NDPredMCratio_trueE.pdf");
271 
272  //plot both ratios N(data)/N(MC) and N(pred)/N(MC) in same canvas.
273  TH1F* ratio_dataByMC_clone = (TH1F*)ratio_dataByMC->Clone("ratio_dataByMC_clone");
274  TH1F* ratio_ND_predictionToMC_clone = (TH1F*)ratio_ND_predictionToMC->Clone("ratio_ND_predictionToMC_clone");
275  new TCanvas;
276  gPad->SetLeftMargin(margin);
277  gPad->SetBottomMargin(margin);
278  ratio_dataByMC_clone->Draw();
279  HistogramAttr1D(ratio_dataByMC_clone, "E_{#nu} (GeV)", "Ratio");
280  ratio_ND_predictionToMC_clone->Draw("same");
281  //legends for ratios N(data)/N(MC) and N(pred)/N(MC)
282  TLegend* leg_ratios = new TLegend(.31,.61,.51,.88);
283  leg_ratios->SetBorderSize(0);
284  leg_ratios->SetFillStyle(0);
285  leg_ratios->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
286  leg_ratios->AddEntry(ratio_dataByMC_clone, (ratio_dataByMC->GetYaxis()->GetTitle()), "lep");
287  leg_ratios->AddEntry(ratio_ND_predictionToMC_clone, (ratio_ND_predictionToMC_clone->GetYaxis()->GetTitle()), "lep");
288  leg_ratios->SetTextSize(.05);
289  leg_ratios->Draw();
290  gPad->Print("Images/ND_doubleratio.pdf");
291 
292  //plot double ratio of N(data)/N(MC) to N(pred)/N(MC)
293  TH1F* double_ratio = ratio(ratio_dataByMC_clone, ratio_ND_predictionToMC_clone);
294  TH1F* double_ratio = ratio(ratio_ND_predictionToMC_clone, ratio_dataByMC_clone);
295  new TCanvas;
296  gPad->SetLeftMargin(margin);
297  gPad->SetBottomMargin(margin);
298  double_ratio->Draw();
299  double_ratio->SetLineColor(41);
300  double_ratio->SetMarkerColor(41);
301  // double_ratio->GetYaxis()->SetRangeUser(0.7, 1.4);
302  HistogramAttr1D(double_ratio, "E_{#nu} (GeV)", "Double Ratio");
303  TLegend* leg_doubleratios = new TLegend(.31,.61,.51,.88);
304  leg_doubleratios->SetBorderSize(0);
305  leg_doubleratios->SetFillStyle(0);
306  leg_doubleratios->AddEntry((TObject*)0, det_type[0]+" "+"Detector"," ");
307  leg_doubleratios->AddEntry(double_ratio, "#frac{N^{data}/N^{MC}}{N^{Pred.}/N^{MC}}_{#nu_{#mu}, S_{#mu}}", "lep");
308  leg_doubleratios->SetTextSize(.05);
309  leg_doubleratios->Draw();
310  gPad->Print("Images/NDdouble_doubleratio.pdf");
311 
312  //take ratios of 2D histograms of N(MC) nu mu events
313  TH2F* ratio_2D_NDMC = ratio2D(hist_recotrue,hist_recotrue_modified);
314  new TCanvas;
315  //gPad->SetLogz();
316  gPad->SetLeftMargin(margin);
317  gPad->SetBottomMargin(margin);
318  gPad->SetRightMargin(1.2*margin);
319  ratio_2D_NDMC->Draw("colz");
320  HistogramAttr2D(ratio_2D_NDMC, "E_{#nu}^{true} (GeV)", "E_{#nu}^{reco} (GeV)", "Ratio N^{MC}_{#nu_{#mu}, S_{#mu}} Events");
321  ratio_2D_NDMC->SetMinimum(0.8);
322  ratio_2D_NDMC->SetMaximum(2.0);
323  gPad->Print("Images/NDratio_weightedTounweighted.pdf");
324 
325  //make the exact same ND prediction/MC for true energy but this time just use the bin content of the ND data and ND MC truereco energy histograms.
326  ND_predictionbyMC(histdata, hist_recotrue);
327 
328  ND_predictionbyMC_usingbincontent(histdata, hist_recotrue);
329 }//end of program ND_pred_MC_nue_signal()
330 //---------------------------------------------------------------------------------//
enum BeamMode kRed
TH1 * ratio(TH1 *h1, TH1 *h2)
TH2F * modify2D(TH2F *h2D, TH1F *h1D, TString prediction)
ND_pred_MC_nue_signal()
string filename
Definition: shutoffs.py:106
interaction_type
Definition: Constants.h:110
void ND_predictionbyMC(TH1F *h1, TH2F *h2)
OStream cout
Definition: OStream.cxx:6
Double_t POTinformationMC()
det_type
Definition: Constants.h:584
TH1F * NDPrediction(TH2F *h)
void ND_predictionbyMC_usingbincontent(TH1F *hdata, TH2F *htruereco)
TH2F * ratio2D(TH2F *h1, TH2F *h2)
TDirectory * dir
Definition: macro.C:5
void HistogramAttr2D(TH2 *h, char *titlx_x, char *titlx_y, char *titlx_z, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy)
Double_t POTinformationData()
file_type
Definition: Constants.h:150
HistogramAttr1D(TH1 *h, char *title_x, char *title_y, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy, Color_t lcolor)
enum BeamMode kBlue