LOverEPlot.C
Go to the documentation of this file.
1 #include <vector>
2 
6 
7 #include "TCanvas.h"
8 #include "TColor.h"
9 #include "TF1.h"
10 #include "TFile.h"
11 #include "TGaxis.h"
12 #include "TH1.h"
13 #include "TH2D.h"
14 #include "TLatex.h"
15 #include "TLegend.h"
16 #include "TLine.h"
17 
18 using namespace ana;
19 
20 const double loemin = 0.5/120.; // 0.5 km / 120 GeV
21 const double loemax = 810./.25; // 810 km / 0.25 GeV
22 
23 // Get a vector of bin edges that are logarithmic
24 // Bins go from loemin to loemax
25 std::vector<double> LogXBins();
26 
27 // Calculate the halfway point on the log axis
28 double HalfAxis();
29 
30 // Default canvas with margins set appropriately for labels
32 
33 // Plot text unrelated to the TLegend on a canvas
34 TLatex* MiscText(double x, double y, double size, TString text);
35 
36 // Set some basic options on a histogram
37 void SetHistOptions(TH1* hist, int col, bool fill = false);
38 void SetHistOptions(TH2* hist);
39 
40 void LOverEPlot()
41 {
42  TH1::AddDirectory(0);
43 
44  std::string folder = "/nova/ana/steriles/bless/";
45  std::string eventsName = "LOverE.root";
46  std::string canFileName = "LOverEPlot.root";
47 
48  TFile *fEvents = new TFile((folder + eventsName).c_str(), "READ");
49 
50  // Histograms of NC flux*xsec at ND and FD
51  TH1* hnd = (TH1*)fEvents->Get("loend/NOvA-ND/loend_numu_allpar_tot_nc_NOvA-ND");
52  TH1* hfd = (TH1*)fEvents->Get("loefd/NOvA-FD/loefd_numu_allpar_tot_nc_NOvA-FD");
53 
54  // Histograms of L, used much later on
55  TH1* hndL = (TH1*)fEvents->Get("Lnd/NOvA-ND/Lnd_numu_allpar_tot_nc_NOvA-ND");
56  TH1* hfdL = (TH1*)fEvents->Get("Lfd/NOvA-FD/Lfd_numu_allpar_tot_nc_NOvA-FD");
57 
58  double maxvalND = hnd->GetBinContent(hnd->GetMaximumBin());
59  double maxvalFD = hfd->GetBinContent(hfd->GetMaximumBin());
60 
61  std::vector<double> bins = LogXBins();
62  int nbins = bins.size() - 1;
63 
64  std::string allLabels = ";"; // Title
65  allLabels += "L/E (km/GeV);"; // X Label
66  allLabels += "1 - P(#nu_{#mu} #rightarrow #nu_{s});"; // Y Label
67  allLabels += "Events (Arbitrary Scale)"; // Z Label
68 
69  TH2D* hNDFlux = new TH2D("hNDFlux", allLabels.c_str(),
70  nbins, &bins[0],
71  1, 0, 1.2);
72  TH2D* hFDFlux = new TH2D("hFDFlux", allLabels.c_str(),
73  nbins, &bins[0],
74  1, 0, 1.2);
75 
76  // Scale histograms to arbitrary units
77  // Scale so max value of both histograms is 1
78  for(int i = 1; i <= nbins; ++i) {
79  hNDFlux->SetBinContent(i,1, hnd->GetBinContent(i)/maxvalND);
80  hFDFlux->SetBinContent(i,1, hfd->GetBinContent(i)/maxvalFD);
81  }
82 
83  SetHistOptions(hNDFlux);
84  SetHistOptions(hFDFlux);
85 
86  // Create grayscale gradient for event histograms
87  Int_t nb = 500;
88  Double_t Stops[2] = {0.50, 0.90};
89  Double_t Red[2] = {0.90, 0.50};
90  Double_t Green[2] = {0.90, 0.50};
91  Double_t Blue[2] = {0.90, 0.50};
92  TColor::CreateGradientColorTable(2, Stops, Red,Green,Blue, nb);
93  hNDFlux->SetContour(nb);
94  hFDFlux->SetContour(nb);
95 
96  std::vector<double> dms = {0.05, 0.5, 5, 50};
97 
98  // Vector of histograms that will contain osc prob curves
99  std::vector<TH1D*> hs;
100 
101  // Set up oscillation calculator that uses default 3 flavor parameters
102  // Set some additional 4 flavor parameters
104  calc->SetAngle(2,4, 10.*M_PI/180.);
105  calc->SetAngle(3,4, 35.*M_PI/180.);
106 
107  for(int i_dm=0, n = dms.size(); i_dm < n; ++i_dm) {
108  hs.push_back(new TH1D(UniqueName().c_str(), "", nbins, &bins[0]));
109 
110  calc->SetDm(4, dms[i_dm]);
111 
112  for(int i = 1; i <= nbins; ++i) {
113  // Oscillations can be very rapid
114  // Calculate average P by creating 1000 trapezoids across the bin
115  // Each bin contributes P(xLo)/2 + P(xHi)/2
116  // Thus, each endpoint contributes P(x)/2;
117  // but each midpoint contributes P(x)
118  // To make an average, multiply each trapezoid by its width,
119  // then divide by the bin width
120  // Since the trapezoid width is the bin width over 1000,
121  // this amounts to dividing by 1000 at the end
122  // (Really, divide by 1001, consider the average of the end points)
123  double loelo = hs[i_dm]->GetXaxis()->GetBinLowEdge(i);
124  double loehi = hs[i_dm]->GetXaxis()->GetBinUpEdge(i);
125  double delta = loehi - loelo;
126  //double loe = hs[i_dm]->GetBinCenter(i);
127 
128  double L = ((hs[i_dm]->GetBinCenter(i) < HalfAxis()) ? 0.7 : 810.) ; // ND vs FD
129  calc->SetL(L);
130 
131  double avgP = 0.;
132 
133  for(double loe = loelo + delta/1000.; loe < loehi; loe += delta/1000.) {
134  double E = L/loe;
135  avgP += calc->P(14, 0, E);
136  }
137 
138  // Don't forget end points
139  avgP += calc->P(14, 0, L/loelo)/2.;
140  avgP += calc->P(14, 0, L/loehi)/2.;
141 
142  // Make this a proper average!
143  avgP /= 1001.;
144 
145  hs[i_dm]->SetBinContent(i, avgP);
146  }
147 
148  hs[i_dm]->GetYaxis()->SetRangeUser(0, 1.2);
149  }
150 
151  SetHistOptions(hs[0], kBlue);
152  SetHistOptions(hs[1], kRed);
153  SetHistOptions(hs[2], kGreen+2);
154  SetHistOptions(hs[3], kMagenta+2);
155 
156  hNDFlux->GetZaxis()->SetRangeUser(0.03,1);
157  hFDFlux->GetZaxis()->SetRangeUser(0.03,1);
158 
159  TCanvas* c = MakeLongCanvas("cLOverE");
160  c->SetLogx();
161 
162  hNDFlux->Draw("colz");
163  hFDFlux->Draw("colz same");
164 
165  // Three-flavor oscillation line
166  TLine* threeF = new TLine(loemin, 1, loemax, 1);
167  threeF->SetLineWidth(3);
168  threeF->Draw("same");
169 
170  /*for(const auto& h : hs) {
171  h->Draw("c same");
172  }*/
173  for(int i = 0; i <= 2; ++i) {
174  hs[i]->Draw("c same");
175  }
176 
177  // Dotted line "separating" ND from FD
178  TLine l;
179  l.SetLineStyle(7);
180  l.SetLineWidth(2);
181  l.DrawLine(HalfAxis(), 0, HalfAxis(), 1.2);
182 
183  TLegend* leg = new TLegend(0.14, 0.35, 0.27, 0.55);
184  leg->SetBorderSize(0);
185  leg->SetFillStyle(0);
186  leg->SetTextFont(42);
187  leg->SetTextSize(0.05);
188  leg->AddEntry(threeF, "3-Flavor Prob.","l");
189  leg->AddEntry(hs[0], "#Deltam^{2}_{41} = 0.05 eV^{2}","l");
190  leg->AddEntry(hs[1], "#Deltam^{2}_{41} = 0.5 eV^{2}","l");
191  leg->AddEntry(hs[2], "#Deltam^{2}_{41} = 5 eV^{2}","l");
192  //leg->AddEntry(hs[3], "#Deltam^{2}_{41} = 50 eV^{2}","l");
193  leg->SetY1(leg->GetY2() - leg->GetNRows()*0.065);
194  leg->Draw();
195 
196  MiscText(0.25, 0.8, 0.06, "ND");
197  MiscText(0.65, 0.8, 0.06, "FD");
198 
199  gPad->SetTicks(0, 1);
200  //gPad->SetTicks(0, 2); // Draw Y-axis tick labels on right side
201 
202  c->RedrawAxis();
203 
204  // Draw energy axis on the top of the canvas
205  // Get bounds for drawing axes
206  double ymax = hNDFlux->GetYaxis()->GetXmax();
207  double ndL = hndL->GetMean();
208  double fdL = hfdL->GetMean();
209  std::cout << "ND Mean L: " << ndL << std::endl;
210  std::cout << "FD Mean L: " << fdL << std::endl;
211 
212  // ND half
213  TF1* f1 = new TF1("f1", "log10(x)", ndL/HalfAxis(), ndL/loemin);
214  TGaxis* gaxis = new TGaxis(HalfAxis(), ymax, loemin, ymax, "f1", 505, "+GS");
215  gaxis->SetLabelOffset(-0.04);
216  gaxis->SetTickSize(0.06);
217  gaxis->SetLabelSize(0.05);
218  gaxis->SetLabelFont(42);
219  gaxis->SetTitle("Neutrino Energy (GeV)");
220  gaxis->SetTitleFont(42);
221  gaxis->SetTitleOffset(1.15);
222  gaxis->SetTitleSize(0.05);
223  gaxis->CenterTitle();
224  gaxis->Draw();
225 
226  // FD half
227  TF1* f2 = new TF1("f2", "log10(x)", fdL/loemax, fdL/HalfAxis());
228  TGaxis* gaxis2 = new TGaxis(loemax, ymax, HalfAxis(), ymax, "f2", 505, "+GS");
229  gaxis2->SetLabelOffset(-0.04);
230  gaxis2->SetTickSize(0.06);
231  gaxis2->SetLabelSize(0.05);
232  gaxis2->SetLabelFont(42);
233  gaxis2->SetTitle("Neutrino Energy (GeV)");
234  gaxis2->SetTitleFont(42);
235  gaxis2->SetTitleOffset(1.15);
236  gaxis2->SetTitleSize(0.05);
237  gaxis2->CenterTitle();
238  gaxis2->Draw();
239 
240  TFile *fOut = new TFile((canFileName).c_str(), "RECREATE");
241  fOut->WriteTObject(c);
242  fOut->Close();
243  c->SaveAs("LOverE.png");
244 }
245 
246 //------------------------------------------------------------------------------
247 std::vector<double> LogXBins()
248 {
249  std::vector<double> bins;
250 
251  for(int i = 0; i <= 100; ++i) {
252  bins.push_back(
253  std::pow(10,
254  std::log10(loemin) +
256  )
257  );
258  }
259 
260  return bins;
261 }
262 
263 //------------------------------------------------------------------------------
264 double HalfAxis()
265 {
266  double minexp = std::log10(loemin);
267  double maxexp = std::log10(loemax);
268 
269  double midexp = minexp + (maxexp - minexp)/2.;
270  return std::pow(10, midexp);
271 }
272 
273 //------------------------------------------------------------------------------
275 {
276  TCanvas *c = new TCanvas(name.c_str(), name.c_str(), 1000, 600);
277  c->SetBottomMargin(0.12);
278  c->SetLeftMargin (0.11);
279  c->SetRightMargin (0.14);
280  c->SetTopMargin (0.12);
281 
282  return c;
283 }
284 
285 //------------------------------------------------------------------------------
286 TLatex* MiscText(double x, double y, double size, TString text)
287 {
288  TLatex *l = new TLatex(x, y, text);
289  l->SetNDC();
290  l->SetTextSize(size);
291  l->SetTextColor(kBlack);
292  l->Draw();
293  return l;
294 }
295 
296 //------------------------------------------------------------------------------
298 {
299  hist->GetXaxis()->CenterTitle();
300  hist->GetYaxis()->CenterTitle();
301  hist->GetZaxis()->CenterTitle();
302 
303  hist->SetLabelSize(0.05,"X");
304  hist->SetLabelSize(0.05,"Y");
305  hist->SetLabelSize(0.05,"Z");
306 
307  hist->GetXaxis()->SetTitleOffset(1.1);
308  hist->GetYaxis()->SetTitleOffset(1.);
309  hist->GetZaxis()->SetTitleOffset(0.8);
310 
311  hist->GetXaxis()->SetTitleSize(0.05);
312  hist->GetYaxis()->SetTitleSize(0.05);
313  hist->GetZaxis()->SetTitleSize(0.05);
314 
315  return;
316 }
317 
318 //------------------------------------------------------------------------------
319 void SetHistOptions(TH1* hist, int col, bool fill)
320 {
321  hist->GetXaxis()->CenterTitle();
322  hist->GetYaxis()->CenterTitle();
323 
324  hist->SetFillColor(fill*col);
325 
326  hist->SetLabelSize(0.05,"X");
327  hist->SetLabelSize(0.05,"Y");
328 
329  hist->SetLineColor(col);
330  hist->SetLineWidth(3);
331 
332  hist->SetMarkerColor(col);
333  hist->SetMarkerSize(0);
334 
335  hist->GetXaxis()->SetTitleOffset(1.1);
336  hist->GetYaxis()->SetTitleOffset(1);
337 
338  hist->GetXaxis()->SetTitleSize(0.05);
339  hist->GetYaxis()->SetTitleSize(0.05);
340 
341  return;
342 }
void SetHistOptions(TH1 *h, double max, std::string title, int ndiv, Color_t col, bool fill)
Set common options for a TLegend.
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double HalfAxis()
Definition: LOverEPlot.C:264
double delta
Definition: runWimpSim.h:98
std::vector< double > LogXBins()
Definition: LOverEPlot.C:247
constexpr T pow(T x)
Definition: pow.h:75
Adapt the PMNS_Sterile calculator to standard interface.
Float_t f2
osc::OscCalcDumb calc
#define M_PI
Definition: SbMath.h:34
Double_t ymax
Definition: plot.C:25
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
static constexpr double L
Float_t E
Definition: plot.C:20
Int_t col[ntarg]
Definition: Style.C:29
Float_t f1
virtual double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
TLatex * MiscText(double x, double y, double size, TString text)
Definition: LOverEPlot.C:286
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
void SetAngle(int i, int j, double th)
const double loemax
Definition: LOverEPlot.C:21
void SetDm(int i, double dm)
virtual void SetL(double L) override
T log10(T number)
Definition: d0nt_math.hpp:120
const double loemin
Definition: LOverEPlot.C:20
void LOverEPlot()
Definition: LOverEPlot.C:40
void fill(std::vector< T > &x, const S &y)
Definition: fill.hpp:22
static const double nb
Definition: Units.h:89
enum BeamMode kGreen
enum BeamMode kBlue
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TCanvas * MakeLongCanvas(std::string name)
Definition: LOverEPlot.C:274
enum BeamMode string