implemented analysis using 3d information (vtx_z, eta, pt)
[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(Char_t* name, Char_t* title) :
24   TNamed(name, title),
25   fData(0),
26   fDataUncorrected(0),
27   fNEvents(0),
28   fVtx(0)
29 {
30   // constructor
31
32   fData  = new TH3F(Form("%s_analysis", name),"dNdEtaAnalysis",80,-20,20,120,-6,6,100, 0, 10);
33   fData->SetXTitle("vtx z [cm]");
34   fData->SetYTitle("#eta");
35   fData->SetZTitle("p_{T}");
36
37   fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
38   fVtx       = dynamic_cast<TH1D*> (fData->Project3D("x"));
39
40   fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
41   for (Int_t i=1; i<kVertexBinning; ++i)
42   {
43     fdNdEta[i]    = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
44     fdNdEta[i]->SetYTitle("dN/d#eta");
45   }
46
47   fData->Sumw2();
48   fVtx->Sumw2();
49 }
50
51 //____________________________________________________________________
52 dNdEtaAnalysis::~dNdEtaAnalysis()
53 {
54   // destructor
55
56   delete fData;
57   fData = 0;
58
59   delete fDataUncorrected;
60   fDataUncorrected = 0;
61
62   delete fVtx;
63   fVtx = 0;
64
65   for (Int_t i=0; i<kVertexBinning; ++i)
66   {
67     delete fdNdEta[i];
68     fdNdEta[i] = 0;
69   }
70 }
71
72 //_____________________________________________________________________________
73 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
74   TNamed(c),
75   fData(0),
76   fDataUncorrected(0),
77   fNEvents(0),
78   fVtx(0)
79 {
80   //
81   // dNdEtaAnalysis copy constructor
82   //
83
84   ((dNdEtaAnalysis &) c).Copy(*this);
85 }
86
87 //_____________________________________________________________________________
88 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
89 {
90   //
91   // Assignment operator
92   //
93
94   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
95   return *this;
96 }
97
98 //_____________________________________________________________________________
99 void dNdEtaAnalysis::Copy(TObject &c) const
100 {
101   //
102   // Copy function
103   //
104
105   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
106
107   target.fData = dynamic_cast<TH3F*> (fData->Clone());
108   target.fDataUncorrected = dynamic_cast<TH3F*> (fDataUncorrected->Clone());
109   target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
110
111   for (Int_t i=0; i<kVertexBinning; ++i)
112     target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
113
114   target.fNEvents = fNEvents;
115
116   TNamed::Copy((TNamed &) c);
117 }
118
119 //____________________________________________________________________
120 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight)
121 {
122   // fills a track into the histograms
123
124   fDataUncorrected->Fill(vtx, eta, pt);
125   fData->Fill(vtx, eta, pt, weight);
126 }
127
128 //____________________________________________________________________
129 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
130 {
131   // fills an event into the histograms
132
133   fVtx->Fill(vtx, weight);
134
135   fNEvents += weight;
136 }
137
138 //____________________________________________________________________
139 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
140 {
141   // correct with correction values if available
142
143   // TODO what do we do with the error?
144   if (!correction)
145     printf("INFO: No correction applied\n");
146
147   // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
148
149   // integrate over pt (with pt cut)
150   fData->GetZaxis()->SetRange(fData->GetZaxis()->FindBin(ptCut), fData->GetZaxis()->GetNbins());
151   TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2"));
152   if (vtxVsEta == 0)
153   {
154     printf("ERROR: pt integration failed\n");
155     return;
156   }
157
158   new TCanvas;
159   vtxVsEta->Draw("COLZ");
160
161   for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
162   {
163     // do we have several histograms for different vertex positions?
164     Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
165     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
166     {
167       Int_t vertexBinBegin = 1;
168       Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
169
170       // the first histogram is always for the whole vertex range
171       if (vertexPos > 0)
172       {
173         vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
174         vertexBinEnd = vertexBinBegin + vertexBinWidth;
175       }
176
177       Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
178       if (totalEvents == 0)
179       {
180         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
181         continue;
182       }
183
184       Float_t sum = 0;
185       Float_t sumError2 = 0;
186       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
187       {
188         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
189         {
190           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
191           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
192         }
193       }
194
195       Float_t ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
196       //ptCutOffCorrection = 1;
197       if (ptCutOffCorrection <= 0)
198       {
199         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
200         continue;
201       }
202
203       Float_t dndeta = sum / totalEvents / ptCutOffCorrection;
204       Float_t error  = TMath::Sqrt(sumError2) / totalEvents / ptCutOffCorrection;
205
206       dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
207       error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
208
209       fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
210       fdNdEta[vertexPos]->SetBinError(iEta, error);
211     }
212   }
213 }
214
215 //____________________________________________________________________
216 void dNdEtaAnalysis::SaveHistograms()
217 {
218   // save the histograms to a directory with the name of this class (retrieved from TNamed)
219
220   gDirectory->mkdir(GetName());
221   gDirectory->cd(GetName());
222
223   fData  ->Write();
224   AliPWG0Helper::CreateProjections(fData);
225   fDataUncorrected->Write();
226   AliPWG0Helper::CreateProjections(fDataUncorrected);
227
228   fVtx       ->Write();
229   for (Int_t i=0; i<kVertexBinning; ++i)
230     fdNdEta[i]    ->Write();
231
232   gDirectory->cd("../");
233 }
234
235 void dNdEtaAnalysis::LoadHistograms()
236 {
237   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
238
239   gDirectory->cd(GetName());
240
241   fData = dynamic_cast<TH3F*> (gDirectory->Get(fData->GetName()));
242   fDataUncorrected = dynamic_cast<TH3F*> (gDirectory->Get(fDataUncorrected->GetName()));
243
244   fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
245
246   for (Int_t i=0; i<kVertexBinning; ++i)
247     fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
248
249   gDirectory->cd("../");
250 }
251
252 //____________________________________________________________________
253 void dNdEtaAnalysis::DrawHistograms()
254 {
255   // draws the histograms
256
257   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
258   canvas->Divide(2, 2);
259
260   canvas->cd(1);
261   if (fData)
262     fData->Draw("COLZ");
263
264   canvas->cd(2);
265   if (fDataUncorrected)
266     fDataUncorrected->Draw("COLZ");
267
268   canvas->cd(3);
269   if (fVtx)
270     fVtx->Draw();
271
272   canvas->cd(4);
273   if (fdNdEta[0])
274     fdNdEta[0]->Draw();
275
276     // histograms for different vertices?
277   if (kVertexBinning > 0)
278   {
279       // doesnt work, but i dont get it, giving up...
280     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
281     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
282     //printf("%d\n", yPads);
283     //canvas2->Divide(2, yPads);
284
285     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
286
287     for (Int_t i=0; i<kVertexBinning; ++i)
288     {
289       //canvas2->cd(i-1);
290       //printf("%d\n", i);
291       if (fdNdEta[i])
292       {
293         fdNdEta[i]->SetLineColor(i+1);
294         fdNdEta[i]->Draw((i == 0) ? "" : "SAME");
295         legend->AddEntry(fdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
296       }
297     }
298
299     legend->Draw();
300   }
301 }
302
303 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
304 {
305   // Merges a list of dNdEtaAnalysis objects with this one.
306   // This is needed for PROOF.
307   // Returns the number of merged objects (including this)
308
309   if (!list)
310     return 0;
311
312   if (list->IsEmpty())
313     return 1;
314
315   TIterator* iter = list->MakeIterator();
316   TObject* obj;
317
318   // sub collections
319   const Int_t nCollections = kVertexBinning + 3;
320   TList* collections[nCollections];
321   for (Int_t i=0; i<nCollections; ++i)
322     collections[i] = new TList;
323
324   Int_t count = 0;
325   while ((obj = iter->Next()))
326   {
327     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
328     if (entry == 0)
329       continue;
330
331     collections[0]->Add(entry->fData);
332     collections[1]->Add(entry->fDataUncorrected);
333     collections[2]->Add(entry->fVtx);
334
335     for (Int_t i=0; i<kVertexBinning; ++i)
336       collections[3+i]->Add(entry->fdNdEta[i]);
337
338     ++count;
339   }
340
341   fData->Merge(collections[0]);
342   fDataUncorrected->Merge(collections[1]);
343   fVtx->Merge(collections[2]);
344   for (Int_t i=0; i<kVertexBinning; ++i)
345     fdNdEta[i]->Merge(collections[3+i]);
346
347   for (Int_t i=0; i<nCollections; ++i)
348     delete collections[i];
349
350   return count+1;
351 }