plotNumuWrongSign.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 // UPDATED RESULTS FOR 2019 & SIMPLER ERROR CALCULATION //
11 /////////////////////////////////////////////////////////////
12 
15 #include "CAFAna/Fit/Fit.h"
17 #include "CAFAna/Analysis/Style.h"
18 #include "CAFAna/Vars/FitVars.h"
19 #include "CAFAna/Core/Spectrum.h"
21 #include "CAFAna/Core/Utilities.h"
22 #include "CAFAna/Analysis/Calcs.h"
29 #include "CAFAna/Core/Progress.h"
36 #include "OscLib/OscCalcPMNSOpt.h"
37 #include "OscLib/OscCalcDumb.h"
38 #include "Utilities/func/MathUtil.h"
40 
41 using namespace ana;
42 
43 #include "TFile.h"
44 #include "TH2.h"
45 #include "TMath.h"
46 #include "TCanvas.h"
47 #include "TGraph.h"
48 #include "TLatex.h"
49 #include "TLegend.h"
50 #include "TBox.h"
51 #include "TLine.h"
52 #include "TSystem.h"
53 #include "TColor.h"
54 #include "TPaletteAxis.h"
55 #include "THStack.h"
56 #include "TString.h"
57 #include "TEllipse.h"
58 #include "TGraphErrors.h"
59 #include "TPaveText.h"
60 #include "TColor.h"
61 
62 #include <iostream>
63 #include <iomanip>
64 #include <string>
65 #include <vector>
66 #include "Utilities/rootlogon.C"
67 
69 {
70  // Setup
71  bool useMultiverse = true;
72  float hMinYE = 0.;
73  float hMaxYE = 1200.;
74  int nBinYE = ceil(hMaxYE-hMinYE);
75 
76  float hMinYM = hMinYE;
77  float hMaxYM = hMaxYE;
78  int nBinYM = ceil(hMaxYM-hMinYM);
79 
80  float hMinYE2 = 0.;
81  float hMaxYE2 = 1200.;
82  int nBinYE2 = ceil(hMaxYE-hMinYE);
83 
84  float hMinXE = 0.;
85  float hMaxXE = 0.45;
86  int nBinXE = 18;
87 
88  float hMinXM = 0.;
89  float hMaxXM = 0.2;
90  int nBinXM = 8;
91 
92  float hMinXE2 = 0.;
93  float hMaxXE2 = 0.45;
94  int nBinXE2 = 18;
95 
96  // MC Multiverse inputs
97  TH1D *multiverseNueLoCVN;
98  TH1D *multiverseNueHiCVN;
99  TH1D *multiverseNueMerge;
100  TH1D *multiverseNumu;
101  if( useMultiverse ){
102  TFile *inFileMultiverse = new TFile("/nova/app/users/mylab/NUE/ws/moreUniv/wsf_ppfxgenieBF_multiU.root","read");
103  multiverseNueLoCVN = (TH1D*)inFileMultiverse->Get("rhc_totNDnue_lowCVN");
104  multiverseNueLoCVN->SetDirectory(0);
105  multiverseNueLoCVN->SetLineColor(kBlue-9);
106  multiverseNueLoCVN->SetFillColor(kBlue-9);
107  multiverseNueHiCVN = (TH1D*)inFileMultiverse->Get("rhc_totNDnue_highCVN");
108  multiverseNueHiCVN->SetDirectory(0);
109  multiverseNueHiCVN->SetLineColor(kBlue-9);
110  multiverseNueHiCVN->SetFillColor(kBlue-9);
111  multiverseNueMerge = (TH1D*)inFileMultiverse->Get("rhc_totNDnue_tot");
112  multiverseNueMerge->SetDirectory(0);
113  multiverseNueMerge->SetLineColor(kBlue-9);
114  multiverseNueMerge->SetFillColor(kBlue-9);
115  multiverseNumu = (TH1D*)inFileMultiverse->Get("rhc_totNDnumu_tot");
116  multiverseNumu->SetDirectory(0);
117  multiverseNumu->SetLineColor(kBlue-9);
118  multiverseNumu->SetFillColor(kBlue-9);
119  inFileMultiverse->Close();
120 
121  // Binning
122  // Nue plots
123  int nBinXEnom = multiverseNueLoCVN->GetXaxis()->GetNbins();
124  if( nBinXEnom!=multiverseNueHiCVN->GetXaxis()->GetNbins() ) std::cout << "Nue plots have different x bins!" << std::endl;
125  if( nBinXEnom == 1 ) std::cout << "ONLY ONE BIN DETECTED!! This will produce nasty results..." << std::endl;
126  hMinXE = multiverseNueLoCVN->GetXaxis()->GetBinLowEdge(1);
127  int nBinXE = nBinXEnom+20;
128  hMaxXE = hMinXE + (nBinXE*(multiverseNueLoCVN->GetXaxis()->GetBinLowEdge(2) - multiverseNueLoCVN->GetXaxis()->GetBinLowEdge(1)));
129  hMaxYE = 1.1*std::max(multiverseNueLoCVN->GetMaximum(),multiverseNueHiCVN->GetMaximum());
130 
131  // Hi CVN
132  nBinXE2 = nBinXEnom*(float(multiverseNueHiCVN->GetMaximumBin())/float(multiverseNueLoCVN->GetMaximumBin())) + 18;
133  hMaxYE2 = hMaxYE;
134  hMinXE2 = multiverseNueHiCVN->GetXaxis()->GetBinLowEdge(1);
135  hMaxXE2 = hMinXE2 + (nBinXE2*(multiverseNueHiCVN->GetXaxis()->GetBinLowEdge(2) - multiverseNueHiCVN->GetXaxis()->GetBinLowEdge(1)));
136 
137  // Numu plots
138  hMaxYM = 1.1*multiverseNumu->GetMaximum();
139  nBinXM = 42;
140  hMinXM = multiverseNumu->GetXaxis()->GetBinLowEdge(1);
141  hMaxXM = hMinXM + (nBinXM*(multiverseNumu->GetXaxis()->GetBinLowEdge(2) - multiverseNumu->GetXaxis()->GetBinLowEdge(1)));
142  }
143  // ---------------------------
144 
145  TH2F *blankHistM = new TH2F("blankHistM",";Wrong-sign fraction in #nu_{#mu} selection;Fake experiments",nBinXM,hMinXM,hMaxXM,nBinYM,hMinYM,hMaxYM);
146  CenterTitles(blankHistM);
147 
148  // this axis has no meaning in the data-driven results, so just offset for visibility
149  float yVal[] = {150.,250.,350.};
150  float yErr[] = { 0., 0., 0.};
151  float yValNumu[] = {float(150.)*float(hMaxYM/hMaxYE),float(250.)*float(hMaxYM/hMaxYE),float(350.)*float(hMaxYM/hMaxYE)};
152 
153  // For using only the one Numu result
154  TPaveText *numu1Text = new TPaveText(0.12,0.195,0.375,0.285,"NB NDC");
155  TText *numu1TextL1 = numu1Text->AddText("Neutron Capture Rate");
156  numu1TextL1->SetTextAlign(12);
157  numu1Text->SetFillStyle(3155);
158 
159  // For the MC Universes
160  TPaveText *addedTextMC1;
161  TPaveText *addedTextMC2;
162  TPaveText *addedTextMC3;
163  if( useMultiverse ){
164  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
165  addedTextMC2 = new TPaveText(0.63, 0.585, 0.88, 0.67, "NB NDC");
166  addedTextMC3 = new TPaveText(0.63, 0.56, 0.875, 0.615, "NB NDC"); // xlow 0.625, ylow 0.54
167  TText *addedText1 = addedTextMC1->AddText("Flux and cross-");
168  TText *addedText2 = addedTextMC2->AddText("section systematic");
169  TText *addedText3 = addedTextMC3->AddText("fake experiments");
170  addedText1->SetTextColor(kBlue-7);
171  addedText2->SetTextColor(kBlue-7);
172  addedText3->SetTextColor(kBlue-7);
173  addedTextMC1->SetFillStyle(3155);
174  addedTextMC2->SetFillStyle(3155);
175  addedTextMC3->SetFillStyle(3155);
176  }
177 
178  // Numu Results: see notes listed at top for more info
179  // n capture study:
180  // Take scale factors from study and recast them as WS Fraction in numu selection
181  // Updated for 7.9e20 POT: 4+6,7d,8b
182  // ------------------------------------------------------------------------------------
183  // -----------------------
184  // --- n capture study ---
185  // -----------------------
186  float scaleWS = 1.05; // UPDATED FOR 2019!!
187  float errorWS = 0.12; // UPDATED FOR 2019!!
188  float scaleNC = 1.0;
189  float errorNC = 0.29;
190 
191  float WSCC = 105176; // all selected CC nubar events
192  float WSNC = 909.995; // all selected NC nubar events
193  float RSCC = 884586; // all selected CC nu events
194  float RSNC = 1737; // all selected NC nu events
195  // All selected numu data at the same POT scaling.
196  float numuData = 983130;
197 
198  // Few other parameters of interest for plots
199  float numuSelNue = 30.8344;
200  float numuSelNumu = 989731;
201  float numuSelNC = 2646.99;
202 
203  float numuFractionNue = numuSelNue / (numuSelNue + numuSelNumu + numuSelNC);
204  float numuFractionNumu = numuSelNumu / (numuSelNue + numuSelNumu + numuSelNC);
205  float numuFractionNC = numuSelNC / (numuSelNue + numuSelNumu + numuSelNC);
206 
207  // Calculate WS Fraction
208  float denomNumu = numuData;
209  float denomsq = denomNumu*denomNumu;
210  float xValNumu[] = { ((scaleWS*WSCC) + (scaleNC*WSNC)) / denomNumu };
211  // covariance 0 because not considering covariance in current handling with pull term
212  float xErrNumu[] = { sqrt( errorWS*errorWS*WSCC*WSCC/(numuData*numuData) + errorNC*errorNC*WSNC*WSNC/(numuData*numuData) ) };
213 
214  TGraphErrors *graphNumu = new TGraphErrors(1,xValNumu,yValNumu,xErrNumu,yErr);
215  graphNumu->SetMarkerStyle(24);
216  graphNumu->SetLineWidth(2);
217  // ------------------------------------------------------------------------------------
218 
219  // Pave Texts for composition fractions
220  TPaveText *typeCompText = new TPaveText(0.24,0.66,0.34,0.84,"NB NDC");
221  TText *typeCompTextL1 = typeCompText->AddText( "#nu_{e} + #bar{#nu}_{e}" );
222  typeCompTextL1->SetTextAlign(12);
223  TText *typeCompTextL2 = typeCompText->AddText( "#nu_{#mu} + #bar{#nu}_{#mu}" );
224  typeCompTextL2->SetTextAlign(12);
225  TText *typeCompTextL3 = typeCompText->AddText( "NC" );
226  typeCompTextL3->SetTextAlign(12);
227  typeCompText->SetFillStyle(3155);
228 
229  TPaveText *numuCompText = new TPaveText(0.14,0.66,0.24,0.84,"NB NDC"); // had used ylo = 0.57
230  TText *numuCompTextL1 = numuCompText->AddText( TString::Format("%.1f%%",100.*numuFractionNue) ); // try 1 decimal place
231  numuCompTextL1->SetTextAlign(32);
232  TText *numuCompTextL2 = numuCompText->AddText( TString::Format("%.1f%%",100.*numuFractionNumu) );
233  numuCompTextL2->SetTextAlign(32);
234  TText *numuCompTextL3 = numuCompText->AddText( TString::Format("%.1f%%",100.*numuFractionNC) );
235  numuCompTextL3->SetTextAlign(32);
236  numuCompText->SetFillStyle(3155);
237 
238  // Draw plots
239  TCanvas *cNumu = new TCanvas();
240  blankHistM->Draw("");
241  if( useMultiverse ){
242  multiverseNumu->Draw("same hist ][");
243  addedTextMC1->Draw();
244  addedTextMC2->Draw();
245  addedTextMC3->Draw();
246  }
247  graphNumu->Draw("same P");
248  numu1Text->Draw();
249  typeCompText->Draw();
250  numuCompText->Draw();
251  blankHistM->Draw("axissame");
252  Preliminary();
253  cNumu->Print("cNumu_Ana2019_WSFraction.pdf");
254 
255  std::cout << std::endl;
256  std::cout << std::endl;
257  std::cout << "===================================================" << std::endl;
258  std::cout << " Wrong-sign Fraction " << std::endl;
259  std::cout << "===================================================" << std::endl;
260  std::cout << " N Capture Numu WS Fraction " << xValNumu[0] << " +/- " << xErrNumu[0] << std::endl;
261  std::cout << "===================================================" << std::endl;
262  std::cout << std::endl;
263  std::cout << std::endl;
264 
265  std::cout << "Done!" << std::endl;
266 }
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 Preliminary()
OStream cout
Definition: OStream.cxx:6
void plotNumuWrongSign()
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