SpectrumParamEffectsAna.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
7 
9 
10 #include "TCanvas.h"
11 #include "TDirectory.h"
12 #include "TFile.h"
13 #include "TH1.h"
14 #include "TLegend.h"
15 
16 #include <string>
17 
18 using namespace ana;
19 
20 double CalcParamEffect(double nom, double shift);
21 
22 void PlotEvents(TDirectory* out, std::string label,
23  TH1* nom, TH1* hup, TH1* hdn);
24 
25 void PrintEvents(FILE* text, std::string label,
26  double nom, double up, double dn);
27 
29 
30 // After running a fit, reset the values in a calculator to preset defaults
32 
34 {
35  TH1::AddDirectory(0);
36 
37  // Set up path to file with filled CAFAna objects and for output of analysis
38  std::string folder = "";
39  std::string filenm = "SpectrumParamEffects";
40 
41  std::string loadLocation = "/nova/ana/steriles/Ana01/FitSystEffects.root";
42  std::string saveLocation = folder + filenm + "Ana.root";
43  std::string textLocation = folder + filenm + "Ana.txt";
44 
45  // Load filled objects from file
46  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
47  TDirectory* tmpL = gDirectory;
48  TDirectory* loadDir = gDirectory;
49  loadDir->cd((loadLocation + ":/pred").c_str());
51 
52  loadDir->cd((loadLocation + ":/sCosmic").c_str());
53  Spectrum sCos = *Spectrum::LoadFrom(gDirectory);
55 
56  rootL->Close(); // Don't forget to close the file!
57 
58  double nCos = sCos.Integral(sCos.POT());
59  std::cout << "Cos spectrum: " << sCos.POT()
60  << " and " << nCos << " events" << std::endl;
61 
63  calc4f->SetNFlavors(4);
64  ResetAngles(calc4f);
65 
66  std::map<std::string, std::pair<double, double> > parammap;
67  parammap["t13"] = {-0.3675*M_PI/180., 0.3714*M_PI/180.};
68  parammap["d13"] = {-0.27, 0.38};
69  parammap["m32"] = {-0.00006, 0.00006};
70 
71  std::map<std::string, TH1*> hups;
72  std::map<std::string, TH1*> hdns;
73  std::map<std::string, double> nups;
74  std::map<std::string, double> ndns;
75 
76  Spectrum snom = pred.Predict(calc4f) + sCos;
77  TH1* hnom = snom.ToTH1(kSecondAnaPOT);
78  double nnom = snom.Integral(kSecondAnaPOT);
79  double scale = 83.509/nnom;
80  nnom *= scale;
81 
82  // Open files for saving analysis output
83  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
84  FILE* textF;
85  textF = fopen(textLocation.c_str(), "w");
86 
87  for(const auto& parampair : parammap) {
88  std::string label = parampair.first;
89  double upshift = parampair.second.second;
90  double dnshift = parampair.second.first;
91  ResetAngles(calc4f);
92 
93  SetCalcParam(calc4f, label, upshift);
94  Spectrum sup = pred.Predict(calc4f) + sCos;
95  hups[label] = sup.ToTH1(kSecondAnaPOT);
96  nups[label] = sup.Integral(kSecondAnaPOT)*scale;
97 
98  ResetAngles(calc4f);
99 
100  SetCalcParam(calc4f, label, dnshift);
101  Spectrum sdn = pred.Predict(calc4f) + sCos;
102  hdns[label] = sdn.ToTH1(kSecondAnaPOT);
103  ndns[label] = sdn.Integral(kSecondAnaPOT)*scale;
104 
105  PrintEvents(textF, label, nnom, nups[label], ndns[label]);
106  PlotEvents( rootF, label, hnom, hups[label], hdns[label]);
107  }
108 
109  fclose(textF); // Text file is done now, close it
110  rootF->Close(); // Close the output root file
111 }
112 
113 //------------------------------------------------------------------------------
114 double CalcParamEffect(double nom, double shift)
115 {
116  return 100.*(shift - nom)/nom;
117 }
118 
119 //------------------------------------------------------------------------------
120 void PlotEvents(TDirectory* out, std::string label,
121  TH1* nom, TH1* hup, TH1* hdn)
122 {
123  double max = nom->GetBinContent(nom->GetMaximumBin());
124  max = std::max(max, hup->GetBinContent(hup->GetMaximumBin()));
125  max = std::max(max, hdn->GetBinContent(hdn->GetMaximumBin()));
126 
127  std::string title = ";Calorimetric Energy (GeV);Events / 0.25 GeV";
128 
129  SetHistOptions(nom, max, title, -404, kBlack, false);
130  SetHistOptions(hup, max, title, -404, kRed, false);
131  SetHistOptions(hdn, max, title, -404, kBlue, false);
132 
133  hup->SetLineStyle(2);
134  hdn->SetLineStyle(3);
135 
136  std::string ptype = label.substr(0,1);
137 
138  std::string param = "";
139  if(ptype.compare("t") == 0) {
140  param = "#theta_{" + label.substr(1,2) + "} ";
141  }
142  if(ptype.compare("d") == 0) {
143  param = "#delta_{" + label.substr(1,2) + "} ";
144  }
145  if(ptype.compare("m") == 0) {
146  param = "#Deltam^{2}_{" + label.substr(1,2) + "} ";
147  }
148 
149  double xL = 0.6, xR = 0.8, yB = 0.65, yT = 0.85;
150  TLegend* leg = new TLegend(xL, yB, xR, yT);
151  SetLegendOptions(leg);
152  leg->AddEntry(nom, "Nominal", "l");
153  leg->AddEntry(hup, (param + "+1#sigma").c_str(), "l");
154  leg->AddEntry(hdn, (param + "-1#sigma").c_str(), "l");
155  leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);
156 
157  std::string ctitle = "c" + label;
158  TCanvas* c = new TCanvas(ctitle.c_str(), ctitle.c_str(), 800, 500);
159  nom->Draw("hist");
160  hup->Draw("hist same");
161  hdn->Draw("hist same");
162  leg->Draw();
163 
164  gPad->RedrawAxis();
165  out->WriteTObject(c);
166  return;
167 }
168 
169 //------------------------------------------------------------------------------
171  double nom, double up, double dn)
172 {
173  fprintf(text, "Number of events for %s:\n", label.c_str());
174  fprintf(text, "Nom: %4.2f, Up: %4.2f, Down: %4.2f\n", nom, up, dn);
175  fprintf(text, "Shifts: Up: %4.2f%%, Down: %4.2f%%\n\n",
176  CalcParamEffect(nom, up),
177  CalcParamEffect(nom, dn));
178  return;
179 }
180 
181 //------------------------------------------------------------------------------
183 {
185  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
186  calc->SetAngle(3, 4, 20.*M_PI/180.);
187  calc->SetAngle(2, 4, 20.*M_PI/180.);
188 
189  return;
190 }
191 
192 //------------------------------------------------------------------------------
194 {
195  std::string ptype = label.substr(0,1);
196  int val1 = std::stoi(label.substr(1,1));
197  int val2 = std::stoi(label.substr(2,1));
198 
199  if(val1 > val2) { std::swap(val1, val2); }
200 
201  if(ptype.compare("t") == 0) {
202  calc->SetAngle(val1, val2, calc->GetAngle(val1, val2) + shift);
203  }
204  else if(ptype.compare("d") == 0) {
205  calc->SetDelta(val1, val2, calc->GetDelta(val1, val2) + shift);
206  }
207  else if(ptype.compare("m") == 0) {
208  if(val2 != 1) {
209  calc->SetDm(val2, calc->GetDm(val2) + shift);
210  }
211  }
212  else {
213  std::cout << "Unknown parameter type." << std::endl;
214  }
215 
216  return;
217 }
void SetHistOptions(TH1 *h, double max, std::string title, int ndiv, Color_t col, bool fill)
Set common options for a TLegend.
T max(const caf::Proxy< T > &a, T b)
enum BeamMode kRed
double GetDelta(int i, int j) const
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
void ResetAngles(osc::OscCalcSterile *calc)
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 OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:233
void SetDelta(int i, int j, double delta)
Adapt the PMNS_Sterile calculator to standard interface.
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
void SpectrumParamEffectsAna()
osc::OscCalcDumb calc
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
fclose(fg1)
double CalcParamEffect(double nom, double shift)
Double_t scale
Definition: plot.C:25
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
void PrintEvents(FILE *text, std::string label, double nom, double up, double dn)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
double POT() const
Definition: Spectrum.h:227
OStream cout
Definition: OStream.cxx:6
void SetAngle(int i, int j, double th)
void SetLegendOptions(TLegend *leg)
Set common options for a TLegend.
void SetDm(int i, double dm)
double GetDm(int i) const
Sum MC predictions from different periods scaled according to data POT targets.
void SetCalcParam(osc::OscCalcSterile *calc, std::string label, double shift)
const double kSecondAnaPOT
Definition: Exposures.h:87
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kBlue
void PlotEvents(TDirectory *out, std::string label, TH1 *nom, TH1 *hup, TH1 *hdn)
double GetAngle(int i, int j) const
Spectrum Predict(osc::IOscCalc *calc) const override
enum BeamMode string