updating makefile to include evgen headers to compile systematics selector
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
1 /* $Id$ */
2
3 #include "dNdEtaAnalysis.h"
4
5 #include <TFile.h>
6 #include <TH3F.h>
7 #include <TH2D.h>
8 #include <TH1D.h>
9 #include <TMath.h>
10 #include <TCanvas.h>
11 #include <TCollection.h>
12 #include <TIterator.h>
13 #include <TList.h>
14 #include <TLegend.h>
15
16 #include "AlidNdEtaCorrection.h"
17 #include "AliPWG0Helper.h"
18
19 //____________________________________________________________________
20 ClassImp(dNdEtaAnalysis)
21
22 //____________________________________________________________________
23 dNdEtaAnalysis::dNdEtaAnalysis() :
24   TNamed(),
25   fData(0),
26   fDataUncorrected(0),
27   fVtx(0),
28   fPtDist(0)
29 {
30   // default constructor
31
32   for (Int_t i=0; i<kVertexBinning; ++i)
33   {
34     fdNdEta[i] = 0;
35     fdNdEtaPtCutOffCorrected[i] = 0;
36   }
37 }
38
39 //____________________________________________________________________
40 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
41   TNamed(name, title),
42   fData(0),
43   fDataUncorrected(0),
44   fVtx(0),
45   fPtDist(0)
46 {
47   // constructor
48
49   fData  = new TH3F(Form("%s_analysis", name),"Input data",80,-20,20,40,-2,2,100, 0, 10);
50   fData->SetXTitle("vtx z [cm]");
51   fData->SetYTitle("#eta");
52   fData->SetZTitle("p_{T}");
53   fData->Sumw2();
54
55   fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
56   fDataUncorrected->SetTitle(Form("%s uncorrected", fData->GetTitle()));
57   fVtx       = dynamic_cast<TH1D*> (fData->Project3D("x"));
58   fVtx->SetName(Form("%s_vtx", name));
59   fVtx->SetTitle("Vertex distribution");
60   fVtx->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
61   //fVtx->Sumw2();
62
63   fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
64   fdNdEta[0]->SetName(Form("%s_dNdEta", name));
65   fdNdEta[0]->SetTitle("dN_{ch}/d#eta");
66   fdNdEta[0]->GetXaxis()->SetTitle(fData->GetYaxis()->GetTitle());
67   fdNdEta[0]->SetYTitle("dN_{ch}/d#eta");
68
69   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
70
71   for (Int_t i=1; i<kVertexBinning; ++i)
72   {
73     fdNdEta[i]    = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
74     fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
75   }
76
77   fPtDist = dynamic_cast<TH1D*> (fData->Project3D("z"));
78   fPtDist->SetName(Form("%s_pt", name));
79   fPtDist->SetTitle("p_{T} [GeV/c]");
80   fPtDist->GetXaxis()->SetTitle(fData->GetZaxis()->GetTitle());
81   fPtDist->SetYTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
82 }
83
84 //____________________________________________________________________
85 dNdEtaAnalysis::~dNdEtaAnalysis()
86 {
87   // destructor
88
89   if (fData)
90   {
91     delete fData;
92     fData = 0;
93   }
94
95   if (fDataUncorrected)
96   {
97     delete fDataUncorrected;
98     fDataUncorrected = 0;
99   }
100
101   if (fVtx)
102   {
103     delete fVtx;
104     fVtx = 0;
105   }
106
107   for (Int_t i=0; i<kVertexBinning; ++i)
108   {
109     if (fdNdEta[i])
110     {
111       delete fdNdEta[i];
112       fdNdEta[i] = 0;
113     }
114     if (fdNdEtaPtCutOffCorrected[i])
115     {
116       delete fdNdEtaPtCutOffCorrected[i];
117       fdNdEtaPtCutOffCorrected[i] = 0;
118     }
119   }
120
121   if (fPtDist)
122   {
123     delete fPtDist;
124     fPtDist = 0;
125   }
126 }
127
128 //_____________________________________________________________________________
129 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
130   TNamed(c),
131   fData(0),
132   fDataUncorrected(0),
133   fVtx(0),
134   fPtDist(0)
135 {
136   //
137   // dNdEtaAnalysis copy constructor
138   //
139
140   ((dNdEtaAnalysis &) c).Copy(*this);
141 }
142
143 //_____________________________________________________________________________
144 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
145 {
146   //
147   // Assignment operator
148   //
149
150   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
151   return *this;
152 }
153
154 //_____________________________________________________________________________
155 void dNdEtaAnalysis::Copy(TObject &c) const
156 {
157   //
158   // Copy function
159   //
160
161   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
162
163   target.fData = dynamic_cast<TH3F*> (fData->Clone());
164   target.fDataUncorrected = dynamic_cast<TH3F*> (fDataUncorrected->Clone());
165   target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
166
167   for (Int_t i=0; i<kVertexBinning; ++i)
168   {
169     target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
170     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[i]->Clone());
171   }
172
173   target.fPtDist = dynamic_cast<TH1D*> (fPtDist->Clone());
174
175   TNamed::Copy((TNamed &) c);
176 }
177
178 //____________________________________________________________________
179 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight)
180 {
181   // fills a track into the histograms
182
183   fDataUncorrected->Fill(vtx, eta, pt);
184   fData->Fill(vtx, eta, pt, weight);
185 }
186
187 //____________________________________________________________________
188 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
189 {
190   // fills an event into the histograms
191
192   fVtx->Fill(vtx, weight);
193 }
194
195 //____________________________________________________________________
196 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
197 {
198   // correct with correction values if available
199
200   // TODO what do we do with the error?
201   if (!correction)
202     printf("INFO: No correction applied\n");
203
204   // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
205
206   // create pt hist
207   {
208     // reset all ranges
209     fData->GetXaxis()->SetRange(0, 0);
210     fData->GetYaxis()->SetRange(0, 0);
211     fData->GetZaxis()->SetRange(0, 0);
212
213     // vtx cut
214     Int_t vertexBinBegin = fData->GetXaxis()->FindBin(-5);
215     Int_t vertexBinEnd = fData->GetXaxis()->FindBin(5);
216     fData->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
217     Float_t nEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd);
218
219     // eta cut
220     fData->GetYaxis()->SetRange(fData->GetYaxis()->FindBin(-0.8), fData->GetYaxis()->FindBin(0.8));
221     Float_t etaWidth = 1.6;
222
223     TH1D* ptHist = dynamic_cast<TH1D*> (fData->Project3D("ze"));
224
225     for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
226     {
227       Float_t binSize = fPtDist->GetBinWidth(i);
228       fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
229       fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
230     }
231
232     delete ptHist;
233   }
234
235   // reset all ranges
236   fData->GetXaxis()->SetRange(0, 0);
237   fData->GetYaxis()->SetRange(0, 0);
238   fData->GetZaxis()->SetRange(0, 0);
239
240   // integrate over pt (with pt cut)
241   Int_t ptLowBin = 1;
242   if (ptCut > 0)
243     ptLowBin = fData->GetZaxis()->FindBin(ptCut);
244
245   fData->GetZaxis()->SetRange(ptLowBin, fData->GetZaxis()->GetNbins());
246   TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2e"));
247   vtxVsEta->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
248   vtxVsEta->GetYaxis()->SetTitle(fData->GetYaxis()->GetTitle());
249
250   if (vtxVsEta == 0)
251   {
252     printf("ERROR: pt integration failed\n");
253     return;
254   }
255
256   //new TCanvas;
257   //vtxVsEta->Draw("COLZ");
258
259   for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
260   {
261     // do we have several histograms for different vertex positions?
262     Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
263     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
264     {
265       Int_t vertexBinBegin = 1;
266       Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
267
268       // the first histogram is always for the whole vertex range
269       if (vertexPos > 0)
270       {
271         vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
272         vertexBinEnd = vertexBinBegin + vertexBinWidth;
273       }
274
275       Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
276       if (totalEvents == 0)
277       {
278         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
279         continue;
280       }
281
282       Float_t sum = 0;
283       Float_t sumError2 = 0;
284       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
285       {
286         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
287         {
288           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
289           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
290         }
291       }
292
293       Float_t ptCutOffCorrection = 1;
294       if (correction)
295         ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
296
297       if (ptCutOffCorrection <= 0)
298       {
299         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
300         continue;
301       }
302
303       Float_t dndeta = sum / totalEvents;
304       Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
305
306       dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
307       error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
308
309       fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
310       fdNdEta[vertexPos]->SetBinError(iEta, error);
311
312       dndeta /= ptCutOffCorrection;
313       error  /= ptCutOffCorrection;
314
315       fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(iEta, dndeta);
316       fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(iEta, error);
317     }
318   }
319 }
320
321 //____________________________________________________________________
322 void dNdEtaAnalysis::SaveHistograms()
323 {
324   // save the histograms to a directory with the name of this class (retrieved from TNamed)
325
326   gDirectory->mkdir(GetName());
327   gDirectory->cd(GetName());
328
329   if (fData)
330   {
331     fData->Write();
332     AliPWG0Helper::CreateProjections(fData);
333   }
334   else
335     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
336
337   if (fDataUncorrected)
338   {
339     fDataUncorrected->Write();
340     AliPWG0Helper::CreateProjections(fDataUncorrected);
341   }
342   else
343     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fDataUncorrected is 0\n");
344
345   if (fVtx)
346     fVtx       ->Write();
347   else
348     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fVtx is 0\n");
349
350   if (fPtDist)
351     fPtDist       ->Write();
352   else
353     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
354
355   for (Int_t i=0; i<kVertexBinning; ++i)
356   {
357     if (fdNdEta[i])
358       fdNdEta[i]->Write();
359     else
360       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
361
362     if (fdNdEtaPtCutOffCorrected[i])
363       fdNdEtaPtCutOffCorrected[i]->Write();
364     else
365       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
366   }
367
368   gDirectory->cd("../");
369 }
370
371 void dNdEtaAnalysis::LoadHistograms()
372 {
373   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
374
375   gDirectory->cd(GetName());
376
377   fData = dynamic_cast<TH3F*> (gDirectory->Get(fData->GetName()));
378   fDataUncorrected = dynamic_cast<TH3F*> (gDirectory->Get(fDataUncorrected->GetName()));
379
380   fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
381
382   for (Int_t i=0; i<kVertexBinning; ++i)
383   {
384     fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
385     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
386   }
387
388   fPtDist = dynamic_cast<TH1D*> (gDirectory->Get(fPtDist->GetName()));
389
390   gDirectory->cd("../");
391 }
392
393 //____________________________________________________________________
394 void dNdEtaAnalysis::DrawHistograms()
395 {
396   // draws the histograms
397
398   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
399   canvas->Divide(2, 2);
400
401   canvas->cd(1);
402   if (fData)
403     fData->Draw("COLZ");
404
405   /*canvas->cd(2);
406   if (fDataUncorrected)
407     fDataUncorrected->Draw("COLZ");*/
408
409   canvas->cd(2);
410   if (fVtx)
411     fVtx->Draw();
412
413   canvas->cd(3);
414   if (fdNdEtaPtCutOffCorrected[0])
415     fdNdEtaPtCutOffCorrected[0]->Draw();
416
417   if (fdNdEta[0])
418   {
419     fdNdEta[0]->SetLineColor(kRed);
420     fdNdEta[0]->Draw("SAME");
421   }
422
423   canvas->cd(4);
424   if (fPtDist)
425     fPtDist->Draw();
426
427     // histograms for different vertices?
428   if (kVertexBinning > 0)
429   {
430       // doesnt work, but i dont get it, giving up...
431     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
432     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
433     //printf("%d\n", yPads);
434     //canvas2->Divide(2, yPads);
435
436     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
437
438     for (Int_t i=0; i<kVertexBinning; ++i)
439     {
440       //canvas2->cd(i-1);
441       //printf("%d\n", i);
442       if (fdNdEtaPtCutOffCorrected[i])
443       {
444         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
445         fdNdEtaPtCutOffCorrected[i]->Draw((i == 0) ? "" : "SAME");
446         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
447       }
448     }
449
450     legend->Draw();
451   }
452 }
453
454 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
455 {
456   // Merges a list of dNdEtaAnalysis objects with this one.
457   // This is needed for PROOF.
458   // Returns the number of merged objects (including this)
459
460   if (!list)
461     return 0;
462
463   if (list->IsEmpty())
464     return 1;
465
466   TIterator* iter = list->MakeIterator();
467   TObject* obj;
468
469   // sub collections
470   const Int_t nCollections = 2 * kVertexBinning + 4; // 4 standalone hists, two arrays of size kVertexBinning
471   TList* collections[nCollections];
472   for (Int_t i=0; i<nCollections; ++i)
473     collections[i] = new TList;
474
475   Int_t count = 0;
476   while ((obj = iter->Next()))
477   {
478     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
479     if (entry == 0)
480       continue;
481
482     collections[0]->Add(entry->fData);
483     if (fDataUncorrected)
484       collections[1]->Add(entry->fDataUncorrected);
485     collections[2]->Add(entry->fVtx);
486     collections[3]->Add(entry->fPtDist);
487
488     for (Int_t i=0; i<kVertexBinning; ++i)
489     {
490       collections[4+i]->Add(entry->fdNdEta[i]);
491       collections[4+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
492     }
493
494     ++count;
495   }
496
497   fData->Merge(collections[0]);
498   if (fDataUncorrected)
499     fDataUncorrected->Merge(collections[1]);
500   fVtx->Merge(collections[2]);
501   fPtDist->Merge(collections[3]);
502   for (Int_t i=0; i<kVertexBinning; ++i)
503   {
504     fdNdEta[i]->Merge(collections[4+i]);
505     fdNdEta[i]->Merge(collections[4+kVertexBinning+i]);
506   }
507
508   for (Int_t i=0; i<nCollections; ++i)
509     delete collections[i];
510
511   return count+1;
512 }