DrawSurface.C
Go to the documentation of this file.
1 /// DrawSurface.C by J. Hewes <jhewes15@fnal.gov>
2 
3 #ifdef __CINT__
4 
5 #include <iostream>
6 
7 void DrawSurface(std::string file, std::string type, double dm41 = -1) {
8  throw std::runtime_error("You must run in compiled mode!");
9 }
10 
11 #else
12 
13 #include "CAFAna/Core/Progress.h"
17 
18 using namespace ana;
19 
20 std::vector<double> LogBins(int nbins,double xlo, double xhi)
21 {
22  double log_spacing = exp( (log(xhi) - log(xlo))/(nbins) );
23  vector<double> binning;
24  for (int i = 0; i <= nbins; ++i) {
25  if (i == 0) binning.push_back(xlo);
26  else binning.push_back(log_spacing*binning[i-1]);
27  }
28  return binning;
29 }
30 
31 double GetBinValue(TFile* f, std::string name, int i, int j)
32 {
33  DontAddDirectory guard;
34  TH2D* surface = (TH2D*)f->Get(name.c_str());
35  double val = surface->GetBinContent(i, j);
36  delete surface;
37  return val;
38 }
39 
40 void DrawSurface(std::string file, std::string type, bool rhc = false,
41  double dm41 = -1, std::string detectors_in = "both", int universe = 1) {
42 
43  DontAddDirectory guard;
44  FormatCanvas();
45 
46  // Deal with options
47  if (type != "th24vsth34" && type != "th24vsdm41" && type != "th34vsdm41")
48  throw std::runtime_error("Plot type \"" + type + "\" not understood!");
49 
50  TFile* surface_file = new TFile(file.c_str(), "read");
51  TH2D* surface = (TH2D*)surface_file->Get("surface");
52 
53  if (type == "th24vsth34" && dm41 == -1)
54  throw std::runtime_error("For a th24 vs th34 surface, you must specify dm41!");
55 
56  // Check on detectors
57  bool use_nd;
58  bool use_fd;
59  if (detectors_in == "nd") {
60  use_nd = true;
61  use_fd = false;
62  }
63  else if (detectors_in == "fd") {
64  use_nd = false;
65  use_fd = true;
66  }
67  else if (detectors_in == "both") {
68  use_nd = true;
69  use_fd = true;
70  }
71  else throw std::runtime_error("The detectors_in string must be \"nd\", \"fd\" or \"both\".");
72 
73  // Polarity flag & exposures
74  std::string polarity = rhc? "RHC" : "FHC";
76  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
77  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
78 
79  // Now put it all together
80  PredictionConcat* sim = GetDefaultSimulation(rhc, detectors_in);
81  double pot = use_fd? fd_pot : nd_pot;
82  double livetime = use_fd? fd_livetime : 0;
83 
84  // Merged cosmic spectrum
85  std::vector<TH1D*> cosmics = GetCosmics(rhc, detectors_in);
86  Spectrum cosmic(sim->MergeHists(cosmics), pot, livetime);
87 
88  // Create sterile oscillation calculator
90  calc->SetNFlavors(4);
91  SetAngles(calc);
92 
93  // Load fake data spectrum and add in cosmics
94  std::string fakedata_filename = rhc? "rhc" : "fhc";
95  fakedata_filename += "_fake_data";
96  Spectrum fake_data(LoadFakeData(fakedata_filename, universe, pot, livetime));
97 
98  // No point adding cosmics if there's no far detector
99  if (use_fd) fake_data += cosmic;
100 
101  // Initialise output file
102  std::string output_name = type + "_ana.root";
103  TFile* output_file = new TFile(output_name.c_str(), "recreate");
104  TCanvas* c = new TCanvas();
105  c->SetLogy();
106 
107  Progress p("Drawing spectrum for each surface point.");
108  double n_loop = 0;
109 
110  // Loop over every bin in surface
111  for (int i = 1; i <= 50; ++i) {
112  for (int j = 1; j <= 50; ++j) {
113 
114  p.SetProgress(n_loop / 2500.);
115 
116  // Reset oscillation calculator
117  SetAngles(calc);
118 
119  // Set binning variables according to surface type
120  double xvar = surface->GetXaxis()->GetBinCenter(i);
121  double yvar = surface->GetYaxis()->GetBinCenter(j);
122  if (type == "th24vsth34") {
123  calc->SetDm(4, dm41);
124  calc->SetAngle(3, 4, xvar / 180. * M_PI);
125  calc->SetAngle(2, 4, yvar / 180. * M_PI);
126  }
127  else if (type == "th24vsdm41") {
128  kFitSinSqTheta24Sterile.SetValue(calc, xvar);
129  kFitDmSq41Sterile.SetValue(calc, yvar);
130  }
131  else if (type == "th34vsdm41") {
132  kFitSinSqTheta34Sterile.SetValue(calc, xvar);
133  kFitDmSq41Sterile.SetValue(calc, yvar);
134  }
135 
136  // Draw simulation before fit
137  TH1D* h_sim_before = sim->Predict(calc).ToTH1(pot);
138  h_sim_before->SetLineColor(kBlue);
139  h_sim_before->Draw("hist");
140 
141  // Set best fit points
142  double delta24 = GetBinValue(surface_file, "delta24(pi)", i, j);
143  calc->SetDelta(2, 4, delta24);
144  double th23 = GetBinValue(surface_file, "th23", i, j);
145  calc->SetAngle(2, 3, th23);
146  double dmsq32 = GetBinValue(surface_file, "dmsq32", i, j);
147  calc->SetDm(3, calc->GetDm(2) + dmsq32);
148 
149  // Draw simulation after fit
150  TH1D* h_sim_after = sim->Predict(calc).ToTH1(pot);
151  h_sim_after->SetLineColor(kRed);
152  h_sim_after->Draw("hist same");
153 
154  // Draw fake data
155  TH1D* h_fake_data = fake_data.ToTH1(pot);
156  h_fake_data->Draw("same");
157 
158  TLegend l(0.55, 0.65, 0.85, 0.85);
159  l.AddEntry(h_sim_before, "Simulation before fit", "l");
160  l.AddEntry(h_sim_after, "Simulation after fit", "l");
161  l.AddEntry(h_fake_data, "Fake data", "p");
162  l.Draw();
163 
164  std::string plot_name = "spectrum_" + std::to_string(i) + "_" + std::to_string(j);
165  output_file->WriteTObject(c, plot_name.c_str());
166 
167  delete h_sim_before;
168  delete h_sim_after;
169  delete h_fake_data;
170  c->Clear();
171 
172  n_loop += 1;
173 
174  }
175  }
176 
177  p.Done();
178 
179 }
180 
181 #endif
const XML_Char * name
Definition: expat.h:151
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
enum BeamMode kRed
string output_name
Pickle.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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:176
double th23
Definition: runWimpSim.h:98
const char * p
Definition: xmltok.h:285
double GetBinValue(TFile *f, std::string name, int i, int j)
Definition: DrawSurface.C:31
Adapt the PMNS_Sterile calculator to standard interface.
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
std::vector< double > LogBins(int nbins, double xlo, double xhi)
Definition: DrawSurface.C:20
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
osc::OscCalcDumb calc
void FormatCanvas()
Canvas formatting utility.
void SetAngles(osc::OscCalcSterile *calc)
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const double kAna2018RHCPOT
Definition: Exposures.h:208
const int nbins
Definition: cellShifts.C:15
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
Sum up livetimes from individual cosmic triggers.
void DrawSurface(std::string file, std::string type, bool rhc=false, double dm41=-1, std::string detectors_in="both", int universe=1)
Definition: DrawSurface.C:40
#define pot
const double j
Definition: BetheBloch.cxx:29
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:16
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
double dmsq32
double livetime
Definition: saveFDMCHists.C:21
const double kAna2018SensitivityRHCNDPOT
Definition: Exposures.h:211
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
TFile * file
Definition: cellShifts.C:17
std::vector< TH1D * > GetCosmics(bool rhc, std::string detector)
Function to get default cosmics.
A simple ascii-art progress bar.
Definition: Progress.h:9
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
enum BeamMode kBlue
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
enum BeamMode string