joint_fit_2017_tools.h
Go to the documentation of this file.
1 #pragma once
2 
15 #include "CAFAna/Systs/XSecSysts.h"
16 #include "CAFAna/Vars/FitVars.h"
17 using namespace ana;
18 
19 #include "OscLib/IOscCalc.h"
20 
21 #include "TFile.h"
22 
23 
26  bool corrSysts,
27  bool GetFromUPS=false){
28 
30  if(decomp == "noextrap") extrap = kNoExtrap;
31  if(decomp == "prop") extrap = kProportional;
32  if(decomp == "combo") extrap = kCombo;
33 
34  PredictionSystNue2017 * pred_fid;
35  // Interactive running
36  if (!GetFromUPS){
37  if(corrSysts) pred_fid = new PredictionSystNue2017 (extrap, calc);
38  else pred_fid = new PredictionSystNue2017 (extrap, calc,{});
39  }
40  // Grid running for FC
41  else{
42  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
43  std::string name = "/pred_nue_reduced_v9.root";
44  if(corrSysts)
45  pred_fid = new PredictionSystNue2017(extrap, calc, kNue2017Systs,
46  upsname+name);
47  else
48  pred_fid = new PredictionSystNue2017(extrap, calc, {},
49  upsname+name);
50  }
51 
52  std::cout << "Nue prediction (no rock)"
53  << pred_fid->Predict(calc).Integral(kAna2017POT) << std::endl;
54  // Load rock
55 
56  double nonsScale = 1./10.45;
57  double swapScale = 1./12.91;
58 
59  std::string dirName = "/nova/ana/nu_e_ana/Ana2017/Predictions/FDRock/";
60  std::string fName2 = "fdrock_nue2017.root";
61  std::string rName = "pred_FDRock";
62 
63  // Point to the rock file in the ups product for FC'ing
64  if (GetFromUPS){
65  dirName = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
66  fName2 = "/fdrock_nue2017.root";
67  rName = "pred_FDRock";
68  }
69 
70  auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).release();
71 
72  std::cout << "Nue prediction (rock only - not scaled) "
73  << rock->Predict(calc).Integral(kAna2017POT) << std::endl;
74 
75  auto nue_pred = new PredictionAddRock (pred_fid, rock, nonsScale, swapScale);
76 
77  std::cout << "Nue prediction (rock added) "
78  << nue_pred->Predict(calc).Integral(kAna2017POT) << std::endl;
79  return nue_pred;
80 }
81 
82 //////////////////////////////////////////////////////////////////////
83 
84 std::pair <TH1D*, double > GetNueCosmics2017 (bool GetFromUPS=false){
86  if (GetFromUPS){
87  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
88  std::string fname = "/nue2017_cospred.root";
89  std::string name = "cos";
90  outtime = LoadFromFile<Spectrum>(upsname+fname,name).release();
91  }
92  else{
93  outtime = LoadFromFile<Spectrum>
94  ("/nova/ana/nu_e_ana/Ana2017/Predictions/cosmic/cosmic_spectra_prediction.root",
95  "All/cosm_merged_4bin" ).release();
96  }
97 
98  double outN = outtime->Integral(outtime->Livetime(), 0, kLivetime);
99 
100  std::cout << "\nAdding " << outtime->Integral(outtime->Livetime(), 0, kLivetime)
101  << " cosmics from out-of-time." << std::endl;
102  std::cout << "Which scale to " << outtime->Integral(kAna2017Livetime, 0,
103  kLivetime)
104  << " in the spill window." << std::endl;
105 
106  std::cout << "uncertainty " << 1/sqrt(outN) << std::endl;
107 
108  TH1 *h = outtime->ToTH1(kAna2017Livetime, kBlack, kSolid, kLivetime);
109  assert( h->Integral() > 4 && "Nue cosmic prediction < 4. Something is wrong");
110 
111  return {outtime->ToTH1(kAna2017Livetime, kBlack, kSolid, kLivetime), 1/sqrt(outN)};
112 }
113 
114 //////////////////////////////////////////////////////////////////////
115 
117  Spectrum* intime = LoadFromFile<Spectrum>
118  ("/nova/ana/nu_e_ana/Ana2017/data/data_spectra_fd.root",
119  "All/spect_merged_4bin").release();
120  std::cout << "Loaded data spectrum:\n"
121  << "POT " << intime->POT()
122  << " (expected " << kAna2017POT << ")\n"
123  << "Livetime " << intime->Livetime()
124  << " (expected " << kAna2017Livetime << ")\n"
125  << "Integral " << intime->Integral(intime->POT())
126  << std::endl;
127 
128  return intime;
129 }
130 
131 //////////////////////////////////////////////////////////////////////
132 
133 std::vector <const IPrediction*> GetNumuPredictions2017(const int nq = 4,
134  bool useSysts = true,
135  bool GetFromUPS=false)
136 {
137  auto calc = DefaultOscCalc();
138 
139  std::vector <const IPrediction*> numu_preds;
140 
141  std::cout << "\nLoading Numu Predictions\n" ;
142  for (int quant = 0; quant < nq; ++quant){
143  if (GetFromUPS){
144  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
145  std::string fname = "/pred_numu_reduced_v5.root";
146  if(!useSysts)
147  numu_preds.push_back(new PredictionSystNumu2017(kNuMu, calc,
148  quant+1,{},
149  upsname+fname));
150  else
151  numu_preds.push_back(new PredictionSystNumu2017(kNuMu, calc,
152  quant+1,kNumu2017Systs,
153  upsname+fname));
154  }
155  else{
156  if(!useSysts)
157  numu_preds.push_back(new PredictionSystNumu2017(kNuMu, calc,
158  quant+1,{}));
159  else
160  numu_preds.push_back(new PredictionSystNumu2017(kNuMu, calc, quant+1));
161  }
162  }
163  return numu_preds;
164 
165 }
166 
167 //////////////////////////////////////////////////////////////////////
168 
169 std::vector <std::pair <TH1D*, double> > GetNumuCosmics2017(const int nq = 4,
170  bool GetFromUPS=false)
171 {
172  //cosmics
173  //the versionn in cafana has float not double and it's a pain
174  TFile *fcosm;
175  if (GetFromUPS){
176  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
177  fcosm = TFile::Open((upsname+"/cos_bkg_hists_v2.root").c_str());
178  }
179  else{
180  std::string dir = "/nova/ana/nu_e_ana/Ana2017/NumuTemp";
181  fcosm = TFile::Open((dir+"/cos_bkg_hists_v2.root").c_str());
182  }
183 
184  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
185  std::vector <std::pair <TH1D*, double > > numu_cosmics;
186 
187  for (int i = 1; i <=nq; ++i) {
188  auto h = (TH1D*) fcosm->Get((TString)"q"+std::to_string(i)+"all")
189  ->Clone(UniqueName().c_str());
190  auto hs = (TH1D*) fcosm->Get((TString)"q"+std::to_string(i)+"allsyst");
191  numu_cosmics.push_back({h,hs->GetBinContent(1)});
192  std::cout << "Cosmics all periods quantile " << i << " "
193  << h->Integral()
194  << " scale error "
195  << hs->GetBinContent(1)
196  << "\n";
197  }//end quantiles
198  return numu_cosmics;
199 }
200 //////////////////////////////////////////////////////////////////////
201 
202 std::vector <Spectrum * > GetNumuData2017(const int nq = 4){
203  TString dir = "/pnfs/nova/persistent/analysis/numu/Ana2017/FD_data/";
204  std::vector <Spectrum *> numu_data;
205 
206  for (int i = 1; i <=nq; ++i) {
207 
208  auto f= new TFile(dir + "fd_data_Quant" + std::to_string(i) + ".root");
209  if(f->IsZombie()) {std::cerr<< "bad data\n"; exit(1);}
210  numu_data.push_back(Spectrum::LoadFrom(f, "fd_data")).release();
211 
212  }//end quantiles
213  return numu_data;
214 }
215 
216 //////////////////////////////////////////////////////////////////////
217 
219  const double pot, TH1D* cosmics = 0)
220 {
221  auto temp = pred->Predict(calc).FakeData(pot);
222  if(cosmics) temp += Spectrum(cosmics, pot,0);
223  return new Spectrum(temp);
224 }
225 
226 //////////////////////////////////////////////////////////////////////
227 
229  double dmsq32 = -5, double dCP_Pi = -5)
230 {
231  if(dCP_Pi > -1) kFitDeltaInPiUnits.SetValue(calc, dCP_Pi);
232  if(ssth23 > -1) kFitSinSqTheta23.SetValue(calc, ssth23);
233  if(dmsq32 > -1) kFitDmSq32.SetValue(calc, dmsq32);
234 }
235 
236 //////////////////////////////////////////////////////////////////////
237 
238 std::vector<const ISyst*> GetJointFitSystematicList(bool corrSysts, bool nueExclusive = false, bool numuExclusive=false, bool intersection = true){
239  std::vector<const ISyst*> ret ={} ;
240  //TODO this is broken
241  if(corrSysts){
242  if(! intersection){
243  if(nueExclusive){
244  ret.push_back(&kRadCorrNue);
245  ret.push_back(&kRadCorrNuebar);
246  ret.push_back(&k2ndClassCurrs);
248  }
249  if(numuExclusive){
251  }
252  }
253  else{
254  if(nueExclusive) ret = getAllAna2017Systs(kNueAna2017);
255  if(numuExclusive) ret = getAllAna2017Systs(kNumuAna2017);
256  if(nueExclusive&& numuExclusive) ret = getAllAna2017Systs(kJointAna2017);
257  }
258  }
259  std::cout << "Return list of systematics: " ;
260  for(auto & s:ret) std::cout << s->ShortName() << ", ";
261  std::cout << std::endl;
262  return ret;
263 }
264 std::vector<std::pair<const ISyst*,const ISyst*> > GetCorrelations (bool isNue){
265  std::vector<std::pair<const ISyst*,const ISyst*> > ret ={};
266 
267  //to first approximation if they are not the same they are not correlated.
268  //might have to be more complicated later
269  auto badsyst = GetJointFitSystematicList(true, !isNue, isNue, false);
270  for (auto & s:badsyst) ret.push_back ({s,NULL});
271 
272  return ret;
273 }
274 
275 //----------------------------------------------------------------------
276 std::map<std::string, double> SumSysts2017 (std::map<std::string, double> systbars)
277 {
278  std::map<std::string, double> summedbars;
279  double summedsyst[kNumSystTypes]={ 0 };
280 
281  for(auto it: systbars){
282  if( it.first == "Bkg. Norm." || it.first == "Sig. Norm." || it.first == "Sig+Bkg Norm." )
283  {
284  summedsyst[kNormSys] += it.second*it.second;
285  continue;
286  }
287 
288  if( it.first == "AbsCalib" ||
289  it.first == "CalibShape" ||
290  it.first == "RelCalib" ||
291  it.first == "Absolute Muon Energy Scale 2017" ||
292  it.first == "Rel Muon Energy Scale 2017" )
293  {
294  summedsyst[kCalibSys] += it.second*it.second;
295  continue;
296  }
297 
298  if( it.first == "DIS vnCC1pi" ||
299  it.first == "MEC q_{0} shape" ||
300  it.first == "MaCCQE" ||
301  it.first == "MaCCRES" ||
302  it.first == "MaNCRES" ||
303  it.first == "MvCCRES" ||
304  it.first == "RPA shape: enh" ||
305  it.first == "RPA shape: resonance events" ||
306  it.first == "Radiative corrections for #bar{#nu_{e}}" ||
307  it.first == "Radiative corrections for #nu_{e}" ||
308  it.first == "Second class currents" ||
309  it.first == "Summed small XSec systs Joint2017" ||
310  it.first == "Extrapolation Bkg" ||
311  it.first == "Extrapolation ND to FD Kinematics Signal" )
312  {
313  summedsyst[kGenieSys] += it.second*it.second;
314  continue;
315  }
316 
317  if( it.first == "Cherenkov" || it.first == "Lightlevel" )
318  {
319  summedsyst[kBirksSys] += it.second*it.second;
320  continue;
321  }
322 
323  if( it.first == "PPFX Universe 00" ||
324  it.first == "PPFX Universe 01" ||
325  it.first == "PPFX Universe 02" ||
326  it.first == "PPFX Universe 03" ||
327  it.first == "PPFX Universe 04" ||
328  it.first == "Combined Beam Transport Systematics" )
329  {
330  summedsyst[kBeamSys] += it.second*it.second;
331  continue;
332  }
333 
334  std::cout<<"Systematic "<<it.first<<" does not belong."<<std::endl;
335  }
336 
337  for ( int i = 0; i < kNumSystTypes; ++i){
338  summedbars[(std::string)systTypes[i]]=sqrt(summedsyst[i]);
339  }
340 
341  return summedbars;
342 }
343 
344 
345 // For slices plots, plot by true octant, and make separate plots of each
347 {
348 public:
349  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
350  return util::sqr(sin(osc->GetTh23()));};
351  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
352  osc->SetTh23(asin(sqrt(Clamp(val))));};
353  virtual std::string ShortName() const {return "ssth23LO";}
354  virtual std::string LatexName() const {return "sin^{2}#theta_{23} LO";}
355  virtual double LowLimit() const {return 0;}
356  virtual double HighLimit() const {return 0.5;}
357 };
358 
360 
362 {
363 public:
364  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
365  return util::sqr(sin(osc->GetTh23()));};
366  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
367  osc->SetTh23(asin(sqrt(Clamp(val))));};
368  virtual std::string ShortName() const {return "ssth23UO";}
369  virtual std::string LatexName() const {return "sin^{2}#theta_{23} UO";}
370  virtual double LowLimit() const {return 0.5;}
371  virtual double HighLimit() const {return 1;}
372 };
374 
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2017(const int nq=4, bool GetFromUPS=false)
const XML_Char * name
Definition: expat.h:151
TCut intime("tave>=217.0 && tave<=227.8")
const NOvARwgtSyst k2ndClassCurrs("2ndclasscurr","Second class currents", novarwgt::kSimpleSecondClassCurrentsSystKnob)
Second-class current syst. See documentation in NOvARwgt.
Definition: XSecSysts.h:71
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double ssth23
set< int >::iterator it
void AddJointAna2017OtherSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana)
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:209
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
std::map< std::string, double > SumSysts2017(std::map< std::string, double > systbars)
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
virtual double HighLimit() const
OStream cerr
Definition: OStream.cxx:7
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
virtual T GetTh23() const
Definition: IOscCalc.h:89
const double kAna2017POT
Definition: Exposures.h:174
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:333
virtual double HighLimit() const
virtual double LowLimit() const
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
const std::vector< const DummyNue2017Syst * > kNue2017Systs
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const NOvARwgtSyst kRadCorrNuebar("radcorrnuebar","Radiative corrections for #bar{#nu}_{e}", novarwgt::kSimpleRadiativeCorrNuebarXsecSystKnob)
Radiative corrections syst (nuebars). See documentation in NOvARwgt.
Definition: XSecSysts.h:67
const std::vector< const DummyNumu2017Syst * > kNumu2017Systs
const XML_Char * s
Definition: expat.h:262
Loads shifted spectra from files.
std::pair< TH1D *, double > GetNueCosmics2017(bool GetFromUPS=false)
const double kAna2017Livetime
Definition: Exposures.h:200
std::vector< const IPrediction * > GetNumuPredictions2017(const int nq=4, bool useSysts=true, bool GetFromUPS=false)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
std::vector< Spectrum * > GetNumuData2017(const int nq=4)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
Loads shifted spectra from files.
Spectrum * GetNueData2017()
#define pot
std::string systTypes[kNumSystTypes]
T Clamp(T val) const
Definition: IFitVar.h:60
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
double POT() const
Definition: Spectrum.h:231
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
virtual std::string LatexName() const
Spectrum Predict(osc::IOscCalc *calc) const override
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
virtual double LowLimit() const
double dmsq32
TDirectory * dir
Definition: macro.C:5
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
exit(0)
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
assert(nhit_max >=nhit_nbins)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
virtual std::string ShortName() const
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
const NOvARwgtSyst kRadCorrNue("radcorrnue","Radiative corrections for #nu_{e}", novarwgt::kSimpleRadiativeCorrNueXsecSystKnob)
Radiative corrections syst (nues). See documentation in NOvARwgt.
Definition: XSecSysts.h:64
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void rock(std::string suffix="full")
Definition: rock.C:28
virtual std::string LatexName() const
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
T asin(T number)
Definition: d0nt_math.hpp:60
virtual std::string ShortName() const