[u/mrichter/AliRoot.git] / PWGJE / macros / PID / AddUpSystematicErrors.C

#include "TH1D.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TString.h"
#include "TLegend.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TGraph.h"
#include "TMath.h"

#include "AliPID.h"

#include "THnSparseDefinitions.h"

#include <iostream>
#include <iomanip>

const Int_t numSpecies = 5;

//________________________________________________________
TCanvas* DrawFractionHistos(TString canvName, TString canvTitle, Double_t pLow, Double_t pHigh, TGraphAsymmErrors** gr,  TH1F** histRef)
{
  TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
  canv->SetGridx(1);
  canv->SetGridy(1);
  canv->SetLogx(1);
  
  for (Int_t i = 0; i < numSpecies; i++) {
    histRef[i]->GetYaxis()->SetRangeUser(0.0, 1.0);
    histRef[i]->GetYaxis()->SetTitle(canvTitle.Data());
    histRef[i]->GetXaxis()->SetRangeUser(pLow, pHigh);
    //histRef[i]->SetFillStyle(3004 + i);
    //histRef[i]->SetFillColor(kGray);
    histRef[i]->SetFillStyle(0);
    histRef[i]->SetFillColor(histRef[i]->GetMarkerColor());
    histRef[i]->SetLineColor(histRef[i]->GetMarkerColor());
  }
  histRef[2]->SetMarkerStyle(20);
  histRef[2]->Draw("e p");
  histRef[0]->SetMarkerStyle(21);
  histRef[0]->Draw("e p same");
  histRef[1]->SetMarkerStyle(22);
  histRef[1]->Draw("e p same");
  histRef[3]->SetMarkerStyle(29);
  histRef[3]->Draw("e p same");
  histRef[4]->SetMarkerStyle(30);
  histRef[4]->Draw("e p same");
  
  gr[0]->GetHistogram()->GetXaxis()->SetRangeUser(pLow, pHigh);
  gr[0]->GetHistogram()->GetYaxis()->SetRangeUser(0., 1.0);
  gr[0]->Draw("2 same");
  gr[1]->Draw("2 same");
  gr[2]->Draw("2 same");
  gr[3]->Draw("2 same");
  gr[4]->Draw("2 same");
  
  TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);    
  legend->SetBorderSize(0);
  legend->SetFillColor(0);
  legend->AddEntry(histRef[2], "#pi", "flp");
  legend->AddEntry(histRef[0], "e", "flp");
  legend->AddEntry(histRef[1], "K", "flp");
  legend->AddEntry(histRef[3], "p", "flp");
  legend->AddEntry(histRef[4], "#mu", "flp");
  legend->Draw();
  
  ClearTitleFromHistoInCanvas(canv);
  
  return canv;
}


//________________________________________________________
TGraphAsymmErrors* loadGraph(const TString graphName, TFile* f)
{
  if (!f) {
    std::cout << "No file. Cannot load graph \"" << graphName.Data() << "\n!" << std::endl;
    return 0x0;
  }
  
  TGraphAsymmErrors* grTemp = dynamic_cast<TGraphAsymmErrors*>(f->Get(graphName.Data()));
  if (!grTemp) {
    std::cout << "Failed to load histo \"" << graphName.Data() << "\"!" << std::endl;
    return 0x0;
  } 
  
  return grTemp;
}


//________________________________________________________
TH1F* loadHisto(const TString histName, TFile* f)
{
  if (!f) {
    std::cout << "No file. Cannot load hist \"" << histName.Data() << "\n!" << std::endl;
    return 0x0;
  }
  
  TH1F* hTemp = dynamic_cast<TH1F*>(f->Get(histName.Data()));
  if (!hTemp) {
    std::cout << "Failed to load histo \"" << histName.Data() << "\"!" << std::endl;
    return 0x0;
  } 
  
  return hTemp;
}


