Plotting_Data2DataComp_SingleCanvas.C
Go to the documentation of this file.
1 /*******************************************************************
2  * A SCRIPT TO READ THE OUTPUT ROOT FILES OF Analyse_Data2DataComp.C,
3  * WILL PLOT SPECTRA FROM THREE DATA SETS AND THE THREE CORRESPONDING
4  * RATIOS ON ON SEPARATE PLOTS ON THE SAME CANVAS. IMPORTANT
5  * PARAMETERS ADDED TO EACH CANVAS ALSO.
6  * OUTPUT CSV FILE CREATED WITH KEY QUANTITIES SUCH AS THE FRACTIONAL
7  * CHANGES IN THE MEANS ETC.
8  * USAGE:
9  * Plotting_Data2DataComp_SingleCanvas.C <Directory to output pdfs>
10  * <String to append outfile names>
11  * <Input root file containing 2 of 3 datasets>
12  * <Input root file containing third dataset>
13  ********************************************************************/
14 #include "CAFAna/Core/Binning.h"
15 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/Spectrum.h"
17 #include "CAFAna/Core/Ratio.h"
18 #include "CAFAna/Analysis/Plots.h"
19 
20 #include <fstream>
21 #include "TCanvas.h"
22 #include "TFile.h"
23 #include "TGraph.h"
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "TH2D.h"
27 #include "TMath.h"
28 #include "TGaxis.h"
29 #include "TMultiGraph.h"
30 #include "TLegend.h"
31 #include "TLatex.h"
32 #include "TStyle.h"
33 #include "THStack.h"
34 #include "TPaveText.h"
35 #include "TLine.h"
36 
37 
38 using namespace ana;
39 
40 double getFractionalError(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
41 {
42  double fractionalError = 0.;
43 
44  double numerator = mean1 - mean0;
45  double denominator = mean0;
46  double numeratorVar = meanErr0*meanErr0 + meanErr1*meanErr1;
47  double denominatorVar = meanErr0*meanErr0;
48 
49  return std::sqrt((numerator/denominator)*(numerator/denominator)*(numeratorVar/(numerator*numerator)+denominatorVar/(denominator*denominator)))*100.;
50 }
51 
52 double is2SigmaBandsTouching(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
53 {
54  std::pair<double,double> largeMean;
55  std::pair<double,double> smallMean;
56 
57  if(mean0 > mean1)
58  {
59  largeMean = {mean0, meanErr0};
60  smallMean = {mean1, meanErr1};
61  }
62  else
63  {
64  largeMean = {mean1, meanErr1};
65  smallMean = {mean0, meanErr0};
66  }
67 
68  double largeMeanShifted = largeMean.first - 2.*largeMean.second;
69  double smallMeanShifted = smallMean.first + 2.*smallMean.second;
70 
71  return (largeMeanShifted <= smallMeanShifted ? 1 : 0);
72 }
73 
74 void addToKeyPlotList(std::string const &s_OutDir, std::string const &s_FileName)
75 {
76  std::ofstream f_KeyPlots;
77  f_KeyPlots.open(s_OutDir+"KeyPlots.txt", std::ios_base::app);
78  f_KeyPlots << s_FileName << std::endl;
79  f_KeyPlots.close();
80 
81  return;
82 }
83 
84 void makeCSVFile(std::string const &s_OutDir, std::string const &s_FileAppend, std::string const &s_FileName, std::string const &s_Title,
85  unsigned int const &quantile, std::vector<double> const &vec_Means)
86 {
87  std::ofstream f_CSV;
88  f_CSV.open(s_OutDir+"MeanFile_"+s_FileAppend+".csv", std::ios_base::app);
89  f_CSV << s_FileName << "=" << s_Title << "=" << quantile << "=";
90  for(unsigned int i = 0; i < vec_Means.size(); i++)
91  {
92  f_CSV << vec_Means.at(i);
93  if(i != vec_Means.size()-1)
94  {
95  f_CSV << "=";
96  }
97  }
98  f_CSV << std::endl;
99  f_CSV.close();
100 
101  return;
102 }
103 
104 //THE ORDER SHOULD BE 7d vs. 8b, 2018 vs 7d.
105 void Plotting_Data2DataComp_SingleCanvas(std::string s_OutDir = "", std::string s_FileAppend = "",
106  std::string s_InFile1 = "", std::string s_InFile2 = "")
107 {
108  gStyle->SetOptStat(0);
109 
110  TFile f_InFile1(s_InFile1.c_str());
111  TFile f_InFile2(s_InFile2.c_str());
112  char s_DirName [256];
113  char s_FileName[256];
114  unsigned int nQuant = 4;
115  unsigned int nSpec = 41;
116 
117  for(unsigned int S = 0; S <= nSpec; ++S) {
118  for(unsigned int Q = 0; Q <= nQuant; ++Q) {
119  sprintf(s_DirName, "dir_Var%d_Q%d_DataSet0",S,Q);
120  std::unique_ptr<Spectrum> f1_spec0 = Spectrum::LoadFrom(&f_InFile1, s_DirName); //7d
121  std::unique_ptr<Spectrum> f2_spec0 = Spectrum::LoadFrom(&f_InFile2, s_DirName); //2018
122  sprintf(s_DirName, "dir_Var%d_Q%d_DataSet1",S,Q);
123  std::unique_ptr<Spectrum> f1_spec1 = Spectrum::LoadFrom(&f_InFile1, s_DirName); //8b.
124 
125  TCanvas *c = new TCanvas("c","c", 1600, 1600);
126  TPad *p_Left = new TPad("p_Left", "", 0.005, 0.005, 0.610, 0.995);
127  TPad *p_Right = new TPad("p_Right", "", 0.620, 0.005, 0.995, 0.995);
128 
129  TPad *p_LeftTop = new TPad("p_LeftTop", "", 0.0, 0.6, 1.0, 1.0);
130  TPad *p_LeftBottom1 = new TPad("p_LeftBottom1", "", 0.0, 0.4, 1.0, 0.6);
131  TPad *p_LeftBottom2 = new TPad("p_LeftBottom2", "", 0.0, 0.2, 1.0, 0.4);
132  TPad *p_LeftBottom3 = new TPad("p_LeftBottom3", "", 0.0, 0.0, 1.0, 0.2);
133 
134  TPad *p_RightTop = new TPad("p_RightTop", "", 0.0, 0.5, 1.0, 1.0);
135  TPad *p_RightBottom = new TPad("p_RightBottom", "", 0.0, 0.0, 1.0, 0.5);
136 
137  p_Left ->Draw();
138  p_Right->Draw();
139 
140  p_Left->cd();
141  p_LeftTop->SetTopMargin(0);
142  p_LeftTop->SetRightMargin(p_LeftTop->GetRightMargin()-0.09);
143  p_LeftTop->Draw();
144  p_LeftBottom1->Draw(); p_LeftBottom2->Draw(); p_LeftBottom3->Draw();
145 
146  p_Right->cd();
147  p_RightTop->Draw(); p_RightBottom->Draw();
148 
149  //EVERYTHING SCALED TO 7d.
150  TH1 *f1_h0 = f1_spec0->ToTH1(f1_spec0->POT()); f1_h0->SetLineColor(kBlue );
151  TH1 *f2_h0 = f2_spec0->ToTH1(f1_spec0->POT()); f2_h0->SetLineColor(kBlack);
152  TH1 *f1_h1 = f1_spec1->ToTH1(f1_spec0->POT()); f1_h1->SetLineColor(kRed );
153 
154  p_LeftTop->cd();
155  f1_h0->Draw();
156  f1_h1->Draw("SAME");
157  f2_h0->Draw("SAME");
158 
159  p_LeftBottom1->cd();
160  p_LeftBottom1->SetTopMargin (p_LeftBottom1->GetTopMargin() -0.08);
161  p_LeftBottom1->SetBottomMargin(p_LeftBottom1->GetBottomMargin()+0.08);
162  p_LeftBottom1->SetRightMargin (p_LeftBottom1->GetRightMargin()-0.09);
163  Ratio rat2(*f1_spec1, *f2_spec0);
164  TH1* hRat2 = rat2.ToTH1();
165  hRat2->SetLineColor(kViolet);
166  hRat2->GetXaxis()->SetTitle("8b / 4+6");
167  hRat2->GetXaxis()->SetTitleSize(0.13);
168  hRat2->GetXaxis()->SetTitleOffset(0.58);
169  hRat2->GetYaxis()->SetTitleSize(0.12);
170  hRat2->GetYaxis()->SetTitleOffset(0.39);
171  hRat2->GetXaxis()->SetLabelSize(0.08);
172  hRat2->GetYaxis()->SetLabelSize(0.08);
173  hRat2->SetAxisRange(0.9,1.1,"Y");
174  hRat2->Draw();
175 
176  double x[2] = {hRat2->GetXaxis()->GetXmin(), hRat2->GetXaxis()->GetXmax()};
177  double y[2] = {1.0, 1.0};
178  TGraph* l_Xaxis = new TGraph(2, x, y);
179  l_Xaxis->SetLineColor(kGray);
180  l_Xaxis->SetLineWidth(3);
181  l_Xaxis->SetLineStyle(2);
182  l_Xaxis->Draw("L SAME");
183 
184  p_LeftBottom2->cd();
185  p_LeftBottom2->SetTopMargin (p_LeftBottom2->GetTopMargin() -0.08);
186  p_LeftBottom2->SetBottomMargin(p_LeftBottom2->GetBottomMargin()+0.08);
187  p_LeftBottom2->SetRightMargin (p_LeftBottom2->GetRightMargin()-0.09);
188  Ratio rat1(*f1_spec1, *f1_spec0);
189  TH1* hRat1 = rat1.ToTH1();
190  hRat1->SetLineColor(kViolet);
191  hRat1->GetXaxis()->SetTitle("8b / 7d");
192  hRat1->GetXaxis()->SetTitleSize(0.13);
193  hRat1->GetXaxis()->SetTitleOffset(0.58);
194  hRat1->GetYaxis()->SetTitleSize(0.12);
195  hRat1->GetYaxis()->SetTitleOffset(0.39);
196  hRat1->GetXaxis()->SetLabelSize(0.08);
197  hRat1->GetYaxis()->SetLabelSize(0.08);
198  hRat1->SetAxisRange(0.9,1.1,"Y");
199  hRat1->Draw();
200  l_Xaxis->Draw("L SAME");
201 
202  p_LeftBottom3->cd();
203  p_LeftBottom3->SetTopMargin (p_LeftBottom3->GetTopMargin() -0.08);
204  p_LeftBottom3->SetBottomMargin(p_LeftBottom3->GetBottomMargin()+0.07);
205  p_LeftBottom3->SetRightMargin (p_LeftBottom3->GetRightMargin()-0.09);
206  Ratio rat3(*f1_spec0, *f2_spec0);
207  TH1* hRat3 = rat3.ToTH1();
208  hRat3->SetLineColor(kViolet);
209  hRat3->GetXaxis()->SetTitle("7d / 4+6");
210  hRat3->GetXaxis()->SetTitleSize(0.13);
211  hRat3->GetXaxis()->SetTitleOffset(0.58);
212  hRat3->GetYaxis()->SetTitleSize(0.12);
213  hRat3->GetYaxis()->SetTitleOffset(0.39);
214  hRat3->GetXaxis()->SetLabelSize(0.08);
215  hRat3->GetYaxis()->SetLabelSize(0.08);
216  hRat3->SetAxisRange(0.9,1.1,"Y");
217  hRat3->Draw();
218  l_Xaxis->Draw("L SAME");
219 
220  p_RightTop->cd();
221  TLegend *leg = new TLegend(0.1, 0.4, 0.9, 0.95);
222  leg->SetTextSize(0.08);
223  std::string s_LegTitle = f1_spec1->GetLabels()[0].substr(0, f1_spec1->GetLabels()[0].find('('));
224  leg->SetHeader(s_LegTitle.c_str(),"C");
225  leg->AddEntry(f1_h1,"8b","L");
226  leg->AddEntry(f1_h0,"7d","L");
227  leg->AddEntry(f2_h0,"4+6","L");
228  leg->Draw();
229  std::string s_Quantile = Q+1 == 5 ? "All" : std::to_string(Q+1);
230  TText *txt_Cut = new TText(0.25, 0.2, ("Cut: "+s_FileAppend).c_str());
231  TText *txt_Quant = new TText(0.25, 0.1, ("Quantile: "+s_Quantile).c_str());
232  txt_Cut->SetTextSize(0.09); txt_Quant->SetTextSize(0.09);
233  txt_Cut->Draw(); txt_Quant->Draw();
234 
235  std::vector<double> vec_Means;
236  p_RightBottom->cd();
237  TText *txt_8b = new TText(0.2, 0.9, "8b:");
238  TText *txt_8b_mean = new TText(0.2, 0.85, Form("Mean: %f +/- %f", f1_h1->GetMean(), f1_h1->GetMeanError()));
239  TText *txt_8b_rms = new TText(0.2, 0.8, Form("RMS: %f +/- %f", f1_h1->GetRMS(), f1_h1->GetRMSError() ));
240  TText *txt_8b_int = new TText(0.2, 0.75, Form("Integral: %f", f1_h1->Integral()));
241  txt_8b->SetTextColor(kRed); txt_8b_mean->SetTextColor(kRed); txt_8b_rms->SetTextColor(kRed); txt_8b_int->SetTextColor(kRed);
242  txt_8b->Draw(); txt_8b_mean->Draw(); txt_8b_rms->Draw(); txt_8b_int->Draw();
243  vec_Means.push_back(f1_h1->GetMean ()); vec_Means.push_back(f1_h1->GetMeanError());
244  vec_Means.push_back(f1_h1->GetRMS ()); vec_Means.push_back(f1_h1->GetRMSError ());
245  vec_Means.push_back(f1_h1->Integral());
246 
247 
248  TText *txt_7d = new TText(0.2, 0.65, "7d:");
249  TText *txt_7d_mean = new TText(0.2, 0.6, Form("Mean: %f +/- %f", f1_h0->GetMean(), f1_h0->GetMeanError()));
250  TText *txt_7d_rms = new TText(0.2, 0.55, Form("RMS: %f +/- %f", f1_h0->GetRMS(), f1_h0->GetRMSError() ));
251  TText *txt_7d_int = new TText(0.2, 0.5, Form("Integral: %f", f1_h0->Integral()));
252  txt_7d->SetTextColor(kBlue); txt_7d_mean->SetTextColor(kBlue); txt_7d_rms->SetTextColor(kBlue); txt_7d_int->SetTextColor(kBlue);
253  txt_7d->Draw(); txt_7d_mean->Draw(); txt_7d_rms->Draw(); txt_7d_int->Draw();
254  vec_Means.push_back(f1_h0->GetMean ()); vec_Means.push_back(f1_h0->GetMeanError());
255  vec_Means.push_back(f1_h0->GetRMS ()); vec_Means.push_back(f1_h0->GetRMSError ());
256  vec_Means.push_back(f1_h0->Integral());
257 
258  TText *txt_Ana2018 = new TText(0.2, 0.4, "4+6:");
259  TText *txt_Ana2018_mean = new TText(0.2, 0.35, Form("Mean: %f +/- %f", f2_h0->GetMean(), f2_h0->GetMeanError()));
260  TText *txt_Ana2018_rms = new TText(0.2, 0.3, Form("RMS: %f +/- %f", f2_h0->GetRMS(), f2_h0->GetRMSError() ));
261  TText *txt_Ana2018_int = new TText(0.2, 0.25, Form("Integral: %f", f2_h0->Integral()));
262  txt_Ana2018->SetTextColor(kBlack); txt_Ana2018_mean->SetTextColor(kBlack); txt_Ana2018_rms->SetTextColor(kBlack); txt_Ana2018_int->SetTextColor(kBlack);
263  txt_Ana2018->Draw(); txt_Ana2018_mean->Draw(); txt_Ana2018_rms->Draw(); txt_Ana2018_int->Draw();
264  vec_Means.push_back(f2_h0->GetMean ()); vec_Means.push_back(f2_h0->GetMeanError());
265  vec_Means.push_back(f2_h0->GetRMS ()); vec_Means.push_back(f2_h0->GetRMSError ());
266  vec_Means.push_back(f2_h0->Integral());
267 
268  vec_Means.push_back((f1_h1->GetMean()-f2_h0->GetMean())/f2_h0->GetMean()*100.);
269  vec_Means.push_back(getFractionalError (f2_h0->GetMean(), f2_h0->GetMeanError(), f1_h1->GetMean(), f1_h1->GetMeanError()));
270  vec_Means.push_back(is2SigmaBandsTouching(f2_h0->GetMean(), f2_h0->GetMeanError(), f1_h1->GetMean(), f1_h1->GetMeanError()));
271  vec_Means.push_back((f1_h1->GetMean()-f1_h0->GetMean())/f1_h0->GetMean()*100.);
272  vec_Means.push_back(getFractionalError (f1_h0->GetMean(), f1_h0->GetMeanError(), f1_h1->GetMean(), f1_h1->GetMeanError()));
273  vec_Means.push_back(is2SigmaBandsTouching(f1_h0->GetMean(), f1_h0->GetMeanError(), f1_h1->GetMean(), f1_h1->GetMeanError()));
274  vec_Means.push_back((f1_h0->GetMean()-f2_h0->GetMean())/f2_h0->GetMean()*100.);
275  vec_Means.push_back(getFractionalError (f2_h0->GetMean(), f2_h0->GetMeanError(), f1_h0->GetMean(), f1_h0->GetMeanError()));
276  vec_Means.push_back(is2SigmaBandsTouching(f2_h0->GetMean(), f2_h0->GetMeanError(), f1_h0->GetMean(), f1_h0->GetMeanError()));
277 
278  sprintf (s_FileName, "Var%d_Q%d_%s.pdf", S, Q+1, s_FileAppend.c_str());
279  c->SaveAs((s_OutDir+(std::string)s_FileName).c_str());
280 
281  delete c;
282 
283  //IF THE PLOT CONTAINS ALL QUANTILES AND QUARTILES, ADD IT TO KEY PLOT LIST.
284  if((Q+1 == 5) && (f1_spec1->GetLabels()[0].find("all slices")!=std::string::npos))
285  {
286  addToKeyPlotList(s_OutDir, s_OutDir+(std::string)s_FileName);
287  }
288 
289  //GENERATE A CSV FILE OF MEAN / RMS ETC DATA.
290  makeCSVFile(s_OutDir, s_FileAppend, s_FileName, f1_spec1->GetLabels()[0], Q+1, vec_Means);
291 
292  }
293  }
294 
295  std::cout << "\n\n\nFINISHED MAKING PLOTS!\n\n\n";
296 
297 }
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
#define S(x, n)
void makeCSVFile(std::string const &s_OutDir, std::string const &s_FileAppend, std::string const &s_FileName, std::string const &s_Title, unsigned int const &quantile, std::vector< double > const &vec_Means)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void Plotting_Data2DataComp_SingleCanvas(std::string s_OutDir="", std::string s_FileAppend="", std::string s_InFile1="", std::string s_InFile2="")
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:563
double is2SigmaBandsTouching(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
Represent the ratio between two spectra.
Definition: Ratio.h:8
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
OStream cout
Definition: OStream.cxx:6
const int nQuant
Definition: varsandcuts.h:4
enum BeamMode kViolet
double getFractionalError(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
app
Definition: demo.py:189
enum BeamMode kBlue
void addToKeyPlotList(std::string const &s_OutDir, std::string const &s_FileName)
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
enum BeamMode string