plot_estimate_energy.C
Go to the documentation of this file.
1 //----------------------------------------
2 // Making Energy Estimator Plots for
3 // Inclusive numu CC xsection Analysis.
4 //
5 // Author : Biswaranjan Behera
6 // email : bbehera@fnal.gov
7 //----------------------------------------
8 
10 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Core/SystShifts.h"
12 
13 #include <vector>
14 #include <string>
15 #include <memory>
16 
17 #include "TH1.h"
18 #include "TH2.h"
19 #include "TFile.h"
20 #include "TCanvas.h"
21 #include "TF1.h"
22 #include "TGraph.h"
23 #include "TGraphAsymmErrors.h"
24 #include "TProfile.h"
25 
26 #include "Utilities/rootlogon.C"
27 
28 using namespace ana;
29 
31 {
35 };
36 
37 TF1 * linposfit = new TF1("linposfit", "[0]+[1]*x", 0., 15.);
38 
39 void MuonFit(Spectrum *plot1, string zaxisname,
40  string label,
41  MuonFitMode fitmode,
42  string outDir=""){
43 
44  double pot = 11E20;
45  TH2* hist = plot1->ToTH2(pot);
46  hist->SetTitle("");
47  switch (fitmode)
48  {
49  case kMuonFitAct:
50  hist->SetTitle("Active Region Only");
51  hist->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
52  hist->GetYaxis()->SetTitle("True Muon Energy (GeV)");
53  break;
54  case kMuonFitActAndCat:
55  if (fitmode == kMuonFitActAndCat){
56  hist->SetTitle("Muons Entering Catcher");
57  hist->GetXaxis()->SetTitle("Reco Muon Track Length in Active (m)");
58  hist->GetYaxis()->SetTitle("True Active Region Muon Energy (GeV)");
59  hist->GetYaxis()->SetTitleSize(0.048);
60  break;
61  case kMuonFitCat:
62  hist->SetTitle("Muons Entering Catcher");
63  hist->GetXaxis()->SetTitle("Reco Muon Track Length in Catcher (m)");
64  hist->GetYaxis()->SetTitle("True Catcher Region Muon Energy (GeV)");
65  break;
66  default:
67  throw std::runtime_error("Invalid fit mode");
68  }
69  }
70 
71  hist->GetZaxis()->SetTitle(zaxisname.c_str());
72  static const int nn = 150; //Number of points on graph - should be equal to or smaller than number of bins in x
73  float xxx [nn];
74  float exxl[nn];
75  float exxh[nn];
76  float yyx [nn];
77  float eyxl[nn];
78  float eyxh[nn];
79 
80  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
81  // make projection in Y view
82  TH1D* py = hist->ProjectionY("py",iw,iw,"");
83  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
84  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
85 
86  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
87 
88  double totalEntries = py->GetEntries();
89  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
90 
91  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
92 
93  //Find the location of width at half max
94  int yyhighbin = maxbin+1;
95  int yylowbin = maxbin-1;
96 
97  // find upper bin of half max
98  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
99  yyhighbin = ihigh;
100  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
101  }
102 
103  // find lower bin of half max
104  for(int ilow = maxbin-1; ilow > 0; --ilow){
105  yylowbin = ilow;
106  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
107  }
108 
109  float multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
110  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
111  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
112 
113  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
114  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
115 
116  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
117 
118  TF1 *f1 = new TF1("customGaus","gaus",0,5);
119  py->Fit("customGaus","OQ","",lowSide,highSide); //O and Q mean do not plot result, quiet mode
120  double mean = f1->GetParameter(1);
121  double error = f1->GetParError(1);
122  double chi = f1->GetChisquare();
123  double ndf = f1->GetNDF();
124 
125  // if (ndf > 0) {chiValues->Fill(chi/ndf);}
126  //else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
127  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
128  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
129  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
130  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
131  //Setting the error on the mean reported by the fit as error in y
132  eyxl[iw-1] = error;
133  eyxh[iw-1] = error;
134 
135  }//End of loop over bins in 2D histogram
136  CenterTitles(hist);
137  // TProfile * profile = hist->ProfileX();
138 
139  //->GetZaxis()->CenterTitle();
140  //Making graph from points
141  TGraphAsymmErrors* customProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
142  // TH1D* customProfile = profile->ProjectionX();
143  customProfile->GetXaxis()->SetNoExponent(kTRUE);
144  customProfile->GetXaxis()->CenterTitle();
145  customProfile->GetYaxis()->CenterTitle();
146  switch (fitmode){
147  case kMuonFitAct:
148  customProfile->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
149  customProfile->GetYaxis()->SetTitle("True Muon Energy (GeV)");
150  customProfile->GetXaxis()->SetLimits(0,15);
151  customProfile->GetHistogram()->SetMaximum(5);
152  customProfile->GetHistogram()->SetMinimum(0);
153  // customProfile->SetMaximum(5);
154  // customProfile->SetMinimum(0);
155  break;
156  case kMuonFitActAndCat:
157  customProfile->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
158  customProfile->GetYaxis()->SetTitle("True Muon E - True Muon E in Catcher (GeV)");
159  customProfile->GetXaxis()->SetLimits(0,15);
160  customProfile->GetHistogram()->SetMaximum(5);
161  customProfile->GetHistogram()->SetMinimum(0);
162  // customProfile->SetMaximum(5);
163  // customProfile->SetMinimum(0);
164  customProfile->GetYaxis()->SetTitleSize(0.048);
165  break;
166  case kMuonFitCat:
167  customProfile->GetXaxis()->SetTitle("Reco Muon Track Length in Catcher (m)");
168  customProfile->GetYaxis()->SetTitle("True Muon Energy in Catcher (GeV)");
169  customProfile->GetXaxis()->SetLimits(0,3);
170  customProfile->GetHistogram()->SetMaximum(3);
171  customProfile->GetHistogram()->SetMinimum(0);
172  // customProfile->SetMaximum(3);
173  // customProfile->SetMinimum(0);
174  break;
175  default:
176  throw std::runtime_error("Invalid fit mode");
177  }
178 
179  //CenterTitles(customProfile);
180  customProfile->SetTitle("");
181 
182  //Overlays 2d plot and graph
183  TCanvas* c1 = new TCanvas("c1","c1");
184  customProfile->SetMarkerColor(18);
185  customProfile->SetLineColor(18);
186  hist->Draw("colz");
187  customProfile->Draw("P");
188  // customProfile->Draw("e same");
189  Simulation();
190  c1->SetRightMargin(0.15);
191  c1->SetLeftMargin(0.15);
192 
193  if (!outDir.empty()){
194  c1->SaveAs((outDir+"/muon_" + label + ".png").c_str());
195  c1->SaveAs((outDir+"/muon_" + label + ".pdf").c_str());
196  }
197 
198  //Overlays 2d loqz plot and graph
199  TCanvas* c2 = new TCanvas("c2","c2");
200  c2->SetLogz();
201  customProfile->SetMarkerColor(13);
202  customProfile->SetLineColor(13);
203  hist->Draw("colz");
204  customProfile->Draw("P");
205  // customProfile->Draw("e same");
206  Simulation();
207  c2->SetRightMargin(0.15);
208  c2->SetLeftMargin(0.15);
209  if (!outDir.empty()){
210  c2->SaveAs((outDir+"/muon_logz" + label + ".png").c_str());
211  c2->SaveAs((outDir+"/muon_logz" + label + ".pdf").c_str());
212  }
213 
214  //Just the graph
215  TCanvas* c3 = new TCanvas("c3","c3");
216  customProfile->SetMarkerColor(1);
217  customProfile->SetLineColor(1);
218  customProfile->Draw("AP");
219  // customProfile->Draw("e same");
220  c3->SetRightMargin(0.15);
221  switch (fitmode){
222  case kMuonFitAct:
223  customProfile->Fit("pol3", "","",0.2, 12.0);
224  break;
225 
226  case kMuonFitActAndCat:
227  customProfile->Fit("pol3", "","",3.75, 12.2);
228  break;
229 
230  case kMuonFitCat:
231  customProfile->Fit("pol1", "","",0.02, 2.3);
232  break;
233 
234  default:
235  throw std::runtime_error("Invalid fit mode");
236  }
237 
238  c3->SetLeftMargin(0.15);
239 
240  if (!outDir.empty()){
241  c3->SaveAs((outDir+ "/fit_" + label + ".png").c_str());
242  c3->SaveAs((outDir+ "/fit_" + label + ".pdf").c_str());
243  }
244 }
245 
246 void plot_estimate_energy(string infile, string outputDir = ".")
247 {
248  TFile *fin = TFile::Open(infile.c_str(),"READ");
249  //------------------------------------------
250  // Accessing Files Here
251  //------------------------------------------
252 
253  Spectrum* MuonE_hist_act = Spectrum::LoadFrom(fin, "MuonE_hist_act" ).release();
254  Spectrum* MuonE_hist_actAndCat = Spectrum::LoadFrom(fin, "MuonE_hist_actAndCat").release();
255  Spectrum* MuonE_hist_cat = Spectrum::LoadFrom(fin, "MuonE_hist_cat" ).release();
256 
257  double pot = 11E20;
258 
259  //------------------------------------------------------
260  // Plots for Energy Estimator
261  //------------------------------------------------------
262 
263  MuonFit(MuonE_hist_act, Form("Events/(%2.2E POT)",pot),
264  "act", kMuonFitAct,outputDir);
265 
266  MuonFit(MuonE_hist_actAndCat, Form("Events/(%2.2E POT)",pot),
267  "act_cat", kMuonFitActAndCat,outputDir);
268 
269  MuonFit(MuonE_hist_cat, Form("Events/(%2.2E POT)",pot),
270  "cat", kMuonFitCat,outputDir);
271 
272 
273  // HadEFit(HadE, Form("Events/(%2.2E POT)",pot),
274  // "hade","hade_Logz",
275  // "hade_fit", true);
276 
277 }// end of void
void Simulation()
Definition: tools.h:16
TString fin
Definition: Style.C:24
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1483
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
c2
Definition: demo5.py:33
string infile
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
#define pot
Float_t f1
void MuonFit(Spectrum *plot1, string zaxisname, string label, MuonFitMode fitmode, string outDir="")
TF1 * linposfit
c1
Definition: demo5.py:24
void plot_estimate_energy(string infile, string outputDir=".")