AnaResultsSpectra.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TDirectory.h"
3 #include "TFile.h"
4 #include "TH1.h"
5 #include "TKey.h"
6 #include "TLatex.h"
7 #include "TLegend.h"
8 #include "TPad.h"
9 
10 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Analysis/Style.h"
14 
15 #include "OscLib/OscCalcSterile.h"
16 
18 
19 using namespace ana;
20 
21 // Definitions for organization
22 
24 
25 void CalcChi2(const Spectrum& sAll, const Spectrum& sData, FILE* textOFS);
26 
27 void LoadMaps(TDirectory* dir,
28  std::map<std::string, PredictionCombinePeriods*>* preds,
29  std::map<std::string, Spectrum*>* cosmics,
30  std::map<std::string, Spectrum*>* dataspecs);
31 
32 // Plot spectra and stack canvases for FD
33 // Sets up call to the workhorse version of PlotStack
34 // \param pred The FD prediction object
35 // \param sCos The FD cosmic prediction spectrum
36 // \param calc Sterile calculator assuming 3 flavor oscillations
37 // \param calcSt Sterile calculator assuming 4 flavor oscillations
38 // \param rootOut Already open file for root output
39 // \param textOFS Already open filestream for text output
40 // \param name A name for the output, only part of the saved output object's name
41 // \param title Title for the histograms. Typically this is the x axis label
42 // \param sterile Pointer to a data spectrum (optional)
44  Spectrum* sCos,
46  TDirectory* rootOut, FILE* textOFS,
47  std::string name, Spectrum* sData = nullptr);
48 
50 {
51  TH1::AddDirectory(0);
52 
53  std::map<std::string, PredictionCombinePeriods*> preds;
54  std::map<std::string, Spectrum*> cosmics;
55  std::map<std::string, Spectrum*> dataspecs;
56 
57 
58  // Load filled objects from file
59  TFile* rootL = new TFile(filename.c_str(), "READ");
60  LoadMaps(rootL, &preds, &cosmics, &dataspecs);
61  rootL->Close(); // Don't forget to close the file!
62 
63  // Set up oscillation calculator that uses default 3 flavor parameters
65 
66  std::string labelEReco = "Calorimetric Energy (GeV)";
67 
68  // Open files for saving analysis output
69  TFile* rootF = new TFile("ana01_fd_spectra.root", "RECREATE");
70  FILE* textF;
71  textF = fopen("ana01_fd_spectra.txt", "w");
72 
73  for(const auto& varobjpair : preds) {
74  std::string label = varobjpair.first;
75 
76  PlotStack(preds[label], cosmics[label],
77  calc3f, rootF, textF, label);
78  PlotStack(preds[label], cosmics[label], calc3f,
79  rootF, textF, label + "DataMC",
80  dataspecs[label]);
81  }
82 
83  fclose(textF); // Text file is done now, close it
84  rootF->Close(); // Close the output root file
85 }
86 
87 //------------------------------------------------------------------------------
89 {
91 
92  int div = -404;
93  double xL = 0.5, xR = 0.7, yB = 0.62, yT = 0.87;
94  double texX = xL;
95  double texY = yB - 0.045;
96 
97  if(name.find("CalE") != std::string::npos) {
98  div = -404;
99  }
100  if(name.find("VtxX") != std::string::npos) {
101  div = -404;
102  xL = 0.125;
103  xR = 0.325;
104  texY = yT - 0.055;
105  }
106  if(name.find("VtxY") != std::string::npos) {
107  div = -404;
108  xL = 0.625;
109  xR = 0.825;
110  texX = xL + 0.025;
111  }
112  if(name.find("VtxZ") != std::string::npos) {
113  div = -803;
114  xL = 0.15;
115  xR = 0.35;
116  texY = yT - 0.055;
117  }
118  if(name.find("CVN") != std::string::npos) {
119  div = -505;
120  xL = 0.125;
121  xR = 0.325;
122  texX = xL + 0.025;
123  }
124  if(name.find("NHit") != std::string::npos) {
125  div = -505;
126  }
127  if(name.find("CosBDT") != std::string::npos) {
128  div = -503;
129  xL = 0.125;
130  xR = 0.325;
131  texX = xL;
132  }
133  if(name.find("PTP") != std::string::npos) {
134  div = -505;
135  xL = 0.625;
136  xR = 0.825;
137  texX = 0.15;
138  texY = yT - 0.055;
139  }
140  if(name.find("EperHit") != std::string::npos) {
141  div = -505;
142  }
143  if(name.find("DistTop") != std::string::npos) {
144  div = -404;
145  xL = 0.125;
146  xR = 0.325;
147  texY = yT - 0.055;
148  }
149 
150 
151  ret.ndiv = div;
152  ret.xL = xL;
153  ret.xR = xR;
154  ret.yT = yT;
155  ret.yB = yB;
156  ret.texX = texX;
157  ret.texY = texY;
158 
159  return ret;
160 }
161 
162 //------------------------------------------------------------------------------
163 void CalcChi2(const Spectrum& sAll, const Spectrum& sData, FILE* textOFS)
164 {
165  TH1* hAll = sAll.ToTH1(kSecondAnaPOT);
166  TH1* hData = sData.ToTH1(sData.POT());
167  double chi2 = LogLikelihood(hAll, hData);
168  fprintf(textOFS, "The chi2 is %.2f\n\n", chi2);
169 }
170 
171 //------------------------------------------------------------------------------
172 /// Load CAFAna objects from a TDirectory
173 /// The maps to contain these objects must already exist
174 void LoadMaps(TDirectory* dir,
175  std::map<std::string, PredictionCombinePeriods*>* preds,
176  std::map<std::string, Spectrum*>* cosmics,
177  std::map<std::string, Spectrum*>* dataspecs)
178 {
179  TDirectory* tmp = gDirectory;
180  TDirectory* keyDir = gDirectory;
181 
182  // Loop over all keys in this directory
183  dir->cd();
184  TIter nextkey(dir->GetListOfKeys());
185  TKey* key;
186  TKey* oldkey = 0;
187 
188  std::string sep = "__"; // This is the token separator
189  std::string cafanaType = "";
190  std::string sampleLabel = "";
191 
192  while((key = (TKey*)nextkey())) {
193  // Keep only the highest cycle number for each key
194  if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
195  if(key->ReadObj()->IsFolder()) {
196  keyDir = (TDirectory*)key->ReadObj();
197  }
198  else {
199  continue;
200  }
201 
202  std::string name(key->GetName());
203  std::vector<std::string> tokens;
204  std::size_t pos = 0, prev = 0;
205  while(pos != std::string::npos) { // Loop over the TObject name
206  pos = name.find(sep, prev); // Find the next separator
207 
208  if(pos != std::string::npos) {
209  tokens.push_back(name.substr(prev, pos-prev));
210  prev = pos+2; // Move past the current separator
211  }
212  else {
213  tokens.push_back(name.substr(prev, name.length()-prev));
214  }
215  }
216 
217  cafanaType = tokens[0];
218  sampleLabel = tokens[1];
219  if(cafanaType.compare("pred") == 0) {
220  (*preds)[sampleLabel] = PredictionCombinePeriods::LoadFrom(keyDir).release();
221  }
222  if(cafanaType.compare("cosmic") == 0) {
223  (*cosmics)[sampleLabel] = Spectrum::LoadFrom(keyDir).release();
224  }
225  if(cafanaType.compare("dataspec") == 0) {
226  (*dataspecs)[sampleLabel] = Spectrum::LoadFrom(keyDir).release();
227  }
228  }
229 
230  tmp->cd();
231  return;
232 }
233 
234 //------------------------------------------------------------------------------
237  TDirectory* rootOut, FILE* textOFS,
238  std::string name, Spectrum* sData)
239 {
240  if(!sCos->POT()) {
241  sCos->OverridePOT(kSecondAnaPOT);
242  }
243 
244  Spectrum spectraFD[MAXSPEC] = {
245  pred->Predict(calc) + *sCos,
249  *sCos
250  };
251 
252  std::string det = "FD";
253  PlotStack(spectraFD, rootOut, textOFS, name, det,
254  kSecondAnaPOT, 6.05, GetPlotOptions(name), sData);
255  if(sData) {
256  CalcChi2(spectraFD[0], *sData, textOFS);
257  }
258 
259  return;
260 }
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string det, double POT, double potEquiv, PlotOptions opt, Spectrum *dataspec)
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double CalcChi2(TH1 *hmc, TH1 *hdata)
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
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
#define MAXSPEC
Definition: NusLoadProd3.h:18
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:225
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
Adapt the PMNS_Sterile calculator to standard interface.
string filename
Definition: shutoffs.py:106
Float_t tmp
Definition: plot.C:36
osc::OscCalcDumb calc
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
double chi2()
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
fclose(fg1)
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
void prev()
Definition: show_event.C:91
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
double POT() const
Definition: Spectrum.h:219
const char sep
PlotOptions GetPlotOptions(std::string name)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void LoadMaps(TDirectory *dir, std::map< std::string, IDecomp * > *nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > *shifted)
Definition: PPFXHelper.h:273
Sum MC predictions from different periods scaled according to data POT targets.
void AnaResultsSpectra(std::string filename)
All neutrinos, any flavor.
Definition: IPrediction.h:26
const double kSecondAnaPOT
Definition: Exposures.h:87
Spectrum Predict(osc::IOscCalc *calc) const override