ApplyOscillations.C
Go to the documentation of this file.
1 #pragma once
2 
3 // root
4 #include "TFile.h"
5 #include "TString.h"
6 #include "TH1.h"
7 #include "TH2.h"
8 #include "TH3.h"
9 
10 // CAFAna
12 #include "CAFAna/Analysis/Calcs.h"
13 #include "CAFAna/Analysis/Style.h"
15 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/Binning.h"
22 #include "CAFAna/Vars/NusVars.h"
23 #include "OscLib/OscCalcSterile.h"
24 
25 // Selection
33 
34 using namespace ana;
35 
36 // struct to hold plot-related stuff
37 struct PlotInfo {
41  double pot;
42  double livetime;
43 
45  name = n;
46  label = l;
47  detector = d;
48  pot = p;
49  livetime = lt;
50  }
51 };
52 
53 // function to setup sterile oscillator
55 
56 // save function for IPrediction objects
57 void Save(TDirectory* outDir,
58  PlotInfo info,
60  HistType ht,
61  IPrediction* pred = nullptr,
62  Spectrum* cosSpec = nullptr,
63  Spectrum* data = nullptr);
64 
65 // save function for CheatDecomp objects
66 void Save(TDirectory* outDir,
67  PlotInfo info,
68  osc::OscCalcSterile* calc,
69  HistType ht,
70  CheatDecomp* pred = nullptr,
71  Spectrum* data = nullptr);
72 
74 
75  // get running options
76  bool isFD = false, isND = false;
77  if (opt.Contains("FD", TString::kIgnoreCase))
78  isFD = true;
79  if (opt.Contains("ND", TString::kIgnoreCase))
80  isND = true;
81 
82  bool isNDData = opt.Contains("nddata", TString::kIgnoreCase);
83  bool isFDData = opt.Contains("fddata", TString::kIgnoreCase);
84  bool isCosmic = opt.Contains("cos" , TString::kIgnoreCase);
85 
86  // setup POT to scale to
87  double kThisPOT = kAna2020FHCPOT;
88  double kThisLivetime = kAna2020FHCLivetime;
89 
90  // setup output file
91  std::string fileName = inFile.substr(0, inFile.find(".root"));
92  std::string saveName = fileName + "_Oscillated.root";
93 
94  TFile* outFile = new TFile(saveName.c_str(), "recreate");
95 
96  // set up sterile oscillator, but only for 3F
98 
99  // container for spectra
100  std::map<std::string, IPrediction*> fdPredictions;
101  std::map<std::string, CheatDecomp*> ndPredictions;
102  std::map<std::string, Spectrum*> fdData, ndData;
103  std::map<std::string, Spectrum*> fdCos;
104 
105  // get predictions, etc. from file
106  if (isFD){
107 
108  for (size_t iCut = 0; iCut < cutsFD.size(); ++iCut){
109  for (size_t iSub = 0; iSub < subCutsFD.size(); ++iSub){
110 
111  bool isCutData = cutsFD.at(iCut).cutdata && subCutsFD.at(iSub).cutdata;
112 
113  for (size_t iVar = 0; iVar < varsFD.size(); ++iVar){
114 
115  std::string mapLabel = varsFD.at(iVar).name +
116  cutsFD.at(iCut).name +
117  subCutsFD.at(iSub).name;
118 
119  fdPredictions[mapLabel] = nullptr;
120  fdCos [mapLabel] = nullptr;
121  fdData [mapLabel] = nullptr;
122 
123  std::cout << "mapLabel: " << mapLabel << std::endl;
124 
125  fdPredictions[mapLabel] =
126  LoadFromFile<IPrediction>(inFile.c_str(),
127  Form("FD%sPred", mapLabel.c_str())).release();
128  if (isCosmic && isCutData)
129  fdCos[mapLabel] =
130  LoadFromFile<Spectrum>(inFile.c_str(),
131  Form("FD%sCosSpec", mapLabel.c_str())).release();
132 
133  if (isFDData && isCutData)
134  fdData[mapLabel] =
135  LoadFromFile<Spectrum>(inFile.c_str(),
136  Form("FD%sSpec", mapLabel.c_str())).release();
137  }
138 
139  for (size_t iVar = 0; iVar < vars2DFD.size(); ++iVar){
140 
141  std::string mapLabel = vars2DFD.at(iVar).name +
142  cutsFD.at(iCut).name +
143  subCutsFD.at(iSub).name;
144 
145  fdPredictions[mapLabel] = nullptr;
146  fdCos [mapLabel] = nullptr;
147  fdData [mapLabel] = nullptr;
148 
149  std::cout << "mapLabel: " << mapLabel << std::endl;
150 
151  fdPredictions[mapLabel] =
152  LoadFromFile<IPrediction>(inFile.c_str(),
153  Form("FD%sPred", mapLabel.c_str())).release();
154  if (isCosmic && isCutData)
155  fdCos[mapLabel] =
156  LoadFromFile<Spectrum>(inFile.c_str(),
157  Form("FD%sCosSpec", mapLabel.c_str())).release();
158 
159  if (isFDData && isCutData)
160  fdData[mapLabel] =
161  LoadFromFile<Spectrum>(inFile.c_str(),
162  Form("FD%sSpec", mapLabel.c_str())).release();
163  }
164 
165  for (size_t iVar = 0; iVar < vars3DFD.size(); ++iVar){
166 
167  std::string mapLabel = vars3DFD.at(iVar).name +
168  cutsFD.at(iCut).name +
169  subCutsFD.at(iSub).name;
170 
171  fdPredictions[mapLabel] = nullptr;
172  fdCos [mapLabel] = nullptr;
173  fdData [mapLabel] = nullptr;
174 
175  std::cout << "mapLabel: " << mapLabel << std::endl;
176 
177  fdPredictions[mapLabel] =
178  LoadFromFile<IPrediction>(inFile.c_str(),
179  Form("FD%sPred", mapLabel.c_str())).release();
180  if (isCosmic && isCutData)
181  fdCos[mapLabel] =
182  LoadFromFile<Spectrum>(inFile.c_str(),
183  Form("FD%sCosSpec", mapLabel.c_str())).release();
184 
185  if (isFDData && isCutData)
186  fdData[mapLabel] =
187  LoadFromFile<Spectrum>(inFile.c_str(),
188  Form("FD%sSpec", mapLabel.c_str())).release();
189  }
190 
191  }
192  }
193  }
194 
195  if (isND){
196 
197  for (size_t iCut = 0; iCut < cutsND.size(); ++iCut){
198  for (size_t iSub = 0; iSub < subCutsND.size(); ++iSub){
199 
200  bool isCutData = cutsND.at(iCut).cutdata && subCutsND.at(iSub).cutdata;
201 
202  for (size_t iVar = 0; iVar < varsND.size(); ++iVar){
203 
204  std::string mapLabel = varsND.at(iVar).name +
205  cutsND.at(iCut).name +
206  subCutsND.at(iSub).name;
207 
208  std::cout << "mapLabel: " << mapLabel << std::endl;
209 
210  ndPredictions[mapLabel] = nullptr;
211  ndData [mapLabel] = nullptr;
212 
213  ndPredictions[mapLabel] =
214  LoadFromFile<CheatDecomp>(inFile.c_str(),
215  Form("ND%sDecomp", mapLabel.c_str())).release();
216  if (isNDData && isCutData)
217  ndData[mapLabel] =
218  LoadFromFile<Spectrum>(inFile.c_str(),
219  Form("ND%sSpec", mapLabel.c_str())).release();
220  }
221 
222  for (size_t iVar = 0; iVar < vars2DND.size(); ++iVar){
223 
224  std::string mapLabel = vars2DND.at(iVar).name +
225  cutsND.at(iCut).name +
226  subCutsND.at(iSub).name;
227 
228  std::cout << "mapLabel: " << mapLabel << std::endl;
229 
230  ndPredictions[mapLabel] = nullptr;
231  ndData [mapLabel] = nullptr;
232 
233  ndPredictions[mapLabel] =
234  LoadFromFile<CheatDecomp>(inFile.c_str(),
235  Form("ND%sDecomp", mapLabel.c_str())).release();
236  if (isNDData && isCutData)
237  ndData[mapLabel] =
238  LoadFromFile<Spectrum>(inFile.c_str(),
239  Form("ND%sSpec", mapLabel.c_str())).release();
240  }
241 
242  for (size_t iVar = 0; iVar < vars3DND.size(); ++iVar){
243 
244  std::string mapLabel = vars3DND.at(iVar).name +
245  cutsND.at(iCut).name +
246  subCutsND.at(iSub).name;
247 
248  std::cout << "mapLabel: " << mapLabel << std::endl;
249 
250  ndPredictions[mapLabel] = nullptr;
251  ndData [mapLabel] = nullptr;
252 
253  ndPredictions[mapLabel] =
254  LoadFromFile<CheatDecomp>(inFile.c_str(),
255  Form("ND%sDecomp", mapLabel.c_str())).release();
256  if (isNDData && isCutData)
257  ndData[mapLabel] =
258  LoadFromFile<Spectrum>(inFile.c_str(),
259  Form("ND%sSpec", mapLabel.c_str())).release();
260  }
261 
262  }
263  }
264  }
265 
266  // --------------------
267  // saving to file
268  // --------------------
269 
270  if (isFD){
271 
272  oscCalc->SetL(810);
273 
274  // first get POT and livetime from data, if it exists
275  std::string caloForPOT = "kCaloEPreselTotal";
276  Spectrum* thisData = nullptr;
277  if (isFDData)
278  thisData = LoadFromFile<Spectrum>(inFile.c_str(),
279  Form("FD%sSpec", caloForPOT.c_str())).release();
280  double thisPOT = thisData ? thisData->POT() : kThisPOT;
281  double thisLivetime = thisData ? thisData->Livetime() : kThisLivetime;
282 
283  std::string detector = "FD";
284  TDirectory* fdDir = outFile->mkdir(detector.c_str(), detector.c_str());
285 
286  for (size_t iCut = 0; iCut < cutsFD.size(); ++iCut){
287  for (size_t iSub = 0; iSub < subCutsFD.size(); ++iSub){
288  for (size_t iVar = 0; iVar < varsFD.size(); ++iVar){
289 
290  std::string totLabel = varsFD.at(iVar).name +
291  cutsFD.at(iCut).name +
292  subCutsFD.at(iSub).name;
293 
294  std::cout << "totLabel: " << totLabel << std::endl;
295 
296  Save(fdDir,
297  PlotInfo(totLabel, totLabel, detector, thisPOT, thisLivetime),
298  oscCalc,
300  fdPredictions[totLabel],
301  fdCos[totLabel],
302  fdData[totLabel]);
303  }
304  for (size_t iVar = 0; iVar < vars2DFD.size(); ++iVar){
305 
306  std::string totLabel = vars2DFD.at(iVar).name +
307  cutsFD.at(iCut).name +
308  subCutsFD.at(iSub).name;
309 
310  std::cout << "totLabel: " << totLabel << std::endl;
311 
312  Save(fdDir,
313  PlotInfo(totLabel, totLabel, detector, thisPOT, thisLivetime),
314  oscCalc,
316  fdPredictions[totLabel],
317  fdCos[totLabel],
318  fdData[totLabel]);
319  }
320  for (size_t iVar = 0; iVar < vars3DFD.size(); ++iVar){
321 
322  std::string totLabel = vars3DFD.at(iVar).name +
323  cutsFD.at(iCut).name +
324  subCutsFD.at(iSub).name;
325 
326  std::cout << "totLabel: " << totLabel << std::endl;
327 
328  Save(fdDir,
329  PlotInfo(totLabel, totLabel, detector, thisPOT, thisLivetime),
330  oscCalc,
332  fdPredictions[totLabel],
333  fdCos[totLabel],
334  fdData[totLabel]);
335  }
336  }
337  }
338  }
339 
340  if (isND){
341 
342  oscCalc->SetL(1);
343 
344  // first get POT and livetime from data, if it exists
345  std::string caloForPOT = "kCaloEPreselTotal";
346  Spectrum* thisData = nullptr;
347  if (isNDData)
348  thisData = LoadFromFile<Spectrum>(inFile.c_str(),
349  Form("ND%sSpec", caloForPOT.c_str())).release();
350  double thisPOT = thisData ? thisData->POT() : kThisPOT;
351  double thisLivetime = thisData ? thisData->Livetime() : kThisLivetime;
352 
353  std::string detector = "ND";
354  TDirectory* ndDir = outFile->mkdir(detector.c_str(), detector.c_str());
355 
356  for (size_t iCut = 0; iCut < cutsND.size(); ++iCut){
357  for (size_t iSub = 0; iSub < subCutsND.size(); ++iSub){
358  for (size_t iVar = 0; iVar < varsND.size(); ++iVar){
359 
360  std::string totLabel = varsND.at(iVar).name +
361  cutsND.at(iCut).name +
362  subCutsND.at(iSub).name;
363 
364  std::cout << "totLabel: " << totLabel << std::endl;
365 
366  Save(ndDir,
367  PlotInfo(totLabel, totLabel, detector, thisPOT, thisLivetime),
368  oscCalc,
370  ndPredictions[totLabel],
371  ndData[totLabel]);
372  }
373  for (size_t iVar = 0; iVar < vars2DND.size(); ++iVar){
374 
375  std::string totLabel = vars2DND.at(iVar).name +
376  cutsND.at(iCut).name +
377  subCutsND.at(iSub).name;
378 
379  std::cout << "totLabel: " << totLabel << std::endl;
380 
381  Save(ndDir,
382  PlotInfo(totLabel, totLabel, detector, thisPOT, thisLivetime),
383  oscCalc,
385  ndPredictions[totLabel],
386  ndData[totLabel]);
387  }
388  for (size_t iVar = 0; iVar < vars3DND.size(); ++iVar){
389 
390  std::string totLabel = vars3DND.at(iVar).name +
391  cutsND.at(iCut).name +
392  subCutsND.at(iSub).name;
393 
394  std::cout << "totLabel: " << totLabel << std::endl;
395 
396  Save(ndDir,
397  PlotInfo(totLabel, totLabel, detector, thisPOT, thisLivetime),
398  oscCalc,
400  ndPredictions[totLabel],
401  ndData[totLabel]);
402  }
403 
404  }
405  }
406  }
407 }
408 
410 
411  osc::OscCalcSterile* calc = DefaultSterileCalc(nFlavours);
412 
413  calc->SetL(810);
414  calc->SetAngle(1, 2, asin(sqrt(.307)));
415  calc->SetAngle(1, 3, asin(sqrt(.021)));
416  calc->SetAngle(2, 3, asin(sqrt(.56)));
417  calc->SetDm(2, 7.53e-5);
418  calc->SetDm(3, 2.48e-3 + 7.53e-5); // Dm31 = Dm32 + Dm21
419  calc->SetDelta(1, 3, 0);
420 
421  if (nFlavours > 3){
422  calc->SetDm (4, 5.0);
423  calc->SetAngle(1, 4, 34);
424  calc->SetAngle(2, 4, 34);
425  calc->SetAngle(3, 4, 34);
426  calc->SetDelta(1, 4, 0);
427  calc->SetDelta(2, 4, 0);
428  }
429 
430  return calc;
431 
432 }
433 
434 void Save(TDirectory* outDir,
435  PlotInfo info,
436  osc::OscCalcSterile* calc,
437  HistType ht,
438  IPrediction* pred ,
439  Spectrum* cosmic ,
440  Spectrum* data ){
441 
442  outDir->cd();
443 
444  double thisPOT = info.pot;
445  double thisLivetime = info.livetime;
446 
447  TH1 *hTotalMC = nullptr;
448  TH1 *hData = nullptr;
449  TH1 *hCosmic = nullptr;
450  if (ht == HistType::kTH1){
451  hTotalMC = pred->Predict(calc).ToTH1(thisPOT);
452  if (data)
453  hData = data->ToTH1(thisPOT);
454  if (cosmic){
455  hCosmic = cosmic->ToTH1(thisLivetime, kLivetime);
456  }
457  }
458  if (ht == HistType::kTH2){
459  hTotalMC = pred->Predict(calc).ToTH2(thisPOT);
460  if (data)
461  hData = data->ToTH2(thisPOT);
462  if (cosmic){
463  hCosmic = cosmic->ToTH2(thisLivetime, kLivetime);
464  }
465  }
466  if (ht == HistType::kTH3){
467  hTotalMC = pred->Predict(calc).ToTH3(thisPOT);
468  if (data)
469  hData = data->ToTH3(thisPOT);
470  if (cosmic){
471  hCosmic = cosmic->ToTH3(thisLivetime, kLivetime);
472  }
473  }
474 
475  hTotalMC->Write((info.name+"MC").c_str());
476  if (cosmic)
477  hCosmic->Write((info.name+"Cosmic").c_str());
478  if (data)
479  hData->Write((info.name+"Data").c_str());
480 
481 }
482 
483 void Save(TDirectory* outDir,
484  PlotInfo info,
485  osc::OscCalcSterile* calc,
486  HistType ht,
488  Spectrum* data ){
489 
490  outDir->cd();
491 
492  double thisPOT = info.pot;
493  double thisLivetime = info.livetime;
494 
495  TH1 *hTotal = nullptr;
496  TH1 *hData = nullptr;
497 
498  if (ht == HistType::kTH1)
499  hTotal = decomp->TotalMC().ToTH1(thisPOT);
500  if (ht == HistType::kTH2)
501  hTotal = decomp->TotalMC().ToTH2(thisPOT);
502  if (ht == HistType::kTH3)
503  hTotal = decomp->TotalMC().ToTH3(thisPOT);
504 
505  hTotal->Write((info.name+"MC").c_str());
506  if (data)
507  hData = data->ToTH1(thisPOT);
508  if (data)
509  hData->Write((info.name+"Data" ).c_str());
510 
511 }
const XML_Char XML_Encoding * info
Definition: expat.h:530
osc::OscCalcSterile * SterileOscillator(int nFlavours)
fileName
Definition: plotROC.py:78
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double lt
std::vector< CutContainer > cutsFD
Definition: CutFlow.h:23
std::vector< CutContainer > subCutsND
Definition: CutFlow.h:46
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:148
void ApplyOscillations(std::string inFile, TString opt)
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
std::string label
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
void Save(TDirectory *outDir, PlotInfo info, osc::OscCalcSterile *calc, HistType ht, IPrediction *pred=nullptr, Spectrum *cosSpec=nullptr, Spectrum *data=nullptr)
void SetDelta(int i, int j, double delta)
Adapt the PMNS_Sterile calculator to standard interface.
Definition: Structs.h:42
std::vector< CutContainer > subCutsFD
Definition: CutFlow.h:40
osc::OscCalcDumb calc
std::string outDir
std::vector< CutContainer > cutsND
Definition: CutFlow.h:32
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
double livetime
std::string name
const XML_Char const XML_Char * data
Definition: expat.h:268
ifstream inFile
Definition: AnaPlotMaker.h:34
Definition: Structs.h:43
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
std::vector< VarContainer3D > vars3DFD
Definition: CutVariables.h:74
Sum up livetimes from individual cosmic triggers.
TFile * outFile
Definition: PlotXSec.C:135
Definition: Structs.h:41
const double kAna2020FHCLivetime
Definition: Exposures.h:237
std::vector< VarContainer2D > vars2DFD
Definition: CutVariables.h:63
Float_t d
Definition: plot.C:236
PlotInfo(std::string n, std::string l, std::string d, double p, double lt)
const double kAna2020FHCPOT
Definition: Exposures.h:233
std::vector< VarContainer > varsND
Definition: CutVariables.h:82
double POT() const
Definition: Spectrum.h:219
OStream cout
Definition: OStream.cxx:6
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:218
void SetDm(int i, double dm)
std::vector< VarContainer3D > vars3DND
Definition: CutVariables.h:109
virtual void SetL(double L) override
std::vector< VarContainer2D > vars2DND
Definition: CutVariables.h:103
std::string detector
Spectrum TotalMC()
Definition: CheatDecomp.h:35
std::vector< VarContainer > varsFD
Definition: CutVariables.h:28
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
HistType
This code creates ratio plots of calibration samples.
Float_t e
Definition: plot.C:35
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:222
T asin(T number)
Definition: d0nt_math.hpp:60