Plotting_Data2DataComp.C
Go to the documentation of this file.
1 /*******************************************************************
2  * A SCRIPT TO READ THE OUTPUT ROOT FILE OF Analyse_Data2DataComp.C,
3  * WILL PLOT SPECTRA FROM BOTH DATA SETS AND THE RATIO OF THE TWO.
4  * OUTPUT CSV FILE CREATED WITH KEY QUANTITIES SUCH AS THE FRACTIONAL
5  * CHANGES IN THE MEANS ETC.
6  * USAGE:
7  * Plotting_Data2DataComp.C <Input root file containing spectra>
8  * <Directory to output pdfs>
9  * <String to append outfile names>
10  * <int color for first spectra>
11  * <int color for second spectra>
12  ********************************************************************/
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  std::cout << numerator << " " << denominator << " " << numeratorVar << " " << denominatorVar << std::endl;
50 
51  return std::sqrt((numerator/denominator)*(numerator/denominator)*(numeratorVar/(numerator*numerator)+denominatorVar/(denominator*denominator)))*100.;
52 }
53 
54 std::string is2SigmaBandsTouching(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
55 {
56  std::pair<double,double> largeMean;
57  std::pair<double,double> smallMean;
58 
59  if(mean0 > mean1)
60  {
61  largeMean = {mean0, meanErr0};
62  smallMean = {mean1, meanErr1};
63  }
64  else
65  {
66  largeMean = {mean1, meanErr1};
67  smallMean = {mean0, meanErr0};
68  }
69 
70  double largeMeanShifted = largeMean.first - 2.*largeMean.second;
71  double smallMeanShifted = smallMean.first + 2.*smallMean.second;
72 
73  return (largeMeanShifted <= smallMeanShifted ? "True" : "False");
74 }
75 
76 void Plotting_Data2DataComp(std::string s_InFile = "", std::string s_OutDir = "", std::string s_OutName = "",
77  unsigned int color1 = 1, unsigned int color2 = 2) {
78 
79  std::cout << "\n\n\nMAKING THE PLOTS...\n\n\n";
80 
81  std::ofstream f_MeansFile;
82  f_MeansFile.open(s_OutDir+"meansFile_"+s_OutName+".csv");
83 
84  TFile f_In(s_InFile.c_str());
85  char s_DirName [256];
86  char s_FileName[256];
87  unsigned int nQuant = 4;
88  unsigned int nSpec = 41;
89 
90  for(unsigned int S = 0; S <= nSpec; ++S) {
91  for(unsigned int Q = 0; Q <= nQuant; ++Q) {
92 
93  sprintf(s_DirName, "dir_Var%d_Q%d_DataSet0",S,Q);
94  std::unique_ptr<Spectrum> spec0 = Spectrum::LoadFrom(f_In, s_DirName);
95  sprintf(s_DirName, "dir_Var%d_Q%d_DataSet1",S,Q);
96  std::unique_ptr<Spectrum> spec1 = Spectrum::LoadFrom(f_In, s_DirName);
97 
98  TCanvas *c = new TCanvas("c", "c", 1000, 800);
99  TPad *p1 = new TPad("p1", "p1", 0, 0, 1, 1);
100  TPad *p2 = new TPad("p2", "p2", 0, 0, 1, 1);
101  SplitCanvas(0.3,p1,p2);
102 
103  p1->cd();
104  //p1->SetLogy();
105  gStyle->SetOptStat(0);
106 
107  TH1 *h0 = spec0->ToTH1(spec1->POT());
108  TH1 *h1 = spec1->ToTH1(spec1->POT());
109 
110  h0->SetLineWidth(2);
111  h0->SetLineColor(color1);
112 
113  h1->SetLineWidth(2);
114  h1->SetLineColor(color2);
115 
116  h0->GetXaxis()->SetTitle("");
117  h0->GetXaxis()->SetLabelSize(0);
118  h1->GetXaxis()->SetTitle("");
119  h1->GetXaxis()->SetLabelSize(0);
120 
121  h0->Draw();
122  h1->Draw("SAME");
123 
124  p2->cd();
125  Ratio rat(*spec1, *spec0);
126  std::cout << "DIVINING " << spec1->GetLabels()[0] << " BY " << spec0->GetLabels()[0] << std::endl;
127 
128  TH1* hRat = rat.ToTH1();
129  hRat->SetAxisRange(0.9,1.1,"Y");
130  hRat->Draw();
131 
132  double x[2] = {hRat->GetXaxis()->GetXmin(), hRat->GetXaxis()->GetXmax()};
133  double y[2] = {1.0, 1.0};
134  TGraph* l_Xaxis = new TGraph(2, x, y);
135  l_Xaxis->SetLineColor(kGray);
136  l_Xaxis->SetLineWidth(3);
137  l_Xaxis->SetLineStyle(2);
138  l_Xaxis->Draw("L SAME");
139 
140  sprintf (s_FileName, "Var%d_Q%d_%s.pdf", S, Q+1, s_OutName.c_str());
141  c->Print((s_OutDir+(std::string)s_FileName).c_str());
142 
143  std::cout << "\n\nMeans for:\t" << hRat->GetXaxis()->GetTitle() << " for quant " << Q+1 << "\n"
144  << "old data: \t" << h0->GetMean() << " +/- " << h0->GetMeanError() << "\n"
145  << "keepup data:\t" << h1->GetMean() << " +/- " << h1->GetMeanError() << "\n\n";
146 
147  f_MeansFile << s_FileName << "=" << hRat->GetXaxis()->GetTitle() << " for quant " << Q+1 << "="
148  << h0->GetMean() << "=" << h0->GetMeanError() << "="
149  << h1->GetMean() << "=" << h1->GetMeanError() << "="
150  << (h1->GetMean()-h0->GetMean())/h0->GetMean()*100. << "="
151  << getFractionalError (h0->GetMean(), h0->GetMeanError(), h1->GetMean(), h1->GetMeanError()) << "="
152  << is2SigmaBandsTouching(h0->GetMean(), h0->GetMeanError(), h1->GetMean(), h1->GetMeanError()) << std::endl;
153 
154  delete l_Xaxis;
155  delete c;
156  }
157  }
158 
159  f_MeansFile.close();
160 
161  std::cout << "\n\n\nFINISHED MAKING PLOTS!\n\n\n";
162 
163 }
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)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SplitCanvas(double ysplit, TPad *&p1, TPad *&p2)
Split the current pad into two vertically stacked pieces, suitable for ratio plots and so on...
Definition: Plots.cxx:1417
void Plotting_Data2DataComp(std::string s_InFile="", std::string s_OutDir="", std::string s_OutName="", unsigned int color1=1, unsigned int color2=2)
double getFractionalError(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
TH1F * h1
const int nQuant
Definition: varsandcuts.h:4
std::string is2SigmaBandsTouching(double const mean0, double const meanErr0, double const mean1, double const meanErr1)
enum BeamMode string