plot_DataMCComp_numu.C
Go to the documentation of this file.
1 //------------------------------------------------------------------------------
2 // This plots the noextrap predictions made by make_DataMCComp_numu.C.
3 // The map of variables are defined at the bottom of the header.
4 // FHC pot-equiv is put in by hand at line 145
5 //
6 // The fakeData bool should be false if you are using real data. I just had that
7 // there for testing while the data files were still blinded
8 //
9 // Andrew Sutton - ats2mk@virginia.edu
10 //------------------------------------------------------------------------------
11 
12 #pragma once
13 
14 #include "DataMCComp_header.h"
15 
16 using namespace ana;
17 
18 TH1D* Load1DHistFromPred(TFile* inFile, std::string name, PredSel predsel, double pot, osc::IOscCalc *calc);
19 
20 TH1D* LoadFakeDataFromPred(TFile* inFile, std::string name, double pot);
21 
22 TH1D* Load1DHistFromSpec(TFile* inFile, std::string name, bool cosmic, double pot, double lt);
23 
24 void plot_DataMCComp_numu(std::string inName, bool fakeData = false)
25 {
26  TFile *inFile = new TFile(inName.c_str(), "READ");
27  //TFile *dataFile = new TFile("dataonly_fhc_period1235.root", "READ");
28  // TFile *dataFile = new TFile("dataonly_fhc_period910.root", "READ");
29 
30  Sign::Sign_t WS;
31  std::string horn, HORN;
32  double pot, lt;
33 
34  if (inName.find("fhc") != std::string::npos) {
35  WS = Sign::kAntiNu;
36  horn = "Neutrino Beam";
37  HORN = "FHC";
38  pot = kAna2020FHCPOT;
40  }
41  else if (inName.find("rhc") != std::string::npos) {
42  WS = Sign::kNu;
43  HORN = "RHC";
44  horn = "Antineutrino Beam";
45  pot = kAna2020RHCPOT;
46  lt = kAna2020RHCLivetime;
47  }
48  else std::abort();
49 
50  // get the BF
51  std::string sHieOct = "NH_UO";
52  std::string sFitDir = "/nova/ana/3flavor/Ana2020/Fits/BestFits/";
53  std::string sFitFile = sFitDir + "bestfits_joint_realData_both_systs.root";
54  TFile *fBestFit = new TFile(sFitFile.c_str(), "READ");
55  auto calc = (osc::IOscCalcAdjustable*)(LoadFrom<osc::IOscCalc>(fBestFit, (sHieOct + "_calc").c_str())).release();
56 
57  // Get POT and LT from the datafile
58  double dpot = pot;
59  double dlt = lt;
60  if (inFile->GetDirectory("dtSpec_muonE_AllQuants")) {
61  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom(inFile, "dtSpec_muonE_AllQuants");
62  dpot = TempSpec->POT();
63  dlt = TempSpec->Livetime();
64  }
65 
66  // Quantiles and desired
67  const std::vector<std::string> quantNames =
68  { "AllQuants", "Quant1", "Quant2", "Quant3", "Quant4"};
69 
70  std::string WSH = "WS" + HORN;
71 
72  std::vector< TH1D* > hPred[5];
73  std::vector< TH1D* > hBkg[5];
74  std::vector< TH1D* > hNumu[5];
75  std::vector< TH1D* > hWS[5];
76  std::vector< TH1D* > hCosm[5];
77  std::vector< TH1D* > hData[5];
78  std::vector< TH1D* > hDataUS[5];
79 
80  for (uint q = 0; q < quantNames.size(); ++q) {
81  for (auto var : vars) {
82 
83  std::string varName = var.first.first;
84 
85  hPred[q].push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at("All") , pot, calc) );
86  hBkg[q] .push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at("All") , pot, calc) );
87  hNumu[q].push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at("Numu"), pot, calc) );
88  hWS[q] .push_back( Load1DHistFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], predSels.at(WSH) , pot, calc) );
89 
90  hCosm[q].push_back( Load1DHistFromSpec(inFile, "csSpec_" + varName + "_" + quantNames[q], true , pot, lt) );
91 
92  if (fakeData) {
93  hData [q].push_back( LoadFakeDataFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], pot) );
94  hDataUS[q].push_back( LoadFakeDataFromPred(inFile, "mcPred_" + varName + "_" + quantNames[q], pot) );
95  }
96  else {
97  hData [q].push_back( Load1DHistFromSpec(inFile, "dtSpec_" + varName + "_" + quantNames[q], false, dpot, dlt) );
98  hDataUS[q].push_back( Load1DHistFromSpec(inFile, "dtSpec_" + varName + "_" + quantNames[q], false, dpot, dlt) );
99  }
100 
101  // adding together for custom hist stack
102  hBkg [q].back()->Add(hNumu[q].back(), -1); // remove the signal
103  hBkg [q].back()->Add(hCosm[q].back(), +1); // add on the cosmics for stack
104  hWS [q].back()->Add(hBkg [q].back(), +1); // add bkg to WS for stack
105  hPred[q].back()->Add(hCosm[q].back(), +1); // add cosmics to full prediction
106 
107  // Beautification
108  PimpHist(hPred [q].back(), kSolid, kFullPredColor , 3, kFullCircle, kFullPredColor , 1);
109  PimpHist(hWS [q].back(), kSolid, kNumuWSColor , 1, kFullCircle, kNumuWSColor , 1);
110  PimpHist(hBkg [q].back(), kSolid, kTotalBkgColor , 1, kFullCircle, kTotalBkgColor , 1);
111  PimpHist(hCosm [q].back(), kSolid, kCosmicBackgroundColor, 1, kFullCircle, kCosmicBackgroundColor, 1);
112 
113  FillWithDimColor(hWS[q].back() , 0);
114  FillWithDimColor(hCosm[q].back(), 0);
115  FillWithDimColor(hBkg[q].back() , 0);
116 
117  if (varName == "recoE"){
118  hData[q].back()->Scale(0.1, "width");
119  hPred[q].back()->Scale(0.1, "width");
120  hCosm[q].back()->Scale(0.1, "width");
121  hBkg [q].back()->Scale(0.1, "width");
122  hWS [q].back()->Scale(0.1, "width");
123  hNumu[q].back()->Scale(0.1, "width");
124  }
125  }// end var loop
126  }// end quant loop
127 
128  // Plotting
129  int v = 0;
130  for (auto var: vars) {
131  std::ofstream f(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.txt").c_str());
132  std::ofstream f_ByQ(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.txt").c_str());
133  TCanvas *c = new TCanvas("c", "c", 960, 800);
134  TCanvas *c_ByQ;
135  TPad *p1, *p2, *p3, *p4;
136  SplitCanvasQuant(c_ByQ, p1, p2, p3, p4);
137  c_ByQ->SetTitle("c_ByQ");
138 
139  TLatex* xTitle = new TLatex(0.5, 0.03, var.first.second.c_str());
140  xTitle->SetNDC();
141  xTitle->SetTextAlign(22);
142  TString eventAppend = Form("%.2f #times 10^{20} POT", dpot/1e20);
143  if (HORN == "FHC") {
144  // put in by hand for equiv
145  dpot = 13.6e20;
146  eventAppend = Form("%.2f #times 10^{20} POT-equiv.", dpot/1e20);
147  }
148  if (var.first.first == "recoE")
149  eventAppend = "0.1 GeV";
150  TLatex* yTitle = new TLatex(0.023, 0.5, "Events / " + eventAppend);
151  yTitle->SetNDC();
152  yTitle->SetTextAlign(22);
153  yTitle->SetTextAngle(90.);
154 
155  // setup the legends
156  TLegend *leg = new TLegend(0.60,0.60,0.88,0.88);
157  TLegend *leg_l = new TLegend(0.10,0.75,0.25,0.88);
158  TLegend *leg_r = new TLegend(0.52,0.75,0.65,0.88);
159  leg->AddEntry(hData[0][v], "FD data" ,"LEP"); leg_l->AddEntry(hData[0][v], "FD data" ,"LEP");
160  leg->AddEntry(hPred[0][v], "2020 Best-fit" ,"L" ); leg_l->AddEntry(hPred[0][v], "2020 Best-fit" ,"L" );
161  leg->AddEntry(hBkg [0][v], "Background" ,"BF" ); leg_l->AddEntry(hBkg [0][v], "Background" ,"BF" );
162  leg ->SetTextSize (0.05); leg ->SetBorderSize(0); leg ->SetFillStyle (0);
163  leg_l->SetTextSize (0.03); leg_l->SetBorderSize(0); leg_l->SetFillStyle (0);
164  leg_r->SetTextSize (0.03); leg_r->SetBorderSize(0); leg_r->SetFillStyle (0);
165 
166  // Get the maximum quant bin for plotting
167  float maxDataByQ = 0;
168  float maxPredByQ = 0;
169  for (uint q = 1; q < quantNames.size(); ++q) {
170  float maxData = hData[q][v]->GetMaximum();
171  float maxPred = hPred[q][v]->GetMaximum();
172 
173  if (maxData > maxDataByQ) maxDataByQ = maxData;
174  if (maxPred > maxPredByQ) maxPredByQ = maxPred;
175  }
176 
177  // quantile loop
178  for (uint q = 0; q < quantNames.size(); ++q) {
179  // Making plots look better and changing pads
180  if (q == 0)
181  MakeHistCanvasReady_Quant(q, hPred[q][v], hData[q][v]->GetMaximum()*2.2);
182  else
183  MakeHistCanvasReady_Quant(q, hPred[q][v], maxDataByQ*1.8);
184 
185  if (q == 0) c ->cd();
186  if (q == 1) p1->cd();
187  if (q == 2) p2->cd();
188  if (q == 3) p3->cd();
189  if (q == 4) p4->cd();
190 
191  // Set the data error bars and draw everything
192  TGraph * grData = GraphWithPoissonErrors2(hData[q][v], hPred[q][v], false, true);
193  if (var.first.first == "recoE") grData = graphAsymmErrorScaled(hData[q][v], hDataUS[q][v]);
194  grData->SetMarkerStyle(kFullCircle);
195 
196  hPred[q][v]->Draw("hist same");
197  hBkg[q][v] ->Draw("hist same");
198  if ( !(var.first.first == "trueE" || var.first.first == "resE") )
199  grData->Draw("P");
200 
201  DrawQuantLabel(q);
202  }// end quant loop
203 
204  c->cd();
205  leg->Draw();
206  xTitle->Draw();
207  yTitle->Draw();
208  Preliminary();
209  DrawBeamLabel(HORN=="FHC");
210  c->Update();
211  c->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.pdf").c_str());
212  c->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.png").c_str());
213  c->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_AllQuants.eps").c_str());
214  f << HORN << " " << var.first.second
215  << " with an un-extrapolated prediction set to the 2020 joint analysis best fit point"
216  << std::endl;
217 
218  c_ByQ->cd();
219  Preliminary();
220  DrawBeamLabel(HORN=="FHC");
221  xTitle->SetTextSize(0.04);
222  yTitle->SetTextSize(0.04);
223  xTitle->Draw();
224  yTitle->Draw();
225  p1->cd();
226  leg_l->Draw();
227  p2->cd();
228  leg_r->Draw();
229  c_ByQ->Update();
230  c_ByQ->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.pdf").c_str());
231  c_ByQ->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.png").c_str());
232  c_ByQ->SaveAs(("dataMC_" + HORN + "_" + var.first.first + "_ByQuant.eps").c_str());
233  f_ByQ << HORN << " " << var.first.second
234  << " split by hadronic energy fraction quartile with an"
235  << " un-extrapolated prediction set to the 2020 joint analysis best fit point"
236  << std::endl;
237 
238  f.close();
239  f_ByQ.close();
240 
241  delete c;
242  delete c_ByQ;
243 
244  ++v;
245  }// end var loop
246 }
247 
248 
249 //*******************************************************************************
251 {
252  // -- Here you can change the best fit point
253  // osc::OscCalcDumb* calc = new osc::OscCalcDumb();
254  // osc::IOscCalcAdjustable* calc = DefaultOscCalc();
255  // calc->SetdCP (2*M_PI);
256  // calc->SetTh23(asin(sqrt(0.565)));
257  // calc->SetDmsq32(2.48e-3);
258 
259  std::unique_ptr<PredictionNoExtrap> TempPred =
260  PredictionNoExtrap::LoadFrom(inFile, (name).c_str());
261 
262  return TempPred->PredictComponent(calc, predsel.flav, predsel.curr, predsel.sign).ToTH1(pot);
263 }
264 
265 //*******************************************************************************
267 {
269 
270  std::unique_ptr<PredictionNoExtrap> TempPred =
271  PredictionNoExtrap::LoadFrom(inFile, (name).c_str());
272 
273  return TempPred->Predict(calc).MockData(pot).ToTH1(pot);
274 
275 }
276 
277 //*******************************************************************************
278 TH1D* Load1DHistFromSpec(TFile* inFile, std::string name, bool cosmic, double pot, double lt)
279 {
280  std::unique_ptr<Spectrum> TempSpec =
281  Spectrum::LoadFrom(inFile, (name).c_str());
282 
283  if(cosmic)
284  return TempSpec->ToTH1(lt, kLivetime);
285  else
286  return TempSpec->ToTH1(pot);
287 }
TH1D * LoadFakeDataFromPred(TFile *inFile, std::string name, double pot)
void plot_DataMCComp_numu(std::string inName, bool fakeData=false)
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double lt
Antineutrinos-only.
Definition: IPrediction.h:50
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
Definition: Plots.cxx:1365
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
const std::map< std::string, PredSel > predSels
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
Flavors::Flavors_t flav
TPaveText * DrawQuantLabel(int quant)
Definition: Plots.cxx:1576
osc::OscCalcDumb calc
Current::Current_t curr
void PimpHist(TH1 *hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize)
Pimp histogram once and for all.
Definition: Plots.cxx:1408
const Color_t kFullPredColor
Definition: Style.h:36
ifstream inFile
Definition: AnaPlotMaker.h:34
Sum up livetimes from individual cosmic triggers.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:980
TPaveText * DrawBeamLabel(bool isFHC)
Put the standardized beam label in the left corner of the active canvas.
Definition: Plots.cxx:1557
TGraphAsymmErrors * GraphWithPoissonErrors2(const TH1 *h, const TH1 *h2, bool noErrorsXaxis, bool drawEmptyBins)
Same as above but use a reference histogram to determine which empty bins to draw.
Definition: Plots.cxx:944
#define pot
TH1D * Load1DHistFromSpec(TFile *inFile, std::string name, bool cosmic, double pot, double lt)
const double kAna2020FHCLivetime
Definition: Exposures.h:237
const Color_t kTotalBkgColor
Definition: Style.h:39
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
const std::map< std::pair< std::string, std::string >, Variable > vars
const double kAna2020FHCPOT
Definition: Exposures.h:233
const double kAna2020RHCPOT
Definition: Exposures.h:235
void Preliminary()
Neutrinos-only.
Definition: IPrediction.h:49
const double kAna2020RHCLivetime
Definition: Exposures.h:238
const Color_t kNumuWSColor
Definition: Style.h:38
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1439
Sign::Sign_t sign
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
TH1D * Load1DHistFromPred(TFile *inFile, std::string name, PredSel predsel, double pot, osc::IOscCalc *calc)
enum BeamMode string
unsigned int uint