T2KToyExperiment.cxx
Go to the documentation of this file.
2 
4 
5 #include "OscLib/OscCalcPMNS.h"
6 #include "Utilities/func/MathUtil.h"
7 
8 #include "TFile.h"
9 #include "TH1.h"
10 
11 #include <cassert>
12 
13 namespace ana
14 {
15  /// Baseline from JPARC to SK in km
16  const int kT2KBaseline = 295;
17 
18  //----------------------------------------------------------------------
20  : fPOT(0), fPOTAnti(0)
21  {
22  // The original version of this code was written by Ryan. Ask him about the
23  // physics assumptions in it.
24 
25  // Scale everything from Neutrino 2013 results
26  fBasePOT = 2.556e20;
27  fBaseSig = 6.60;
28  fBaseBkg = 2.47;
29 
30  fTotSyst = 0.08; // 0.103 at Nu2012
31 
32  // using 3.8x reduction (15% from flux: 11/11/2012 IPMU talk, Nakaya;
33  // from an old 2-deg off-axis simulation. xsec reduction from those
34  // usual xsec plots from Sam Z., taken at 0.8 GeV)
35  // For BG: wrong-sign goes up by 3.8
36  // For Signal: nu to nubar adjustment to fixed pot obtained by hacking
37  // GetSignalScaleFactorT2K() to compare nubar to nu. Result: 0.797151
38  fBaseSigAnti = (fBaseSig / 3.8) * 0.797151;
39  fBaseBkgAnti = 2.34 / 3.8 + 0.13 * 3.8;
40 
42  osc->SetL(kT2KBaseline);
43  osc->SetRho(3.2);
44  osc->SetDmsq21(7.6e-5);
45  osc->SetDmsq32(2.4e-3);
46  osc->SetTh12(.601);
47  osc->SetTh13(asin(sqrt(.1))/2);
48  osc->SetTh23(M_PI/4);
49  osc->SetdCP(0);
50  fOscBase = osc;
51 
52  LoadSpectra();
53 
54  // By default, fit to Neutrino results
55  SetTruthParams(osc);
56  }
57 
58  //----------------------------------------------------------------------
60  {
61  // It seems very likely people are going to call us with NOvA
62  // baseline. Ignore them.
63  const double tmpL = osc->GetL();
64  osc->SetL(kT2KBaseline);
65 
68 
69  osc->SetL(tmpL);
70  }
71 
72  //----------------------------------------------------------------------
74  {
75  delete fOscBase;
76  }
77 
78  //----------------------------------------------------------------------
80  {
81  TFile flux((FindCAFAnaDir()+"/data/beam/flux_t2k.root").c_str());
82  assert(!flux.IsZombie());
83  fSpect = (TH1*)flux.Get("flux_t2k")->Clone();
84  assert(fSpect);
85  fSpect->SetDirectory(0);
86  flux.Close();
87 
88  // cross section and osc weights
89  for(int i = 1; i <= fSpect->GetNbinsX(); ++i){
90  const double E = fSpect->GetXaxis()->GetBinCenter(i);
91  double val = fSpect->GetBinContent(i);
92  val *= E; // take QE xsec proportional to E
93 
94  // kill tiny bins at low E
95  if(E < 0.4) val = 0;
96 
97  // kill bins above the E_reco cut at 1.25
98  if(E > 1.25) val = 0;
99 
100  fSpect->SetBinContent(i, val);
101  }
102  // Normalize
103  fSpect->Scale(1./fSpect->Integral());
104  }
105 
106  //----------------------------------------------------------------------
107  TH1* T2KToyExperiment::GetReweightedSignal(TH1* h, TH1* ret, bool anti,
108  osc::IOscCalc* calc) const
109  {
110  if(!ret || ret->GetNbinsX() != h->GetNbinsX())
111  ret = (TH1*)h->Clone();
112  // else, assume the histo has been provided
113 
114  // Average over the bin since the ratio can be quite spikey at any given E
115  const int kNumSteps = 30;
116 
117  const int from = anti ? -14 : +14;
118  const int to = anti ? -12 : +12;
119  for(int i = 1; i <= h->GetNbinsX(); ++i){
120  if(h->GetBinContent(i) == 0){
121  ret->SetBinContent(i, 0);
122  continue;
123  }
124 
125  double pTarget = 0;
126  double pBase = 0;
127  for(double E = h->GetXaxis()->GetBinLowEdge(i);
128  E < h->GetXaxis()->GetBinUpEdge(i);
129  E += h->GetXaxis()->GetBinWidth(i)/kNumSteps){
130  pTarget += calc->P(from, to, E);
131  pBase += fOscBase->P(from, to, E);
132  }
133 
134  ret->SetBinContent(i, h->GetBinContent(i)*pTarget/pBase);
135  }
136  return ret;
137  }
138 
139  //----------------------------------------------------------------------
141  osc::IOscCalc* calc) const
142  {
143  static TH1* hW = 0;
144  hW = GetReweightedSignal(fSpect, hW, anti, calc);
145  return hW->Integral()/ fSpect->Integral();
146  }
147 
148  //----------------------------------------------------------------------
149  double T2KToyExperiment::LogL(double data, double pred, double syst) const
150  {
151  // Scale as if it were a simple chi^2 to get some systematics in there
152  const double syst_factor = pred ? pred/(pred + util::sqr(pred*syst)) : 1;
153 
154  return 2 * (pred-data+((data && pred) ? (data*log(data/pred)) : 0)) * syst_factor;
155  }
156 
157  //----------------------------------------------------------------------
159  const SystShifts& /*shift*/) const
160  {
161  // It seems very likely people are going to call us with NOvA
162  // baseline. Ignore them.
163  const double tmpL = calc->GetL();
164  calc->SetL(kT2KBaseline);
165 
166  double chi = 0;
167 
168  if(fPOT){
169  // Scale for POTs
170  const double scaledSig = (fPOT/fBasePOT) * fBaseSig;
171  const double scaledBkg = (fPOT/fBasePOT) * fBaseBkg;
172 
173  // And then for two different sets of osc parameters
174  const double data = scaledSig * fSigFactorTrue + scaledBkg;
175  const double pred = scaledSig * GetSignalScaleFactor(false, calc) + scaledBkg;
176 
177  chi += LogL(data, pred, fTotSyst);
178  }
179  if(fPOTAnti){
180  const double scaledSigAnti = (fPOTAnti/fBasePOT) * fBaseSigAnti;
181  const double scaledBkgAnti = (fPOTAnti/fBasePOT) * fBaseBkgAnti;
182 
183  const double dataAnti = scaledSigAnti * fSigFactorAntiTrue + scaledBkgAnti;
184  const double predAnti = scaledSigAnti * GetSignalScaleFactor(true, calc) + scaledBkgAnti;
185 
186  chi += LogL(dataAnti, predAnti, fTotSyst);
187  }
188 
189  calc->SetL(tmpL);
190 
191  return chi;
192  }
193 } // end namespace
virtual void SetL(double L)=0
void SetDmsq21(const T &dmsq21) override
Definition: OscCalcPMNS.h:35
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh12(const T &th12) override
Definition: OscCalcPMNS.h:37
void SetTh23(const T &th23) override
Definition: OscCalcPMNS.h:39
const SpillVar fPOT([](const caf::SRSpillProxy *spill){return spill->spillpot;})
Adapt the PMNS calculator to standard interface.
Definition: StanTypedefs.h:28
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SetL(double L) override
Definition: OscCalcPMNS.h:33
T sqrt(T number)
Definition: d0nt_math.hpp:156
const int kT2KBaseline
Baseline from JPARC to SK in km.
virtual double GetL() const
Definition: IOscCalc.h:83
osc::OscCalcDumb calc
std::string FindCAFAnaDir()
Definition: Utilities.cxx:429
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Loaders::FluxType flux
#define M_PI
Definition: SbMath.h:34
double LogL(double data, double pred, double syst) const
const XML_Char const XML_Char * data
Definition: expat.h:268
double GetSignalScaleFactor(bool anti, osc::IOscCalc *calc) const
osc::IOscCalc * fOscBase
Float_t E
Definition: plot.C:20
void SetdCP(const T &dCP) override
Definition: OscCalcPMNS.h:40
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcPMNS.h:36
void SetTh13(const T &th13) override
Definition: OscCalcPMNS.h:38
TH1 * fSpect
Unoscillated, normalized.
void SetRho(double rho) override
Definition: OscCalcPMNS.h:34
Oscillation probability calculators.
Definition: Calcs.h:5
virtual T P(int flavBefore, int flavAfter, double E)=0
E in GeV; flavors as PDG codes (so, neg==>antinu)
_OscCalcPMNS< double > OscCalcPMNS
Definition: OscCalcPMNS.h:57
assert(nhit_max >=nhit_nbins)
void SetTruthParams(osc::IOscCalcAdjustable *osc)
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &shift=SystShifts::Nominal()) const override
Float_t e
Definition: plot.C:35
TH1 * GetReweightedSignal(TH1 *h, TH1 *ret, bool anti, osc::IOscCalc *calc) const
T asin(T number)
Definition: d0nt_math.hpp:60