plotWrongSignFraction.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////
2 // //
3 // WRONG SIGN FRACTION PLOT //
4 // USES RESULT OF CVN AND BDT BASED WS ESTIMATES IN NUE WG //
5 // (see note, docdb29919) //
6 // //
7 // USES RESULT OF NEUTRON CAPTURE STUDY FOR NUMU AS WELL //
8 // (see note, docdb22955, esp v9) //
9 /////////////////////////////////////////////////////////////
10 
13 #include "CAFAna/Fit/Fit.h"
15 #include "CAFAna/Analysis/Style.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "CAFAna/Core/Spectrum.h"
19 #include "CAFAna/Core/Utilities.h"
20 #include "CAFAna/Analysis/Calcs.h"
27 #include "CAFAna/Core/Progress.h"
34 #include "OscLib/OscCalcPMNSOpt.h"
35 #include "OscLib/OscCalcDumb.h"
36 #include "Utilities/func/MathUtil.h"
38 
39 using namespace ana;
40 
41 #include "TFile.h"
42 #include "TH2.h"
43 #include "TMath.h"
44 #include "TCanvas.h"
45 #include "TGraph.h"
46 #include "TLatex.h"
47 #include "TLegend.h"
48 #include "TBox.h"
49 #include "TLine.h"
50 #include "TSystem.h"
51 #include "TColor.h"
52 #include "TPaletteAxis.h"
53 #include "THStack.h"
54 #include "TString.h"
55 #include "TEllipse.h"
56 #include "TGraphErrors.h"
57 #include "TPaveText.h"
58 #include "TColor.h"
59 
60 #include <iostream>
61 #include <iomanip>
62 #include <string>
63 #include <vector>
64 #include "Utilities/rootlogon.C"
65 
67 {
68  // Setup
69  bool useMultiverse = true;
70  float hMinYE = 0.;
71  float hMaxYE = 1200.;
72  int nBinYE = ceil(hMaxYE-hMinYE);
73 
74  float hMinYM = hMinYE;
75  float hMaxYM = hMaxYE;
76  int nBinYM = ceil(hMaxYM-hMinYM);
77 
78  float hMinYE2 = 0.;
79  float hMaxYE2 = 1200.;
80  int nBinYE2 = ceil(hMaxYE-hMinYE);
81 
82  float hMinXE = 0.;
83  float hMaxXE = 0.45;
84  int nBinXE = 18;
85 
86  float hMinXM = 0.;
87  float hMaxXM = 0.2;
88  int nBinXM = 8;
89 
90  float hMinXE2 = 0.;
91  float hMaxXE2 = 0.45;
92  int nBinXE2 = 18;
93 
94  // MC Multiverse inputs
95  TH1D *multiverseNueLoCVN;
96  TH1D *multiverseNueHiCVN;
97  TH1D *multiverseNueMerge;
98  TH1D *multiverseNumu;
99  if( useMultiverse ){
100  TFile *inFileMultiverse = new TFile("/nova/app/users/mylab/NUE/ws/moreUniv/wsf_ppfxgenieBF_multiU.root","read");
101  multiverseNueLoCVN = (TH1D*)inFileMultiverse->Get("rhc_totNDnue_lowCVN");
102  multiverseNueLoCVN->SetDirectory(0);
103  multiverseNueLoCVN->SetLineColor(kBlue-9);
104  multiverseNueLoCVN->SetFillColor(kBlue-9);
105  multiverseNueHiCVN = (TH1D*)inFileMultiverse->Get("rhc_totNDnue_highCVN");
106  multiverseNueHiCVN->SetDirectory(0);
107  multiverseNueHiCVN->SetLineColor(kBlue-9);
108  multiverseNueHiCVN->SetFillColor(kBlue-9);
109  multiverseNueMerge = (TH1D*)inFileMultiverse->Get("rhc_totNDnue_tot");
110  multiverseNueMerge->SetDirectory(0);
111  multiverseNueMerge->SetLineColor(kBlue-9);
112  multiverseNueMerge->SetFillColor(kBlue-9);
113  multiverseNumu = (TH1D*)inFileMultiverse->Get("rhc_totNDnumu_tot");
114  multiverseNumu->SetDirectory(0);
115  multiverseNumu->SetLineColor(kBlue-9);
116  multiverseNumu->SetFillColor(kBlue-9);
117  inFileMultiverse->Close();
118 
119  // Binning
120  // Nue plots
121  int nBinXEnom = multiverseNueLoCVN->GetXaxis()->GetNbins();
122  if( nBinXEnom!=multiverseNueHiCVN->GetXaxis()->GetNbins() ) std::cout << "Nue plots have different x bins!" << std::endl;
123  if( nBinXEnom == 1 ) std::cout << "ONLY ONE BIN DETECTED!! This will produce nasty results..." << std::endl;
124  hMinXE = multiverseNueLoCVN->GetXaxis()->GetBinLowEdge(1);
125  int nBinXE = nBinXEnom+20;
126  hMaxXE = hMinXE + (nBinXE*(multiverseNueLoCVN->GetXaxis()->GetBinLowEdge(2) - multiverseNueLoCVN->GetXaxis()->GetBinLowEdge(1)));
127  hMaxYE = 1.1*std::max(multiverseNueLoCVN->GetMaximum(),multiverseNueHiCVN->GetMaximum());
128 
129  // Hi CVN
130  nBinXE2 = nBinXEnom*(float(multiverseNueHiCVN->GetMaximumBin())/float(multiverseNueLoCVN->GetMaximumBin())) + 18;
131  hMaxYE2 = hMaxYE;
132  hMinXE2 = multiverseNueHiCVN->GetXaxis()->GetBinLowEdge(1);
133  hMaxXE2 = hMinXE2 + (nBinXE2*(multiverseNueHiCVN->GetXaxis()->GetBinLowEdge(2) - multiverseNueHiCVN->GetXaxis()->GetBinLowEdge(1)));
134 
135  // Numu plots
136  hMaxYM = 1.1*multiverseNumu->GetMaximum();
137  nBinXM = 42;
138  hMinXM = multiverseNumu->GetXaxis()->GetBinLowEdge(1);
139  hMaxXM = hMinXM + (nBinXM*(multiverseNumu->GetXaxis()->GetBinLowEdge(2) - multiverseNumu->GetXaxis()->GetBinLowEdge(1)));
140  }
141  // ---------------------------
142 
143  TH2F *blankHistE = new TH2F("blankHistE",";Wrong-sign fraction, low CVN #nu_{e} selection;Fake experiments",nBinXE,hMinXE,hMaxXE,nBinYE,hMinYE,hMaxYE);
144  CenterTitles(blankHistE);
145  TH2F *blankHistE2=new TH2F("blankHistE2",";Wrong-sign fraction, high CVN #nu_{e} selection;Fake experiments",nBinXE2,hMinXE2,hMaxXE2,nBinYE2,hMinYE2,hMaxYE2);
146  CenterTitles(blankHistE2);
147 
148  TH2F *blankHistM = new TH2F("blankHistM",";Wrong-sign fraction in #nu_{#mu} selection;Fake experiments",nBinXM,hMinXM,hMaxXM,nBinYM,hMinYM,hMaxYM);
149  CenterTitles(blankHistM);
150 
151  // this axis has no meaning in the data-driven results, so just offset for visibility
152  float yVal[] = {150.,250.,350.};
153  float yErr[] = { 0., 0., 0.};
154  float yValNumu[] = {float(150.)*float(hMaxYM/hMaxYE),float(250.)*float(hMaxYM/hMaxYE),float(350.)*float(hMaxYM/hMaxYE)};
155 
156  // Text boxes
157  TPaveText *nueText = new TPaveText(0.12,0.195,0.375,0.465,"NB NDC");
158  TText *nueTextL1 = nueText->AddText("Wrong-sign BDT");
159  nueTextL1->SetTextAlign(12);
160  TText *nueTextL2 = nueText->AddText("Prong CVN Proton ID");
161  nueTextL2->SetTextAlign(12);
162  TText *nueTextL3 = nueText->AddText("Event CVN Proton ID");
163  nueTextL3->SetTextAlign(12);
164  nueText->SetFillStyle(3155);
165 
166  // For using only the one Numu result
167  TPaveText *numu1Text = new TPaveText(0.12,0.195,0.375,0.285,"NB NDC");
168  TText *numu1TextL1 = numu1Text->AddText("Neutron Capture Rate");
169  numu1TextL1->SetTextAlign(12);
170  numu1Text->SetFillStyle(3155);
171 
172  // For the MC Universes
173  TPaveText *addedTextMC1;
174  TPaveText *addedTextMC2;
175  TPaveText *addedTextMC3;
176  if( useMultiverse ){
177  addedTextMC1 = new TPaveText(0.63, 0.644, 0.8725, 0.69625, "NB NDC"); // xlow 0.6275, ylow 0.65, xhigh 0.8725, yhigh 0.705
178  addedTextMC2 = new TPaveText(0.63, 0.585, 0.88, 0.67, "NB NDC");
179  addedTextMC3 = new TPaveText(0.63, 0.56, 0.875, 0.615, "NB NDC"); // xlow 0.625, ylow 0.54
180  TText *addedText1 = addedTextMC1->AddText("Flux and cross-");
181  TText *addedText2 = addedTextMC2->AddText("section systematic");
182  TText *addedText3 = addedTextMC3->AddText("fake experiments");
183  addedText1->SetTextColor(kBlue-7);
184  addedText2->SetTextColor(kBlue-7);
185  addedText3->SetTextColor(kBlue-7);
186  addedTextMC1->SetFillStyle(3155);
187  addedTextMC2->SetFillStyle(3155);
188  addedTextMC3->SetFillStyle(3155);
189  }
190 
191  // Nue results: see docDB-27899 and technote above
192  // These are found by running the calculate macro
193  // x-value is WS Fraction
194  // ------------------------------------------------------------------------------------
195  // Low CVN electron bin (>=0.89 and <0.98)
196  float xValNueLo[] = {0.2866, 0.2754, 0.3499};
197  float xErrNueLo[] = {0.0364, 0.0730, 0.0334};
198  // High CVN electron bin (>=0.98)
199  float xValNueHi[] = {0.2423, 0.1837, 0.2308};
200  float xErrNueHi[] = {0.0284, 0.0551, 0.0271};
201 
202  // Some values from running script
203  // NUE Selection nominal (CV Weights) composition {LoCVN, HiCVN}
204  /* HICVN WSCC: 236.266 WSNC: 20.5481 RSCC: 874.598 RSNC: 43.2865 */
205  /* LOCVN WSCC: 556.763 WSNC: 165.225 RSCC: 1105.84 RSNC: 356.042 */
206  float nueSelNue[] = {1475.42,1069.5};
207  float nueSelNumu[] = {187.188,41.3593};
208  float nueSelNC[] = {521.267,63.8345};
209 
210  float nueLoFractionNue = nueSelNue[0] / (nueSelNue[0] + nueSelNumu[0] + nueSelNC[0] );
211  float nueLoFractionNumu = nueSelNumu[0] / (nueSelNue[0] + nueSelNumu[0] + nueSelNC[0] );
212  float nueLoFractionNC = nueSelNC[0] / (nueSelNue[0] + nueSelNumu[0] + nueSelNC[0] );
213  float nueHiFractionNue = nueSelNue[1] / (nueSelNue[1] + nueSelNumu[1] + nueSelNC[1] );
214  float nueHiFractionNumu = nueSelNumu[1] / (nueSelNue[1] + nueSelNumu[1] + nueSelNC[1] );
215  float nueHiFractionNC = nueSelNC[1] / (nueSelNue[1] + nueSelNumu[1] + nueSelNC[1] );
216 
217  TGraphErrors *graphNueLoCVN = new TGraphErrors(3,xValNueLo,yVal,xErrNueLo,yErr);
218  graphNueLoCVN->SetMarkerStyle(24);
219  graphNueLoCVN->SetLineWidth(2);
220 
221  TGraphErrors *graphNueHiCVN = new TGraphErrors(3,xValNueHi,yVal,xErrNueHi,yErr);
222  graphNueHiCVN->SetMarkerStyle(24);
223  graphNueHiCVN->SetLineWidth(2);
224  // ------------------------------------------------------------------------------------
225 
226  // Numu Results: see notes listed at top for more info
227  // n capture study:
228  // Take scale factors from study and recast them as WS Fraction in numu selection
229  // ------------------------------------------------------------------------------------
230  // -----------------------
231  // --- n capture study ---
232  // -----------------------
233  float scaleWS = 1.08;
234  float errorWS = 0.14;
235  float scaleNC = 1.0;
236  float errorNC = 0.29;
237  // Then taking numu prediction numbers from info used to make docDB-27899
238  // But you get basically the same numbers (up to rounding) by using
239  // committed macros for each quantile and summing
240  float WSCC = 41208.9; // all selected CC nubar events
241  float WSNC = 356.544; // all selected NC nubar events
242  float RSCC = 346589; // all selected CC nu events
243  float RSNC = 680.572; // all selected NC nu events
244  // All selected numu data at the same POT scaling.
245  float numuData = 386485;
246 
247  // Define scaleRS as the value which makes the predictions work with the given scaleWS and scaleNC
248  // data = scaleWS*WSCC + scaleNC*(WSNC+RSNC) + scaleRS*RSCC
249  float scaleRS = (numuData - scaleWS*WSCC - scaleNC*(WSNC+RSNC))/RSCC;
250  std::cout << "For neutron capture study with current MC predictions, scaleRS is found to be " << scaleRS << std::endl;
251 
252  // Few other parameters of interest for plots
253  // as above, you get the same to within rounding errors
254  // if you use committed macros for each quantile and sum
255  float numuSelNue = 12.0812;
256  float numuSelNumu = 387786;
257  float numuSelNC = 1037.12;
258 
259  float numuFractionNue = numuSelNue / (numuSelNue + numuSelNumu + numuSelNC);
260  float numuFractionNumu = numuSelNumu / (numuSelNue + numuSelNumu + numuSelNC);
261  float numuFractionNC = numuSelNC / (numuSelNue + numuSelNumu + numuSelNC);
262 
263  // Calculate WS Fraction
264  float denomNumu = (scaleWS*WSCC) + (scaleNC*(WSNC+RSNC)) + (scaleRS*RSCC);
265  float denomsq = denomNumu*denomNumu;
266  float xValNumu[] = { ((scaleWS*WSCC) + (scaleNC*WSNC)) / denomNumu };
267  // covariance 0 because not considering covariance in current handling with pull term
268  float xErrNumu[] = { sqrt( errorWS*errorWS*(WSCC*(scaleNC*RSNC + scaleRS*RSCC)/denomsq)*(WSCC*(scaleNC*RSNC + scaleRS*RSCC)/denomsq) +
269  errorNC*errorNC*((scaleRS*RSCC*WSNC - scaleWS*WSCC*RSNC)/denomsq)*((scaleRS*RSCC*WSNC - scaleWS*WSCC*RSNC)/denomsq) ) };
270  // ------------------------------------------------------------------------------------
271 
272  // Put together plots
273  // ------------------------------------------------------------------------------------
274  TGraphErrors *graphNumu = new TGraphErrors(1,xValNumu,yValNumu,xErrNumu,yErr);
275  graphNumu->SetMarkerStyle(24);
276  graphNumu->SetLineWidth(2);
277 
278  // Pave Texts for composition fractions
279  TPaveText *typeCompText = new TPaveText(0.24,0.66,0.34,0.84,"NB NDC");
280  TText *typeCompTextL1 = typeCompText->AddText( "#nu_{e} + #bar{#nu}_{e}" );
281  typeCompTextL1->SetTextAlign(12);
282  TText *typeCompTextL2 = typeCompText->AddText( "#nu_{#mu} + #bar{#nu}_{#mu}" );
283  typeCompTextL2->SetTextAlign(12);
284  TText *typeCompTextL3 = typeCompText->AddText( "NC" );
285  typeCompTextL3->SetTextAlign(12);
286  typeCompText->SetFillStyle(3155);
287 
288  TPaveText *nueLoCompText = new TPaveText(0.14,0.66,0.24,0.84,"NB NDC"); // had used ylo = 0.57
289  TText *nueLoCompTextL1 = nueLoCompText->AddText( TString::Format("%.0f%%",100.*nueLoFractionNue) );
290  nueLoCompTextL1->SetTextAlign(32); //right-align
291  TText *nueLoCompTextL2 = nueLoCompText->AddText( TString::Format("%.0f%%",100.*nueLoFractionNumu) );
292  nueLoCompTextL2->SetTextAlign(32);
293  TText *nueLoCompTextL3 = nueLoCompText->AddText( TString::Format("%.0f%%",100.*nueLoFractionNC) );
294  nueLoCompTextL3->SetTextAlign(32);
295  nueLoCompText->SetFillStyle(3155);
296 
297  TPaveText *nueHiCompText = new TPaveText(0.14,0.66,0.24,0.84,"NB NDC"); // had used ylo = 0.57
298  TText *nueHiCompTextL1 = nueHiCompText->AddText( TString::Format("%.0f%%",100.*nueHiFractionNue) );
299  nueHiCompTextL1->SetTextAlign(32);
300  TText *nueHiCompTextL2 = nueHiCompText->AddText( TString::Format("%.0f%%",100.*nueHiFractionNumu) );
301  nueHiCompTextL2->SetTextAlign(32);
302  TText *nueHiCompTextL3 = nueHiCompText->AddText( TString::Format("%.0f%%",100.*nueHiFractionNC) );
303  nueHiCompTextL3->SetTextAlign(32);
304  nueHiCompText->SetFillStyle(3155);
305 
306  TPaveText *numuCompText = new TPaveText(0.14,0.66,0.24,0.84,"NB NDC"); // had used ylo = 0.57
307  TText *numuCompTextL1 = numuCompText->AddText( TString::Format("%.1f%%",100.*numuFractionNue) ); // try 1 decimal place
308  numuCompTextL1->SetTextAlign(32);
309  TText *numuCompTextL2 = numuCompText->AddText( TString::Format("%.1f%%",100.*numuFractionNumu) );
310  numuCompTextL2->SetTextAlign(32);
311  TText *numuCompTextL3 = numuCompText->AddText( TString::Format("%.1f%%",100.*numuFractionNC) );
312  numuCompTextL3->SetTextAlign(32);
313  numuCompText->SetFillStyle(3155);
314 
315  // Draw plots
316  TCanvas *cNueLoCVN = new TCanvas();
317  blankHistE->Draw("");
318  if( useMultiverse ){
319  multiverseNueLoCVN->Draw("same hist ][");
320  addedTextMC1->Draw();
321  addedTextMC2->Draw();
322  addedTextMC3->Draw();
323  }
324  graphNueLoCVN->Draw("same P");
325  nueText->Draw();
326  typeCompText->Draw();
327  nueLoCompText->Draw();
328  blankHistE->Draw("axissame");
329  Preliminary();
330  cNueLoCVN->Print("cNueLoCVN_Ana2018_WSFraction.pdf");
331 
332  TCanvas *cNueHiCVN = new TCanvas();
333  blankHistE2->Draw("");
334  if( useMultiverse ){
335  multiverseNueHiCVN->Draw("same hist ][");
336  addedTextMC1->Draw();
337  addedTextMC2->Draw();
338  addedTextMC3->Draw();
339  }
340  graphNueHiCVN->Draw("same P");
341  nueText->Draw();
342  typeCompText->Draw();
343  nueHiCompText->Draw();
344  blankHistE2->Draw("axissame");
345  Preliminary();
346  cNueHiCVN->Print("cNueHiCVN_Ana2018_WSFraction.pdf");
347 
348  TCanvas *cNumu = new TCanvas();
349  blankHistM->Draw("");
350  if( useMultiverse ){
351  multiverseNumu->Draw("same hist ][");
352  addedTextMC1->Draw();
353  addedTextMC2->Draw();
354  addedTextMC3->Draw();
355  }
356  graphNumu->Draw("same P");
357  numu1Text->Draw();
358  typeCompText->Draw();
359  numuCompText->Draw();
360  blankHistM->Draw("axissame");
361  Preliminary();
362  cNumu->Print("cNumu_Ana2018_WSFraction.pdf");
363 
364  std::cout << std::endl;
365  std::cout << std::endl;
366  std::cout << "===================================================" << std::endl;
367  std::cout << " Wrong-sign Fraction " << std::endl;
368  std::cout << "===================================================" << std::endl;
369  std::cout << " N Capture Numu WS Fraction " << xValNumu[0] << " +/- " << xErrNumu[0] << std::endl;
370  std::cout << "===================================================" << std::endl;
371  std::cout << std::endl;
372  std::cout << std::endl;
373 
374  std::cout << "Done!" << std::endl;
375 }
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
T sqrt(T number)
Definition: d0nt_math.hpp:156
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
void plotWrongSignFraction()
void Preliminary()
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
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11