joint_fit_make_experiments.C
Go to the documentation of this file.
1 //Creates a file with nue and numu SingleSampleExperiments
2 //Todo: these would be better as functions,
3 //optionally taking a calc for fake data, no cosmics, etc
4 //Todo: all POT from exposures.h
5 
12 #include "CAFAna/Systs/NumuSysts.h"
14 using namespace ana;
15 
17 
18 #include "TFile.h"
19 
20 
22  std::string decomp = "prop")
23 {
24  //////////////////////////////////////////////////
25  // Nue experiment
26  //////////////////////////////////////////////////
27 
29  std::string pidStr = "cvn";
30  if(PIDStr == "LID"){pid = kSystLID; pidStr = "lid";}
31  if(PIDStr == "LEM"){pid = kSystLEM; pidStr = "lem";}
32 
34  if(decomp == "noextrap") extrap = kNoExtrap;
35  if(decomp == "combo") extrap = kCombo;
36 
38  // pred.DebugPlots(DefaultOscCalc(), "debug_"+PIDStr+"_%s.eps");
39 
40  Spectrum* intime = LoadFromFile<Spectrum>("/nova/ana/nu_e_ana/SecondAna/NueSA_numi_data.root",
41  "spec_"+pidStr+"_"+pidStr+"2d_numi").release();
42  Spectrum* outtime = LoadFromFile<Spectrum>("/nova/ana/nu_e_ana/SecondAna/NueSA_numi_data.root",
43  "spec_"+pidStr+"_"+pidStr+"2d_cosmic").release();
44  double outN = outtime->Integral(outtime->Livetime(), 0, kLivetime);
45  double innN = outtime->Integral(intime->Livetime(), 0,
46  kLivetime);
47  std::cout << "Adding " << outtime->Integral(outtime->Livetime(), 0, kLivetime)
48  << " cosmics from out-of-time." << std::endl;
49  std::cout << "Which scale to " << outtime->Integral(intime->Livetime(), 0,
50  kLivetime)
51  << " in the spill window." << std::endl;
52  std::cout << "uncertainty " << innN/sqrt(outN) << std::endl;
53 
54  const SingleSampleExperiment exptNue(&pred, *intime, *outtime, 1/sqrt(outN));
55 
56  //////////////////////////////////////////////////
57  // Numu experiment
58  //////////////////////////////////////////////////
59 
60  // from numu/SA/plot_sa_results.C
61  std::cout << "\nLoading Numu experiment\n\n";
62 
63  TFile inputdata((FindCAFAnaDir()+"/numu/SA/fd_spec.root").c_str(), "READ");
64  // load in real data spectra. Precalculate since tons of FD files, takes forever until concats. I made the spectra with the grid and used hadd_cafa for the file
65  std::unique_ptr<Spectrum> FDSpectrumP1 = Spectrum::LoadFrom((TDirectory*)inputdata.Get("FDSpecP1"));
66  std::unique_ptr<Spectrum> FDSpectrumP2 = Spectrum::LoadFrom((TDirectory*)inputdata.Get("FDSpecP2"));
67  std::unique_ptr<Spectrum> FDSpectrumP3b = Spectrum::LoadFrom((TDirectory*)inputdata.Get("FDSpecP3b"));
68  std::unique_ptr<Spectrum> FDSpectrumP3c = Spectrum::LoadFrom((TDirectory*)inputdata.Get("FDSpecP3c"));
69  std::unique_ptr<Spectrum> FDSpectrumP3d = Spectrum::LoadFrom((TDirectory*)inputdata.Get("FDSpecP3d"));
70 
71  float potp1 = FDSpectrumP1->POT();
72  float potp2 = FDSpectrumP2->POT();
73  float potp3b = FDSpectrumP3b->POT();
74  float potp3c = FDSpectrumP3c->POT();
75  float potp3d = FDSpectrumP3d->POT();
76  float potp3 = potp3b + potp3c + potp3d;
77  float pot = potp1 + potp2 + potp3;
78 
79  // don't trust spectrum adding to do the right thing with POT
80  TH1 *hSpectrumP1 = FDSpectrumP1->ToTH1(potp1);
81  TH1 *hSpectrumP2 = FDSpectrumP2->ToTH1(potp2);
82  TH1 *hSpectrumP3b = FDSpectrumP3b->ToTH1(potp3b);
83  TH1 *hSpectrumP3c = FDSpectrumP3c->ToTH1(potp3c);
84  TH1 *hSpectrumP3d = FDSpectrumP3d->ToTH1(potp3d);
85 
86  TH1 *hSpectrumP3 = FDSpectrumP3b->ToTH1(potp3b);
87  hSpectrumP3->Add(hSpectrumP3c);
88  hSpectrumP3->Add(hSpectrumP3d);
89 
90  TH1 *hSpectrum = FDSpectrumP1->ToTH1(potp1);
91  hSpectrum->Add(hSpectrumP2);
92  hSpectrum->Add(hSpectrumP3);
93 
94  float livetimep1 = FDSpectrumP1->Livetime();
95  float livetimep2 = FDSpectrumP2->Livetime();
96  float livetimep3 = FDSpectrumP3b->Livetime() + FDSpectrumP3c->Livetime() + FDSpectrumP3d->Livetime();
97  float livetime = livetimep1 + livetimep2 + livetimep3;
98 
99  std::cout <<"Loaded numu data: "
100  << "For a total of " << hSpectrum->GetSumOfWeights() << " events and " << pot << " POT and " << livetime << " s of livetime. " << std::endl;
101 
102  Spectrum * FDSpectrum = new Spectrum(hSpectrum, pot, livetime);
103 
104  TFile cosfile((FindCAFAnaDir()+"/numu/SA/SACosmics.root").c_str(),"read"); // cosmic estimate, spectrum from cosmic triggers, normalized by numi out of time
105  TH1F *hcosp1 = (TH1F*)cosfile.Get("hp1_all");
106  TH1F *hcosp2 = (TH1F*)cosfile.Get("hp2_all");
107  TH1F *hcosp3 = (TH1F*)cosfile.Get("hp3_all");
108 
109  TH1F *hcosmics = (TH1F*)cosfile.Get("hp1_all");
110  hcosmics->Add(hcosp2);
111  hcosmics->Add(hcosp3);
112 
113  std::cout << "Taking cosmic trigger spectrum, normalized to numi out of time data, which estimates: " << hcosmics->GetSumOfWeights() << " events total " << std::endl;
114  Spectrum cosmics(hcosmics,pot,livetime);
115 
116  //nue style
117  std::cout << "Loading Numu prediction p123b+p3c \n";
118 
119  TFile input("/nova/ana/users/ecatanom/SecondAna/Predictions/pred2_numu_sa_123b_3c.root", "READ");
120 
121  auto predSACCp1 = PredictionInterp::LoadFrom((TDirectory*)input.Get("numu_pred_p123b")).release();
122  auto predSACCp3c = PredictionInterp::LoadFrom((TDirectory*)input.Get("numu_pred_p3c")).release();
123 
124  PredictionCombinePeriods predSACC ({
125  (std::make_pair(predSACCp1,potp1+potp2+potp3b)),
126  (std::make_pair(predSACCp3c,potp3c+potp3d))});
127 
128  std::cout << "creating single sample experiment " << std::endl;
129  SingleSampleExperiment exptNumu( &(predSACC), *FDSpectrum, cosmics, 0.33/2.81); // todo generalize this
130  std::cout << "done" << std::endl;
131 
132 
133  ////////////////////////////////////////////////////////////
134 
135  TFile * outfile = new TFile (("expts_sa_numu_nue_"+pidStr+"_"+decomp+".root").c_str(),"recreate");
136 
137  exptNue.SaveTo(outfile->mkdir("exptNue"));
138  exptNumu.SaveTo(outfile->mkdir("exptNumu"));
139 
140  std::cout << "output saved to " << outfile->GetName() << std::endl;
141  return;
142 
143 /*
144  //////////////////////////////////////////////////
145  // Add constraints, make experiments
146  //////////////////////////////////////////////////
147 
148  std::cout << "\nCreating multiexperiment\n" << std::endl;
149 
150  const ReactorExperiment react(0.085, 0.005);
151  MultiExperiment exptJoint({&exptNue,&exptNumu,&react});
152 
153  MultiExperiment * exptThis = &exptJoint;
154 
155  ////////////////////////////////////////////////////////////
156  // Systematics
157  std::vector<const ISyst*> surfSysts;
158  GenieRwgtTableSyst kNumuReweightMaNCEL(rwgt::fReweightMaNCEL);
159  GenieRwgtTableSyst kNumuReweightMaCCQEshape(rwgt::fReweightMaCCQEshape);
160  GenieRwgtTableSyst kNumuReweightNormCCQE(rwgt::fReweightNormCCQE);
161  GenieRwgtTableSyst kNumuReweightMaNCRES(rwgt::fReweightMaNCRES);
162  GenieRwgtTableSyst kNumuReweightMvNCRES(rwgt::fReweightMvNCRES);
163  GenieRwgtTableSyst kNumuReweightCCQEPauliSupViaKF(rwgt::fReweightCCQEPauliSupViaKF);
164  GenieRwgtTableSyst kNumuReweightMaCCRES(rwgt::fReweightMaCCRES);
165  GenieRwgtTableSyst kNumuReweightMvCCRES(rwgt::fReweightMvCCRES);
166 
167 
168  surfSysts={
169  &kNueMaCCQEShapeSyst,
170  &kNueCCQENormSyst,
171  &kNueMaCOHSyst,
172  &kNuePauliSyst,
173  &kNueMaRESSyst,
174  &kNueMvRESSyst,
175  &kNueAllPiSyst,
176  &kNueMECSyst,
177  &kNueRPASyst,
178  &kNueBirksCSyst,
179  &kBeamAll,
180  &kNueAbsCalibSyst,
181  &kNueRelCalibSyst,
182  &kNueShiftXSyst,
183  &kNueShiftYSyst,
184  &kNueNormSystBoth,
185  &kNueNormSystSig,
186  &kNueNormSystBkg,
187  &kNumuRelNormSyst,
188  &kNumuNCScaleSyst,
189  &kSAMuEnergyScaleSyst,
190  &kFDSAMuEnergyScaleSyst,
191  &kNumuSummedSmallGENIESyst,
192  &kNumuReweightMaNCEL
193  };
194 
195  exptThis->SetSystCorrelations(0, {
196  {&kBeamAll, &kNueNA49Syst},
197  {&kBeamAll, &kNueAllTransportSyst},
198  {&kNumuRelNormSyst, NULL},
199  {&kNumuNCScaleSyst, NULL},
200  {&kSAMuEnergyScaleSyst, NULL},
201  {&kFDSAMuEnergyScaleSyst, NULL},
202  {&kNumuSummedSmallGENIESyst, NULL},
203  {&kNumuReweightMaNCEL , NULL}
204  });
205 
206  exptThis->SetSystCorrelations(1, {
207  {&kNueMaCCQEShapeSyst, &kNumuReweightMaCCQEshape},
208  {&kNueCCQENormSyst, &kNumuReweightNormCCQE},
209  {&kNuePauliSyst, &kNumuReweightCCQEPauliSupViaKF},
210  {&kNueMaRESSyst, &kNumuReweightMaCCRES},
211  {&kNueMaRESSyst, &kNumuReweightMaNCRES},
212  {&kNueMvRESSyst, &kNumuReweightMvCCRES},
213  {&kNueMvRESSyst, &kNumuReweightMvNCRES},
214  {&kNueBirksCSyst, &kSABirksSyst},
215  {&kNueRelCalibSyst, &kSARelHadESyst},
216  {&kNueMECSyst, &kMECScaleSystSA},
217  {&kNueRPASyst, &kRPACCQESystSA},
218  {&kNueShiftYSyst, &kSACalibYFuncSyst},
219  {&kNueAbsCalibSyst, &kSACalibXYSyst},
220  {&kNueNormSystBoth, &kNumuNormSyst},
221  {&kNueMaCOHSyst, NULL},
222  {&kNueAllPiSyst, NULL},
223  {&kNueShiftXSyst, NULL},
224  {&kNueNormSystSig, NULL},
225  {&kNueNormSystBkg, NULL }
226  });
227  */
228 
229 }
TCut intime("tave>=217.0 && tave<=227.8")
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
void joint_fit_make_experiments(std::string PIDStr="CVN", std::string decomp="prop")
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:335
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
std::string FindCAFAnaDir()
Definition: Utilities.cxx:550
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
#define pot
NuESystPID
Definition: NueSysts.h:12
Loads shifted spectra from files.
virtual void SaveTo(TDirectory *dir) const override
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir)
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
Sum MC predictions from different periods scaled according to data POT targets.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1055
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
FILE * outfile
Definition: dump_event.C:13
Compare a single data spectrum to the MC + cosmics expectation.