first_ana_proj.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Cuts/Cuts.h"
4 #include "CAFAna/Fit/Fit.h"
5 #include "CAFAna/Vars/FitVars.h"
11 
12 #include "TCanvas.h"
13 #include "TGraph.h"
14 #include "TLegend.h"
15 
16 #include <iostream>
17 
18 #include "Utilities/func/MathUtil.h"
19 
20 using namespace ana;
21 
22 const Var& pidVar = kLEM;
23 const TString pidName = "LEM";
24 const double cutVal = 0.4;
25 
27 {
28  const std::string fname = "/nfs/raid6/beam.mc.fhc.nonswap.loose.caf.root";
29  const std::string fnameSwap = "/nfs/raid6/beam.mc.fhc.swap.loose.caf.root";
30 
31  SpectrumLoader loader(fname);
32  SpectrumLoader loaderSwap(fnameSwap);
33 
34  // Make a calculator. This is the fastest variant
36 
37  calc.SetL(810);
38  calc.SetRho(2.75);
39  calc.SetDmsq21(7.6e-5);
40  calc.SetDmsq32(2.35e-3);
41  calc.SetTh12(asin(sqrt(.87))/2);
42  calc.SetTh13(asin(sqrt(.10))/2);
43  calc.SetTh23(TMath::Pi()/4);
44  calc.SetdCP(3*TMath::Pi()/2); // Biggest nue appearance
45 
46  // Set these parameters again this way to make sure we do it right
47  kFitSinSq2Theta13.SetValue(&calc, 0.1);
48  kFitSinSq2Theta23.SetValue(&calc, 1);
49 
50  Binning bins = Binning::Simple(1, 1, 2);
51 
52  const int kNumPreds = 51;
53  const double kCutStep = .02;
54  IPrediction* preds[kNumPreds];
55 
56  for(int i = 0; i < kNumPreds; ++i){
57  preds[i] = new PredictionNoExtrap(loader, loaderSwap,
58  pidName.Data(), bins,
59  kUnweighted, pidVar > i*kCutStep);
60  }
61 
62  IPrediction* predNominalCut = new PredictionNoExtrap(loader, loaderSwap,
63  pidName.Data(), bins,
64  kUnweighted, pidVar > cutVal);
65 
66  loader.Go();
67  loaderSwap.Go();
68 
69  TGraph* gSig = new TGraph;
70  TGraph* gBkg = new TGraph;
71  TGraph* gFOM = new TGraph;
72  TGraph* gSigma = new TGraph;
73 
74  TGraph* g2Sig = new TGraph;
75  TGraph* g3Sig = new TGraph;
76  TGraph* g4Sig = new TGraph;
77  TGraph* g5Sig = new TGraph;
78 
79  for(int i = 0; i < kNumPreds; ++i){
80  // What we expect to observe at those oscillation parameters
81 
82  const double nsig = preds[i]->PredictComponent(&calc,
85  Sign::kBoth).ToTH1(1e20)->Integral();
86 
87  const double nbkg = preds[i]->Predict(&calc).ToTH1(1e20)->Integral()-nsig;
88 
89  std::cout << nsig << " " << nbkg << std::endl;
90 
91  gSig->SetPoint(i, kCutStep*i, nsig);
92  gBkg->SetPoint(i, kCutStep*i, nbkg);
93 
94  gFOM->SetPoint(i, kCutStep*i, nsig ? nsig/sqrt(nsig+nbkg) : 0);
95  gSigma->SetPoint(i, kCutStep*i, nbkg ? nsig/sqrt(nbkg) : 0);
96 
97  double prob2 = 1;
98  double prob3 = 1;
99  double prob4 = 1;
100  double prob5 = 1;
101  for(int nobs = 0; ; ++nobs){
102  // What we would quote as our significance. Probability of seeing this
103  // many or more, given H0 of background-only.
104  double quote = 1;
105  for(int j = 0; j < nobs; ++j)
106  quote -= TMath::Poisson(j, nbkg);
107 
108  if(quote > TMath::Erfc(2/sqrt(2))) prob2 -= TMath::Poisson(nobs, nsig+nbkg);
109  if(quote > TMath::Erfc(3/sqrt(2))) prob3 -= TMath::Poisson(nobs, nsig+nbkg);
110  if(quote > TMath::Erfc(4/sqrt(2))) prob4 -= TMath::Poisson(nobs, nsig+nbkg);
111  if(quote > TMath::Erfc(5/sqrt(2))) prob5 -= TMath::Poisson(nobs, nsig+nbkg);
112  if(quote < TMath::Erfc(5/sqrt(2))) break;
113 
114  // LL formula
115  // double sigma = 2*(nbkg-nobs);
116  // if(nobs > 0) sigma += 2*nobs*log(nobs/nbkg);
117  // sigma = sqrt(sigma);
118 
119  // const double sigma = (nobs-nbkg)/sqrt(nbkg);
120 
121  /*
122  if(sigma < 2) prob2 -= TMath::Poisson(nobs, nsig+nbkg);
123  if(sigma < 3) prob3 -= TMath::Poisson(nobs, nsig+nbkg);
124  if(sigma < 4) prob4 -= TMath::Poisson(nobs, nsig+nbkg);
125  if(sigma < 5) prob5 -= TMath::Poisson(nobs, nsig+nbkg);
126  if(sigma > 5) break;
127  */
128  }
129  if(nsig == 0) prob2 = prob3 = prob4 = prob5 = 0;
130 
131  g2Sig->SetPoint(i, kCutStep*i, 100*prob2);
132  g3Sig->SetPoint(i, kCutStep*i, 100*prob3);
133  g4Sig->SetPoint(i, kCutStep*i, 100*prob4);
134  g5Sig->SetPoint(i, kCutStep*i, 100*prob5);
135  }
136 
137  (new TH2F("", ";"+pidName+" cut;Events", 100, 0, 1, 100, 0, 10))->Draw();
138 
139  gSig->SetLineWidth(2);
140  gSig->SetLineColor(kRed);
141  gSig->Draw("l same");
142  gBkg->SetLineWidth(2);
143  gBkg->SetLineColor(kBlue);
144  gBkg->Draw("l same");
145 
146  gFOM->SetLineWidth(2);
147  gFOM->SetLineColor(kGray);
148  gFOM->Draw("l same");
149  gSigma->SetLineWidth(2);
150  gSigma->Draw("l same");
151 
152  TLegend* leg = new TLegend(.6, .65, .85, .85);
153  leg->SetFillStyle(0);
154  leg->AddEntry(gSig, "Signal", "l");
155  leg->AddEntry(gBkg, "Background", "l");
156  leg->AddEntry(gFOM, "s/#sqrt{s+b}", "l");
157  leg->AddEntry(gSigma, "s/#sqrt{b}", "l");
158  leg->Draw();
159 
160  gPad->Print("plots/first_ana_counts.png");
161  gPad->Print("plots/first_ana_counts.eps");
162 
163  new TCanvas;
164 
165  (new TH2F("", ";"+pidName+" cut;Odds of rejecting #theta_{13}=0 (%)", 100, 0, 1, 100, 0, 100))->Draw();
166  g2Sig->SetLineColor(kBlue);
167  g3Sig->SetLineColor(kBlack);
168  g4Sig->SetLineColor(kGreen+2);
169  g5Sig->SetLineColor(kRed);
170 
171  g2Sig->SetLineWidth(2);
172  g2Sig->Draw("l same");
173  g3Sig->SetLineWidth(2);
174  g3Sig->Draw("l same");
175  g4Sig->SetLineWidth(2);
176  g4Sig->Draw("l same");
177  g5Sig->SetLineWidth(2);
178  g5Sig->Draw("l same");
179 
180  leg = new TLegend(.15, .8, .85, .85);
181  leg->SetNColumns(4);
182  leg->SetFillStyle(0);
183  leg->AddEntry(g2Sig, "2#sigma", "l");
184  leg->AddEntry(g3Sig, "3#sigma", "l");
185  leg->AddEntry(g4Sig, "4#sigma", "l");
186  leg->AddEntry(g5Sig, "5#sigma", "l");
187  leg->Draw();
188 
189  gPad->Print("plots/first_ana_odds.png");
190  gPad->Print("plots/first_ana_odds.eps");
191 
192  new TCanvas;
193 
194  const double nnom = predNominalCut->Predict(&calc).ToTH1(1e20)->Integral();
195 
196  TH1* hPoi = new TH1F("", ";Number of candidates;Probability (%)", 10, 0, 10);
197  for(int i = 0; i < 10; ++i){
198  hPoi->Fill(i, 100*TMath::Poisson(i, nnom));
199  }
200 
201  hPoi->Draw();
202 
203  gPad->Print("plots/poi.png");
204  gPad->Print("plots/poi.eps");
205 
206  for(int nobs = 0; nobs < 10; ++nobs){
207  new TCanvas;
208 
209  Spectrum data("Count", 1e20, 0, bins);
210  data.Fill(1, nobs);
211 
212  SingleSampleExperiment expt(predNominalCut, data);
213 
214  // Normal hierarchy
215  calc.SetDmsq32(+2.35e-3);
216 
217  // The log-likelihood surface
218  FrequentistSurface surf(&expt, &calc,
219  &kFitSinSq2Theta13, 25, 0, .5,
220  &kFitDeltaInPiUnits, 25, 0, 2);
221  surf.SetTitle(TString::Format("%d candidate%s %d%% chance", nobs, nobs != 1 ? "s" : "", int(100*TMath::Poisson(nobs, nnom)+.5)));
222 
223  surf.DrawContour(Gaussian68Percent1D(surf), 7, kRed); // dashed
224  surf.DrawContour(Gaussian90Percent1D(surf), 1, kRed); // solid
225 
226  // Inverted hierarchy
227  calc.SetDmsq32(-2.35e-3);
228 
229  FrequentistSurface surfInv(&expt, &calc,
230  &kFitSinSq2Theta13, 25, 0, .5,
231  &kFitDeltaInPiUnits, 25, 0, 2);
232 
233  surfInv.DrawContour(Gaussian68Percent1D(surf), 7, kBlue); // dashed
234  surfInv.DrawContour(Gaussian90Percent1D(surf), 1, kBlue); // solid
235 
236  // Daya Bay
237  TGraph* gdb1 = new TGraph;
238  gdb1->SetPoint(0, .084-.005, 0);
239  gdb1->SetPoint(1, .084-.005, 2);
240  gdb1->SetLineWidth(2);
241  gdb1->Draw("l same");
242 
243  TGraph* gdb2 = new TGraph;
244  gdb2->SetPoint(0, .084+.005, 0);
245  gdb2->SetPoint(1, .084+.005, 2);
246  gdb2->SetLineWidth(2);
247  gdb2->Draw("l same");
248 
249  gPad->Print(TString::Format("plots/cont_%d.png", nobs));
250  gPad->Print(TString::Format("plots/cont_%d.eps", nobs));
251  } // end for nobs
252 }
TH2 * Gaussian90Percent1D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 1D in gaussian approximation.
const Var kLEM
PID
Definition: Vars.cxx:24
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
const TString pidName
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
void SetDmsq21(const T &dmsq21) override
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
const FitSinSq2Theta23 kFitSinSq2Theta23
string fnameSwap
Definition: demo4.py:6
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetTitle(const char *str)
Definition: ISurface.cxx:192
void Fill(double x, double w=1)
Definition: Spectrum.cxx:774
const double cutVal
void SetTh23(const T &th23) override
const Var & pidVar
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
void SetDmsq32(const T &dmsq32) override
const XML_Char const XML_Char * data
Definition: expat.h:268
Log-likelihood scan across two parameters.
void SetdCP(const T &dCP) override
Charged-current interactions.
Definition: IPrediction.h:39
expt
Definition: demo5.py:34
std::string quote(const T &val)
Definition: utility.hpp:156
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:45
void SetL(double L) override
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
void SetTh13(const T &th13) override
loader
Definition: demo0.py:10
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:46
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void SetTh12(const T &th12) override
recTree Draw("energy.numu.trkccE>>precosmics","fdcuts&&preshutdown")
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void SetRho(double rho) override
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Prediction that just uses FD MC, with no extrapolation.
void first_ana_proj()
void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const override
Definition: FitVars.cxx:103
Float_t e
Definition: plot.C:35
loaderSwap
Definition: demo4.py:9
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:100
TH2 * Gaussian68Percent1D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 1D in gaussian approximation.
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35