adding nsigma study
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / drawSystematics.C
1 /* $Id$ */
2
3 #if !defined(__CINT__) || defined(__MAKECINT__)
4
5 #include "AliPWG0Helper.h"
6 #include "dNdEtaAnalysis.h"
7 #include "AlidNdEtaCorrection.h"
8
9 #include <TCanvas.h>
10 #include <TFile.h>
11 #include <TH1.h>
12 #include <TH2F.h>
13 #include <TH3F.h>
14 #include <TLine.h>
15 #include <TSystem.h>
16
17 #endif
18
19 extern TPad* gPad;
20
21 void SetRanges(TAxis* axis)
22 {
23   if (strcmp(axis->GetTitle(), "#eta") == 0)
24     axis->SetRangeUser(-1.7999, 1.7999);
25   if (strcmp(axis->GetTitle(), "p_{T} [GeV/c]") == 0)
26     axis->SetRangeUser(0, 9.9999);
27   if (strcmp(axis->GetTitle(), "vtx z [cm]") == 0)
28     axis->SetRangeUser(-15, 14.9999);
29   if (strcmp(axis->GetTitle(), "Ntracks") == 0)
30     axis->SetRangeUser(0, 99.9999);
31 }
32
33 void SetRanges(TH1* hist)
34 {
35   SetRanges(hist->GetXaxis());
36   SetRanges(hist->GetYaxis());
37   SetRanges(hist->GetZaxis());
38 }
39
40 void Prepare3DPlot(TH3* hist)
41 {
42   hist->GetXaxis()->SetTitleOffset(1.5);
43   hist->GetYaxis()->SetTitleOffset(1.5);
44   hist->GetZaxis()->SetTitleOffset(1.5);
45
46   hist->SetStats(kFALSE);
47 }
48
49 void Prepare2DPlot(TH2* hist)
50 {
51   hist->SetStats(kFALSE);
52   hist->GetYaxis()->SetTitleOffset(1.4);
53
54   SetRanges(hist);
55 }
56
57 void Prepare1DPlot(TH1* hist)
58 {
59   hist->SetLineWidth(2);
60   hist->SetStats(kFALSE);
61
62   SetRanges(hist);
63 }
64
65 void InitPad()
66 {
67   if (!gPad)
68     return;
69
70   gPad->Range(0, 0, 1, 1);
71   gPad->SetLeftMargin(0.15);
72   //gPad->SetRightMargin(0.05);
73   //gPad->SetTopMargin(0.13);
74   //gPad->SetBottomMargin(0.1);
75
76   //gPad->SetGridx();
77   //gPad->SetGridy();
78 }
79
80 void InitPadCOLZ()
81 {
82   gPad->Range(0, 0, 1, 1);
83   gPad->SetRightMargin(0.15);
84   gPad->SetLeftMargin(0.12);
85
86   gPad->SetGridx();
87   gPad->SetGridy();
88 }
89
90 void Secondaries()
91 {
92   TFile* file = TFile::Open("systematics.root");
93
94   TH3F* secondaries = dynamic_cast<TH3F*> (file->Get("fSecondaries"));
95   if (!secondaries)
96   {
97     printf("Could not read histogram\n");
98     return;
99   }
100
101   for (Int_t ptBin=1; ptBin<=secondaries->GetNbinsZ(); ptBin++)
102   //for (Int_t ptBin = 1; ptBin<=2; ptBin++)
103   {
104     TH1F* hist = new TH1F(Form("secondaries_%d", ptBin), Form("secondaries_%d", ptBin), secondaries->GetNbinsY(), secondaries->GetYaxis()->GetXmin(), secondaries->GetYaxis()->GetXmax());
105
106     hist->SetTitle(Form("%f < p_{T} < %f", secondaries->GetZaxis()->GetBinLowEdge(ptBin), secondaries->GetZaxis()->GetBinUpEdge(ptBin)));
107
108     for (Int_t cBin=1; cBin<=secondaries->GetNbinsY(); ++cBin)
109     {
110       if (secondaries->GetBinContent(0, cBin, ptBin) > 0)
111         printf("WARNING: Underflow bin not empty!");
112       if (secondaries->GetBinContent(secondaries->GetNbinsX()+1, cBin, ptBin) > 0)
113         printf("WARNING: Overflow bin not empty!");
114
115       Double_t sum = 0;
116       Double_t count = 0;
117       for (Int_t nBin=1; nBin<=secondaries->GetNbinsX(); ++nBin)
118       {
119         //printf("%f %f\n", secondaries->GetXaxis()->GetBinCenter(nBin), secondaries->GetBinContent(nBin, cBin, ptBin));
120         sum += secondaries->GetXaxis()->GetBinCenter(nBin) * secondaries->GetBinContent(nBin, cBin, ptBin);
121         count += secondaries->GetBinContent(nBin, cBin, ptBin);
122       }
123
124       printf("%f %f\n", sum, count);
125
126       if (count > 0)
127       {
128         hist->SetBinContent(cBin, sum / count);
129         hist->SetBinError(cBin, TMath::Sqrt(sum) / count);
130       }
131     }
132
133     hist->Scale(1.0 / hist->GetBinContent(hist->GetXaxis()->FindBin(1)));
134     hist->Add(new TF1("flat", "-1", 0, 2));
135
136     new TCanvas;
137     hist->SetMarkerStyle(21);
138     hist->Draw("");
139   }
140 }
141
142 void Track2Particle1DComposition(const char* fileName = "correction_map.root", Int_t folderCount, const char** folderNames, Float_t upperPtLimit = 9.9)
143 {
144   gSystem->Load("libPWG0base");
145
146   TString canvasName;
147   canvasName.Form("Track2Particle1DComposition");
148   TCanvas* canvas = new TCanvas(canvasName, canvasName, 1200, 400);
149   canvas->Divide(3, 1);
150
151   TLegend* legend = new TLegend(0.7, 0.7, 0.95, 0.95);
152
153   for (Int_t i=0; i<folderCount; ++i)
154   {
155     Track2Particle1DCreatePlots(fileName, folderNames[i], upperPtLimit);
156
157     TH1* corrX = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_x_div_meas_nTrackToNPart_x"));
158     TH1* corrY = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_y_div_meas_nTrackToNPart_y"));
159     TH1* corrZ = dynamic_cast<TH1*> (gROOT->FindObject("gene_nTrackToNPart_z_div_meas_nTrackToNPart_z"));
160
161     Prepare1DPlot(corrX);
162     Prepare1DPlot(corrY);
163     Prepare1DPlot(corrZ);
164
165     const char* title = "Track2Particle Correction";
166     corrX->SetTitle(title);
167     corrY->SetTitle(title);
168     corrZ->SetTitle(title);
169
170     corrZ->GetXaxis()->SetRangeUser(0, upperPtLimit);
171
172     corrX->SetLineColor(i+1);
173     corrY->SetLineColor(i+1);
174     corrZ->SetLineColor(i+1);
175
176     canvas->cd(1);
177     InitPad();
178     corrX->DrawCopy(((i>0) ? "SAME" : ""));
179
180     canvas->cd(2);
181     InitPad();
182     corrY->DrawCopy(((i>0) ? "SAME" : ""));
183
184     canvas->cd(3);
185     InitPad();
186     corrZ->DrawCopy(((i>0) ? "SAME" : ""));
187
188     legend->AddEntry(corrZ, folderNames[i]);
189   }
190
191   legend->Draw();
192
193   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.gif", fileName, gMax, upperPtLimit));
194   //canvas->SaveAs(Form("Track2Particle1D_%s_%d_%f.eps", fileName, gMax, upperPtLimit));
195 }
196
197 void DrawDifferentSpecies()
198 {
199   gROOT->ProcessLine(".L drawPlots.C");
200
201   const char* folderNames[] = { "correction_0", "correction_1", "correction_2", "correction_3" };
202
203   Track2Particle1DComposition("systematics.root", 4, folderNames);
204 }
205
206 void ScalePtDependent(TH3F* hist, TF1* function)
207 {
208   // assumes that pt is the third dimension of hist
209   // scales with function(pt)
210
211   for (Int_t z=1; z<=hist->GetNbinsZ(); ++z)
212   {
213     Double_t factor = function->Eval(hist->GetZaxis()->GetBinCenter(z));
214     printf("z = %d, pt = %f, scaling with %f\n", z, hist->GetZaxis()->GetBinCenter(z), factor);
215
216     for (Int_t x=1; x<=hist->GetNbinsX(); ++x)
217       for (Int_t y=1; y<=hist->GetNbinsY(); ++y)
218         hist->SetBinContent(x, y, z, hist->GetBinContent(x, y, z) * factor);
219   }
220 }
221
222 void ChangeComposition(void** correctionPointer, Int_t index)
223 {
224   AlidNdEtaCorrection** fdNdEtaCorrection = (AlidNdEtaCorrection**) correctionPointer;
225
226   switch (index)
227   {
228     case 0:
229       break;
230
231     case 1:
232       fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Scale(10);
233       fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Scale(10);
234       break;
235
236     case 2:
237       TF1* ptDependence = new TF1("simple", "x", 0, 100);
238       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetGeneratedHistogram(), ptDependence);
239       ScalePtDependent(fdNdEtaCorrection[0]->GetTrack2ParticleCorrection()->GetMeasuredHistogram(), ptDependence);
240       break;
241
242   }
243 }
244
245 void Composition()
246 {
247   gSystem->Load("libPWG0base");
248
249   gSystem->Unlink("new_compositions.root");
250
251   for (Int_t comp = 0; comp < 3; ++comp)
252   {
253     AlidNdEtaCorrection* fdNdEtaCorrection[4];
254
255     TFile::Open("systematics.root");
256
257     for (Int_t i=0; i<4; ++i)
258     {
259       TString name;
260       name.Form("correction_%d", i);
261       fdNdEtaCorrection[i] = new AlidNdEtaCorrection(name, name);
262       fdNdEtaCorrection[i]->LoadHistograms("systematics.root", name);
263     }
264
265     ChangeComposition(fdNdEtaCorrection, comp);
266
267     Double_t geneCount[5];
268     Double_t measCount[5];
269     geneCount[4] = 0;
270     measCount[4] = 0;
271
272     for (Int_t i=0; i<4; ++i)
273     {
274       geneCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetGeneratedHistogram()->Integral();
275       measCount[i] = fdNdEtaCorrection[i]->GetTrack2ParticleCorrection()->GetMeasuredHistogram()->Integral();
276
277       geneCount[4] += geneCount[i];
278       measCount[4] += measCount[i];
279
280       printf("Particle %s: %d gene, %d meas\n", ((i == 0) ? "pi" : (i == 1) ? "K" : (i == 2) ? "p" : "others"), (Int_t) geneCount[i], (Int_t) measCount[i]);
281     }
282
283     printf("Generated ratios are:     %f pi, %f K, %f p, %f others\n", geneCount[0] / geneCount[4], geneCount[1] / geneCount[4], geneCount[2] / geneCount[4], geneCount[3] / geneCount[4]);
284
285     printf("Reconstructed ratios are: %f pi, %f K, %f p, %f others\n", measCount[0] / measCount[4], measCount[1] / measCount[4], measCount[2] / measCount[4], measCount[3] / measCount[4]);
286
287     TList* collection = new TList;
288
289     for (Int_t i=0; i<4; ++i)
290       collection->Add(fdNdEtaCorrection[i]);
291
292     TString correctionName;
293     correctionName.Form("new_composition_%d", comp);
294
295     AlidNdEtaCorrection* newComposition = new AlidNdEtaCorrection(correctionName, correctionName);
296     newComposition->Merge(collection);
297     newComposition->Finish();
298
299     delete collection;
300
301     TFile* file = TFile::Open("new_compositions.root", "UPDATE");
302     newComposition->SaveHistograms();
303     //file->Write();
304     file->Close();
305   }
306
307   gROOT->ProcessLine(".L drawPlots.C");
308
309   const char* folderNames[] = { "new_composition_0", "new_composition_1", "new_composition_2" };
310
311   Track2Particle1DComposition("new_compositions.root", 3, folderNames);
312 }
313
314 Double_t ConvSigma1To2D(Double_t sigma)
315 {
316   return TMath::Sqrt( - TMath::Log( 1 - TMath::Erf(sigma / TMath::Sqrt(2)) ) * 2);
317 }
318
319 Double_t ConvDistance1DTo2D(Double_t distance)
320 {
321   return TMath::ErfInverse(1 - TMath::Exp(-distance * distance / 2)) * TMath::Sqrt(2);
322 }
323
324 Double_t Sigma2VertexCount(TH2F* tracks, Double_t nSigma)
325 {
326   Double_t count = 0;
327
328   //nSigma = ConvSigma1To2D(nSigma);
329
330   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
331     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
332     {
333       Double_t impactX = tracks->GetXaxis()->GetBinCenter(x);
334       Double_t impactY = tracks->GetYaxis()->GetBinCenter(y);
335
336       Float_t d = TMath::Sqrt(impactX*impactX + impactY*impactY);
337
338       d = ConvDistance1DTo2D(d);
339
340       if (d < nSigma)
341         count += tracks->GetBinContent(x, y);
342     }
343
344   return count;
345 }
346
347 TH2F* Sigma2VertexGaussianTracksHist()
348 {
349   TH2F* tracks = new TH2F("Sigma2Vertex_tracks", "Sigma2Vertex_tracks", 200, -5, 5, 200, -5, 5);
350
351   TF2* gaussian2D = new TF2("gaussian2D", "xgausn(0) * ygausn(3)", -5, 5, -5, 5);
352   gaussian2D->SetParameters(1, 0, 1, 1, 0, 1);
353
354   for (Int_t x=1; x<=tracks->GetNbinsX(); ++x)
355     for (Int_t y=1; y<=tracks->GetNbinsY(); ++y)
356       tracks->SetBinContent(x, y, gaussian2D->Eval(tracks->GetXaxis()->GetBinCenter(x), tracks->GetYaxis()->GetBinCenter(y)));
357
358   //normalize
359   tracks->Scale(1.0 / tracks->Integral());
360
361   return tracks;
362 }
363
364 void Sigma2VertexGaussian()
365 {
366   TH2F* tracks = Sigma2VertexGaussianTracksHist();
367
368   TCanvas* canvas = new TCanvas("Sigma2VertexGaussian", "Sigma2VertexGaussian", 1000, 1000);
369   canvas->Divide(2, 2);
370
371   canvas->cd(1);
372   tracks->Draw("COLZ");
373
374   TH1F* ratio = new TH1F("Sigma2Vertex_ratio", "Sigma2Vertex_ratio;n sigma;included", 10, 0.25, 5.25);
375   for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
376     ratio->Fill(nSigma, Sigma2VertexCount(tracks, nSigma));
377   ratio->SetMarkerStyle(21);
378
379   canvas->cd(2);
380   ratio->Draw("P");
381
382   TH1F* ratio2 = new TH1F("Sigma2Vertex_ratio2", "Sigma2Vertex_ratio2;nSigma;% included 3 sigma / % included n sigma", 10, 0.25, 5.25);
383   Double_t sigma3 = Sigma2VertexCount(tracks, 3);
384   for (Double_t nSigma = 0.5; nSigma < 5.1; nSigma += 0.5)
385     ratio2->Fill(nSigma, sigma3 / Sigma2VertexCount(tracks, nSigma));
386   ratio2->SetMarkerStyle(21);
387
388   canvas->cd(3);
389   ratio2->Draw("P");
390
391   canvas->SaveAs("Sigma2Vertex.eps");
392 }
393
394 void Sigma2VertexSimulation()
395 {
396   TFile* file = TFile::Open("systematics.root");
397
398   TH1F* sigmavertex = dynamic_cast<TH1F*> (file->Get("fSigmaVertex"));
399   if (!sigmavertex)
400   {
401     printf("Could not read histogram\n");
402     return;
403   }
404
405   TH1F* ratio = new TH1F("sigmavertex_ratio", "sigmavertex_ratio;Nsigma;% included 3 sigma / % included n sigma", sigmavertex->GetNbinsX(), sigmavertex->GetXaxis()->GetXmin(), sigmavertex->GetXaxis()->GetXmax());
406
407   for (Int_t i=1; i<=sigmavertex->GetNbinsX(); ++i)
408     ratio->SetBinContent(i, sigmavertex->GetBinContent(sigmavertex->GetXaxis()->FindBin(3)) / sigmavertex->GetBinContent(i));
409
410   TCanvas* canvas = new TCanvas("Sigma2VertexSimulation", "Sigma2VertexSimulation", 1000, 500);
411   canvas->Divide(2, 1);
412
413   canvas->cd(1);
414   sigmavertex->SetMarkerStyle(21);
415   sigmavertex->Draw("P");
416
417   canvas->cd(2);
418   ratio->SetMarkerStyle(21);
419   ratio->Draw("P");
420 }
421
422
423 void drawSystematics()
424 {
425   //Secondaries();
426   //DrawDifferentSpecies();
427   //Composition();
428
429   Sigma2VertexSimulation();
430 }