//________________________________________________________
Int_t AddUpSystematicErrors(const TString path, const TString outFileTitle, const TString* fileNames, const Int_t numFiles,
                            const TString fileNameReference) 
{ 
  if (!fileNames || numFiles < 1)
    return -1;
  
  TFile* f[numFiles];
  TGraphAsymmErrors** grSysError[numSpecies];
  TString grNames[numSpecies] = {"systematicError_electron", "systematicError_kaon", "systematicError_pion", "systematicError_proton",
                                 "systematicError_muon" };
  
  const TString histNames[numSpecies] = {"hFractionElectrons", "hFractionKaons", "hFractionPions", "hFractionProtons", "hFractionMuons" };
  
  
  TGraphAsymmErrors** grSysErrorYields[numSpecies];
  TString grNamesYields[numSpecies] = {"systematicErrorYields_electron", "systematicErrorYields_kaon", "systematicErrorYields_pion",
                                       "systematicErrorYields_proton", "systematicErrorYields_muon" };
                                 
  const TString histNamesYields[numSpecies] = {"hYieldElectrons", "hYieldKaons", "hYieldPions", "hYieldProtons", "hYieldMuons" };
  
  for (Int_t i = 0; i < numSpecies; i++) {
    grSysError[i] = new TGraphAsymmErrors*[numFiles];
    grSysErrorYields[i] = new TGraphAsymmErrors*[numFiles];
  }
    
  for (Int_t iFile = 0; iFile < numFiles; iFile++) {
    f[iFile] = TFile::Open(fileNames[iFile].Data());
    if (!f[iFile])  {
      std::cout << "Failed to open file \"" << fileNames[iFile].Data() << "\"!" << std::endl;
      return -1;
    }
        
    // Extract the data graphs
    for (Int_t i = 0; i < numSpecies; i++) {
      grSysError[i][iFile] = loadGraph(grNames[i], f[iFile]);
      if (!grSysError[i][iFile])
        return -1;
      
      grSysError[i][iFile]->SetName(Form("%s_fileID%d", grSysError[i][iFile]->GetName(), iFile));
      
      
      grSysErrorYields[i][iFile] = loadGraph(grNamesYields[i], f[iFile]);
      if (!grSysErrorYields[i][iFile])
        return -1;
      
      grSysErrorYields[i][iFile]->SetName(Form("%s_fileID%d", grSysErrorYields[i][iFile]->GetName(), iFile));
    }
  }
  
  // Load data points with statistical errors
  TFile* fReferenceData = TFile::Open(fileNameReference.Data());
  if (!fReferenceData) {
    std::cout << "Failed to open file \"" << fileNameReference.Data() << "\"!" << std::endl;
    return -1;
  }
    
  TH1F* hReferenceFractions[numSpecies];
  TH1F* hReferenceYields[numSpecies];
  
  for (Int_t i = 0; i < numSpecies; i++) {
    hReferenceFractions[i] = loadHisto(histNames[i], fReferenceData);
    if (!hReferenceFractions[i])
      return -1;
    
    hReferenceYields[i] = loadHisto(histNamesYields[i], fReferenceData);
    if (!hReferenceYields[i])
      return -1;
  }
  
  const Int_t reference = 0;
  
  // The x,y coordinates should be those of the reference graph. Then, all corresponding sys. errors are added in quadrature.
  TGraphAsymmErrors* grTotSysError[numSpecies];
  TGraphAsymmErrors* grTotSysErrorYields[numSpecies];
  Double_t sysErrorTotalSquared = 0;
  Double_t temp = 0;
  
  for (Int_t i = 0; i < numSpecies; i++) {
    grTotSysError[i] = new TGraphAsymmErrors(*grSysError[i][reference]);
    grTotSysError[i]->SetName(grNames[i]);
    
    for (Int_t iPoint = 0; iPoint < grTotSysError[i]->GetN(); iPoint++) {
      sysErrorTotalSquared = 0;
      for (Int_t iFile = 0; iFile < numFiles; iFile++) {
        // Already averages high and low value -> Since they are by now the same, this is ok.
        temp = grSysError[i][iFile]->GetErrorY(iPoint); 
        if (temp > 0) {
          sysErrorTotalSquared += temp * temp;
        }
      }
      
      grTotSysError[i]->SetPointEYlow(iPoint, TMath::Sqrt(sysErrorTotalSquared));
      grTotSysError[i]->SetPointEYhigh(iPoint, TMath::Sqrt(sysErrorTotalSquared));
    }
    
    // Same for the yields
    grTotSysErrorYields[i] = new TGraphAsymmErrors(*grSysErrorYields[i][reference]);
    grTotSysErrorYields[i]->SetName(grNamesYields[i]);
    
    for (Int_t iPoint = 0; iPoint < grTotSysErrorYields[i]->GetN(); iPoint++) {
      sysErrorTotalSquared = 0;
      for (Int_t iFile = 0; iFile < numFiles; iFile++) {
        // Already averages high and low value -> Since they are by now the same, this is ok.
        temp = grSysErrorYields[i][iFile]->GetErrorY(iPoint); 
        if (temp > 0) {
          sysErrorTotalSquared += temp * temp;
        }
      }
      
      grTotSysErrorYields[i]->SetPointEYlow(iPoint, TMath::Sqrt(sysErrorTotalSquared));
      grTotSysErrorYields[i]->SetPointEYhigh(iPoint, TMath::Sqrt(sysErrorTotalSquared));
    }
  }
  
  const Double_t pLow = 0.15;
  const Double_t pHigh = 50.;
  TCanvas* cFractionsWithTotalSystematicError = DrawFractionHistos("cFractionsWithTotalSystematicError", "Particle fractions", pLow, pHigh,
                                                                   grTotSysError, hReferenceFractions);
  
    
  // Output file
  TFile* fSave = 0x0;
  TDatime daTime;
  TString saveFileName;
    
  saveFileName = Form("outputSystematicsTotal_%s__%04d_%02d_%02d.root", outFileTitle.Data(), daTime.GetYear(),
                      daTime.GetMonth(), daTime.GetDay());
    
  fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");
  if (!fSave) {
    std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;
    return -1;
  }
  
  // Save final results
  fSave->cd();
  
  if (cFractionsWithTotalSystematicError)
      cFractionsWithTotalSystematicError->Write();
    
    for (Int_t i = 0; i < numSpecies; i++) {
      if (grTotSysError[i])
        grTotSysError[i]->Write();
      if (hReferenceFractions[i])
        hReferenceFractions[i]->Write();
      
      if (grTotSysErrorYields[i])
        grTotSysErrorYields[i]->Write();
      if (hReferenceYields[i])
        hReferenceYields[i]->Write();
  }
  
  // Save list of file names in output file
  TString listOfFileNames = "";
  for (Int_t i = 0; i < numFiles; i++) {
    listOfFileNames.Append(Form("%s%d: %s", i == 0 ? "" : ", ", i, fileNames[i].Data()));
  }
  
  TNamed* settings = new TNamed(Form("Used files for systematics: %s\n", listOfFileNames.Data()),
                                Form("Used files for systematics: %s\n", listOfFileNames.Data()));
  settings->Write();
  
  delete cFractionsWithTotalSystematicError;
    
  fSave->Close();
     
  return 0;
}
Commit	Line	Data
896d3200	1	#include "TH1D.h"
	2	#include "TCanvas.h"
	3	#include "TStyle.h"
	4	#include "TString.h"
	5	#include "TLegend.h"
	6	#include "TFile.h"
	7	#include "TGraphErrors.h"
	8	#include "TGraphAsymmErrors.h"
	9	#include "TGraph.h"
	10	#include "TMath.h"
	11
	12	#include "AliPID.h"
	13
	14	#include "THnSparseDefinitions.h"
	15
	16	#include <iostream>
	17	#include <iomanip>
	18
	19	const Int_t numSpecies = 5;
	20
	21	//________________________________________________________
	22	TCanvas* DrawFractionHistos(TString canvName, TString canvTitle, Double_t pLow, Double_t pHigh, TGraphAsymmErrors gr, TH1F histRef)
	23	{
	24	TCanvas* canv = new TCanvas(canvName.Data(), canvTitle.Data(),100,10,1200,800);
	25	canv->SetGridx(1);
	26	canv->SetGridy(1);
	27	canv->SetLogx(1);
	28
	29	for (Int_t i = 0; i < numSpecies; i++) {
	30	histRef[i]->GetYaxis()->SetRangeUser(0.0, 1.0);
	31	histRef[i]->GetYaxis()->SetTitle(canvTitle.Data());
	32	histRef[i]->GetXaxis()->SetRangeUser(pLow, pHigh);
	33	//histRef[i]->SetFillStyle(3004 + i);
	34	//histRef[i]->SetFillColor(kGray);
	35	histRef[i]->SetFillStyle(0);
	36	histRef[i]->SetFillColor(histRef[i]->GetMarkerColor());
	37	histRef[i]->SetLineColor(histRef[i]->GetMarkerColor());
	38	}
	39	histRef[2]->SetMarkerStyle(20);
	40	histRef[2]->Draw("e p");
	41	histRef[0]->SetMarkerStyle(21);
	42	histRef[0]->Draw("e p same");
	43	histRef[1]->SetMarkerStyle(22);
	44	histRef[1]->Draw("e p same");
	45	histRef[3]->SetMarkerStyle(29);
	46	histRef[3]->Draw("e p same");
	47	histRef[4]->SetMarkerStyle(30);
	48	histRef[4]->Draw("e p same");
	49
	50	gr[0]->GetHistogram()->GetXaxis()->SetRangeUser(pLow, pHigh);
	51	gr[0]->GetHistogram()->GetYaxis()->SetRangeUser(0., 1.0);
	52	gr[0]->Draw("2 same");
	53	gr[1]->Draw("2 same");
	54	gr[2]->Draw("2 same");
	55	gr[3]->Draw("2 same");
	56	gr[4]->Draw("2 same");
	57
	58	TLegend* legend = new TLegend(0.622126, 0.605932, 0.862069, 0.855932);
	59	legend->SetBorderSize(0);
	60	legend->SetFillColor(0);
	61	legend->AddEntry(histRef[2], "#pi", "flp");
	62	legend->AddEntry(histRef[0], "e", "flp");
	63	legend->AddEntry(histRef[1], "K", "flp");
	64	legend->AddEntry(histRef[3], "p", "flp");
65	legend->AddEntry(histRef[4], "#mu", "flp");
66	legend->Draw();
67
68	ClearTitleFromHistoInCanvas(canv);
69
70	return canv;
71	}
72
73
74	//________________________________________________________
75	TGraphAsymmErrors* loadGraph(const TString graphName, TFile* f)
76	{
77	if (!f) {
78	std::cout << "No file. Cannot load graph \"" << graphName.Data() << "\n!" << std::endl;
79	return 0x0;
80	}
81
82	TGraphAsymmErrors* grTemp = dynamic_cast<TGraphAsymmErrors*>(f->Get(graphName.Data()));
83	if (!grTemp) {
84	std::cout << "Failed to load histo \"" << graphName.Data() << "\"!" << std::endl;
85	return 0x0;
86	}
87
88	return grTemp;
89	}
90
91
92	//________________________________________________________
93	TH1F* loadHisto(const TString histName, TFile* f)
94	{
95	if (!f) {
96	std::cout << "No file. Cannot load hist \"" << histName.Data() << "\n!" << std::endl;
97	return 0x0;
98	}
99
100	TH1F* hTemp = dynamic_cast<TH1F*>(f->Get(histName.Data()));
101	if (!hTemp) {
102	std::cout << "Failed to load histo \"" << histName.Data() << "\"!" << std::endl;
103	return 0x0;
104	}
105
106	return hTemp;
107	}
108
109
110	//________________________________________________________
111	Int_t AddUpSystematicErrors(const TString path, const TString outFileTitle, const TString* fileNames, const Int_t numFiles,
112	const TString fileNameReference)
113	{
114	if (!fileNames \|\| numFiles < 1)
115	return -1;
116
117	TFile* f[numFiles];
118	TGraphAsymmErrors** grSysError[numSpecies];
119	TString grNames[numSpecies] = {"systematicError_electron", "systematicError_kaon", "systematicError_pion", "systematicError_proton",
120	"systematicError_muon" };
121
122	const TString histNames[numSpecies] = {"hFractionElectrons", "hFractionKaons", "hFractionPions", "hFractionProtons", "hFractionMuons" };
123
124
125	TGraphAsymmErrors** grSysErrorYields[numSpecies];
126	TString grNamesYields[numSpecies] = {"systematicErrorYields_electron", "systematicErrorYields_kaon", "systematicErrorYields_pion",
127	"systematicErrorYields_proton", "systematicErrorYields_muon" };
128
129	const TString histNamesYields[numSpecies] = {"hYieldElectrons", "hYieldKaons", "hYieldPions", "hYieldProtons", "hYieldMuons" };
130
131	for (Int_t i = 0; i < numSpecies; i++) {
132	grSysError[i] = new TGraphAsymmErrors*[numFiles];
133	grSysErrorYields[i] = new TGraphAsymmErrors*[numFiles];
134	}
135
136	for (Int_t iFile = 0; iFile < numFiles; iFile++) {
137	f[iFile] = TFile::Open(fileNames[iFile].Data());
138	if (!f[iFile]) {
139	std::cout << "Failed to open file \"" << fileNames[iFile].Data() << "\"!" << std::endl;
140	return -1;
141	}
142
143	// Extract the data graphs
144	for (Int_t i = 0; i < numSpecies; i++) {
145	grSysError[i][iFile] = loadGraph(grNames[i], f[iFile]);
146	if (!grSysError[i][iFile])
147	return -1;
148
149	grSysError[i][iFile]->SetName(Form("%s_fileID%d", grSysError[i][iFile]->GetName(), iFile));
150
151
152	grSysErrorYields[i][iFile] = loadGraph(grNamesYields[i], f[iFile]);
153	if (!grSysErrorYields[i][iFile])
154	return -1;
155
156	grSysErrorYields[i][iFile]->SetName(Form("%s_fileID%d", grSysErrorYields[i][iFile]->GetName(), iFile));
157	}
158	}
159
160	// Load data points with statistical errors
161	TFile* fReferenceData = TFile::Open(fileNameReference.Data());
162	if (!fReferenceData) {
163	std::cout << "Failed to open file \"" << fileNameReference.Data() << "\"!" << std::endl;
164	return -1;
165	}
166
167	TH1F* hReferenceFractions[numSpecies];
168	TH1F* hReferenceYields[numSpecies];
169
170	for (Int_t i = 0; i < numSpecies; i++) {
171	hReferenceFractions[i] = loadHisto(histNames[i], fReferenceData);
172	if (!hReferenceFractions[i])
173	return -1;
174
175	hReferenceYields[i] = loadHisto(histNamesYields[i], fReferenceData);
176	if (!hReferenceYields[i])
177	return -1;
178	}
179
180	const Int_t reference = 0;
181
182	// The x,y coordinates should be those of the reference graph. Then, all corresponding sys. errors are added in quadrature.
183	TGraphAsymmErrors* grTotSysError[numSpecies];
184	TGraphAsymmErrors* grTotSysErrorYields[numSpecies];
185	Double_t sysErrorTotalSquared = 0;
186	Double_t temp = 0;
187
188	for (Int_t i = 0; i < numSpecies; i++) {
189	grTotSysError[i] = new TGraphAsymmErrors(*grSysError[i][reference]);
190	grTotSysError[i]->SetName(grNames[i]);
191
192	for (Int_t iPoint = 0; iPoint < grTotSysError[i]->GetN(); iPoint++) {
193	sysErrorTotalSquared = 0;
194	for (Int_t iFile = 0; iFile < numFiles; iFile++) {
195	// Already averages high and low value -> Since they are by now the same, this is ok.
196	temp = grSysError[i][iFile]->GetErrorY(iPoint);
197	if (temp > 0) {
198	sysErrorTotalSquared += temp * temp;
199	}
200	}
201
202	grTotSysError[i]->SetPointEYlow(iPoint, TMath::Sqrt(sysErrorTotalSquared));
203	grTotSysError[i]->SetPointEYhigh(iPoint, TMath::Sqrt(sysErrorTotalSquared));
204	}
205
206	// Same for the yields
207	grTotSysErrorYields[i] = new TGraphAsymmErrors(*grSysErrorYields[i][reference]);
208	grTotSysErrorYields[i]->SetName(grNamesYields[i]);
209
210	for (Int_t iPoint = 0; iPoint < grTotSysErrorYields[i]->GetN(); iPoint++) {
211	sysErrorTotalSquared = 0;
212	for (Int_t iFile = 0; iFile < numFiles; iFile++) {
213	// Already averages high and low value -> Since they are by now the same, this is ok.
214	temp = grSysErrorYields[i][iFile]->GetErrorY(iPoint);
215	if (temp > 0) {
216	sysErrorTotalSquared += temp * temp;
217	}
218	}
219
220	grTotSysErrorYields[i]->SetPointEYlow(iPoint, TMath::Sqrt(sysErrorTotalSquared));
221	grTotSysErrorYields[i]->SetPointEYhigh(iPoint, TMath::Sqrt(sysErrorTotalSquared));
222	}
223	}
224
225	const Double_t pLow = 0.15;
226	const Double_t pHigh = 50.;
227	TCanvas* cFractionsWithTotalSystematicError = DrawFractionHistos("cFractionsWithTotalSystematicError", "Particle fractions", pLow, pHigh,
228	grTotSysError, hReferenceFractions);
229
230
231	// Output file
232	TFile* fSave = 0x0;
233	TDatime daTime;
234	TString saveFileName;
235
236	saveFileName = Form("outputSystematicsTotal_%s__%04d_%02d_%02d.root", outFileTitle.Data(), daTime.GetYear(),
237	daTime.GetMonth(), daTime.GetDay());
238
239	fSave = TFile::Open(Form("%s/%s", path.Data(), saveFileName.Data()), "recreate");
240	if (!fSave) {
241	std::cout << "Failed to open save file \"" << Form("%s/%s", path.Data(), saveFileName.Data()) << "\"!" << std::endl;
242	return -1;
243	}
244
245	// Save final results
246	fSave->cd();
247
248	if (cFractionsWithTotalSystematicError)
249	cFractionsWithTotalSystematicError->Write();
250
251	for (Int_t i = 0; i < numSpecies; i++) {
252	if (grTotSysError[i])
253	grTotSysError[i]->Write();
254	if (hReferenceFractions[i])
255	hReferenceFractions[i]->Write();
256
257	if (grTotSysErrorYields[i])
258	grTotSysErrorYields[i]->Write();
259	if (hReferenceYields[i])
260	hReferenceYields[i]->Write();
261	}
262
263	// Save list of file names in output file
264	TString listOfFileNames = "";
265	for (Int_t i = 0; i < numFiles; i++) {
266	listOfFileNames.Append(Form("%s%d: %s", i == 0 ? "" : ", ", i, fileNames[i].Data()));
267	}
268
269	TNamed* settings = new TNamed(Form("Used files for systematics: %s\n", listOfFileNames.Data()),
270	Form("Used files for systematics: %s\n", listOfFileNames.Data()));
271	settings->Write();
272
273	delete cFractionsWithTotalSystematicError;
274
275	fSave->Close();
276
277	return 0;
278	}