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