DrawCCOscillations.C
Go to the documentation of this file.
1 /// DrawCCOscillations by J. Hewes <jhewes15@fnal.gov>
2 
3 // CAFAna headers
4 #include "CAFAna/Core/Progress.h"
7 
8 // ROOT headers
9 #include "TCanvas.h"
10 #include "TGraph.h"
11 #include "TMarker.h"
12 #include "TLegend.h"
13 
14 using namespace ana;
15 
17 
18  DontAddDirectory guard;
19  gErrorIgnoreLevel = kError;
20  std::vector<float> margin = { 0.15, 0.05 };
21 
22  auto calc = DefaultSterileCalc(4);
24  calc->SetL(1);
25 
26  // Flavours
27  std::vector<int> initFlav = { 12, 14 };
28  std::vector<int> finalFlav = { 12, 14, 16 };
29 
30  // Flavour names
31  std::map<int, std::string> flavString = { {12, "#nu_{e}"}, {14, "#nu_{#mu}"},
32  {16, "#nu_{#tau}"}, {-12, "#bar{#nu}_{e}"},
33  {-14, "#bar{#nu}_{#mu}"}, {-16, "#bar{#nu}_{#tau}"} };
34 
35 
36  Progress* p = new Progress("Drawing CC backgrounds");
37  size_t counter = 0;
38 
39  Binning xBins = Binning::Simple(50, 0, 1);
40  Binning yBins = Binning::LogUniform(50, 1e-2, 100);
41 
42  TH2D* hSurface = new TH2D(UniqueName().c_str(), ";sin^{2} #theta_{24};#Delta m_{41}^{2}",
43  xBins.NBins(), &xBins.Edges()[0], yBins.NBins(), &yBins.Edges()[0]);
44 
45  std::vector<float> th34Vals = { 0., 20., 45. };
46  std::vector<size_t> colours = { kRed-6, kGreen-5, kBlue-6 };
47 
48  // Loop over flavours
49  size_t nPoints = 10000;
50  for (int init : initFlav) {
51  for (int final : finalFlav) {
52  for (int sign = 0; sign <= 1; ++sign) {
53  int i = sign? init : -1 * init;
54  int f = sign? final : -1 * final;
55 
56  // Get the three-flavour oscillation probability
58  std::vector<float> x3f(nPoints);
59  std::vector<float> y3f(nPoints);
60  for (size_t it = 1; it <= nPoints; ++it) {
61  float energy = (20./(float)nPoints) * (float)it;
62  x3f[it] = energy;
63  y3f[it] = calc->P(i, f, energy);
64  }
65  TGraph g3f(x3f.size(), &x3f[0], &y3f[0]);
66  std::ostringstream title;
67  title << "P(" << flavString[i] << " #rightarrow " << flavString[f] << ")";
68  g3f.SetTitle(title.str().c_str());
69  g3f.GetXaxis()->SetTitle("True E (GeV)");
70  g3f.GetYaxis()->SetTitle("Osc. prob.");
71  g3f.SetLineWidth(3);
72  g3f.SetLineColor(kGray);
73 
74  // Loop over surface points
75  for (int yVal = -2; yVal < 3; ++yVal) {
76  double dm = pow(10, yVal);
78  std::vector<std::string> images;
79  for (size_t xBin = 1; xBin <= (size_t)xBins.NBins(); ++xBin) {
80  double ssth24 = hSurface->GetXaxis()->GetBinCenter(xBin);
82  std::vector<float> x(nPoints);
83  std::vector<std::vector<float>> y(th34Vals.size(), std::vector<float>(nPoints));
84  std::vector<std::vector<float>> yDiff(th34Vals.size(), std::vector<float>(nPoints));
85  for (size_t it = 1; it <= nPoints; ++it) {
86  float energy = (20./(float)nPoints) * (float)it;
87  x[it] = energy;
88  for (size_t itTh34 = 0; itTh34 < th34Vals.size(); ++itTh34) {
89  kFitTheta34InDegreesSterile.SetValue(calc, th34Vals[itTh34]);
90  y[itTh34][it] = calc->P(i, f, energy);
91  yDiff[itTh34][it] = y[itTh34][it] - y3f[it];
92  }
93  }
94  TCanvas c("c", "c", 1600, 900);
95  c.Divide(3);
96  // Draw surface pad
97  c.cd(1);
98  gPad->SetPad(0.01, 0.01, 0.49, 0.99);
99  gPad->SetLogy();
100  hSurface->Draw("colz");
101  TMarker m(ssth24, dm, 29);
102  m.SetMarkerSize(3);
103  m.Draw("same");
104  // Draw oscillation plot
105  c.cd(2);
106  gPad->SetPad(0.51, 0.51, 0.99, 0.99);
107  g3f.Draw("al");
108  std::vector<TGraph*> g;
109  TLegend l1(0.6, 0.15, 0.85, 0.4);
110  for (size_t it = 0; it < y.size(); ++it) {
111  g.push_back(new TGraph(x.size(), &x[0], &y[it][0]));
112  g[it]->SetLineWidth(3);
113  g[it]->SetLineColor(colours[it]);
114  g[it]->Draw("l same");
115  std::ostringstream name;
116  name << "#theta_{34} = " << th34Vals[it] << "#circ";
117  l1.AddEntry(g[it], name.str().c_str(), "l");
118  }
119  l1.Draw();
120  // Draw ratio plot
121  c.cd(3);
122  gPad->SetPad(0.51, 0.01, 0.99, 0.49);
123  std::vector<TGraph*> gDiff;
124  TLegend l2(0.6, 0.15, 0.85, 0.4);
125  for (size_t it = 0; it < yDiff.size(); ++it) {
126  gDiff.push_back(new TGraph(x.size(), &x[0], &yDiff[it][0]));
127  gDiff[it]->SetLineWidth(3);
128  gDiff[it]->SetLineColor(colours[it]);
129  gDiff[it]->SetTitle("");
130  gDiff[it]->GetXaxis()->SetTitle("True E (GeV)");
131  gDiff[it]->GetYaxis()->SetTitle("P(4flav) - P(3flav)");
132  std::string drawOpts = it == 0 ? "al" : "l same";
133  gDiff[it]->Draw(drawOpts.c_str());
134  std::ostringstream name;
135  name << "#theta_{34} = " << th34Vals[it] << "#circ";
136  l2.AddEntry(gDiff[it], name.str().c_str(), "l");
137  }
138  l2.Draw();
139  // Save plot
140  std::ostringstream fileName;
141  fileName << "plots/cc_oscprob/prob_flav" << i << "to" << f
142  << "_dm" << dm << "_x" << xBin << ".png";
143  c.cd(0);
144  c.SaveAs(fileName.str().c_str());
145  images.push_back(fileName.str());
146  c.Clear();
147  for (TGraph* tmp : g) delete tmp;
148  for (TGraph* tmp : gDiff) delete tmp;
149  } // for ssth24 value
150  std::ostringstream gifName;
151  gifName << "plots/cc_oscprob/prob_flav" << i << "to" << f
152  << "_dm" << dm << ".gif";
153  MakeGif(images, gifName.str(), 0.1);
154  for (std::string pngFile : images) {
155  std::ostringstream delCmd;
156  delCmd << "rm " << pngFile;
157  gSystem->Exec(delCmd.str().c_str());
158  }
159  ++counter;
160  p->SetProgress(float(counter)/float(initFlav.size()*finalFlav.size()*10));
161  } // for dm41 value
162  } // for sign
163  } // for final flavour
164  } // for initial flavour
165 
166  p->Done();
167 
168 } // Main macro function
const XML_Char * name
Definition: expat.h:151
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
fileName
Definition: plotROC.py:78
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
size_t nPoints
set< int >::iterator it
const char * p
Definition: xmltok.h:285
constexpr T pow(T x)
Definition: pow.h:75
Float_t tmp
Definition: plot.C:36
virtual double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
Definition: OscCalcDumb.h:21
osc::OscCalcDumb calc
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
double energy
Definition: plottest35.C:25
const std::vector< double > & Edges() const
Definition: Binning.h:34
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
static Binning LogUniform(int n, double lo, double hi)
Definition: Binning.cxx:134
int NBins() const
Definition: Binning.h:29
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
A simple ascii-art progress bar.
Definition: Progress.h:9
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const FitDmSq41Sterile kFitDmSq41Sterile
def sign(x)
Definition: canMan.py:197
void DrawCCOscillations()