o) new way of calculating the final dndeta
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
1 /* $Id$ */
2
3 #include "dNdEtaAnalysis.h"
4
5 #include <TFile.h>
6 #include <TH2F.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 "dNdEtaCorrection.h"
16
17 //____________________________________________________________________
18 ClassImp(dNdEtaAnalysis)
19
20 //____________________________________________________________________
21 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
22 TNamed(name, title)
23 {
24   hEtaVsVtx  = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
25   hEtaVsVtx->SetXTitle("vtx z [cm]");
26   hEtaVsVtx->SetYTitle("#eta");
27
28   hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
29   hEtaVsVtxCheck = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_check", name)));
30   hVtx       = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
31   for (Int_t i=0; i<kVertexBinning; ++i)
32   {
33     hdNdEta[i]    = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
34     hdNdEta[i]->SetYTitle("dN/d#eta");
35   }
36
37   hEtaVsVtx->Sumw2();
38   hVtx->Sumw2();
39 }
40
41 //____________________________________________________________________
42 void
43 dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t c) {
44   hEtaVsVtxUncorrected->Fill(vtx,eta);
45   hEtaVsVtxCheck->Fill(vtx, eta, c);
46 }
47
48 //____________________________________________________________________
49 void
50 dNdEtaAnalysis::FillEvent(Float_t vtx) {
51   hVtx->Fill(vtx);
52 }
53
54 //____________________________________________________________________
55 void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction)
56 {
57   // correct with correction values if available
58   // TODO what do we do with the error?
59   if (!correction)
60     printf("INFO: No correction applied\n");
61
62   // this can be replaced by TH2F::Divide if we agree that the binning will be always the same
63   for (Int_t iVtx=0; iVtx<=hEtaVsVtxUncorrected->GetNbinsX(); iVtx++)
64   {
65     for (Int_t iEta=0; iEta<=hEtaVsVtxUncorrected->GetNbinsY(); iEta++)
66     {
67       Float_t correctionValue = 1;
68       if (correction)
69         correctionValue = correction->GetCorrection(hEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), hEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta));
70
71       Float_t value = hEtaVsVtxUncorrected->GetBinContent(iVtx, iEta);
72       Float_t error = hEtaVsVtxUncorrected->GetBinError(iVtx, iEta);
73
74       Float_t correctedValue = value * correctionValue;
75       Float_t correctedError = error * correctionValue;
76
77       if (correctedValue != 0)
78       {
79         hEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue);
80         hEtaVsVtx->SetBinError(iVtx, iEta, correctedError);
81       }
82     }
83   }
84
85   // normalize with n events (per vtx)
86 /*  for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) {
87     Float_t nEvents      = hVtx->GetBinContent(iVtx);
88     Float_t nEventsError = hVtx->GetBinError(iVtx);
89
90     if (nEvents==0) continue;
91
92     for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
93       Float_t value = hEtaVsVtx->GetBinContent(iVtx, iEta) / nEvents;
94       if (value==0) continue;
95       Float_t error = hEtaVsVtx->GetBinError(iVtx, iEta)/nEvents;
96       error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta)/
97                                        hEtaVsVtx->GetBinContent(iVtx, iEta),2) +
98                           TMath::Power(nEventsError/nEvents,2));
99       hEtaVsVtx->SetBinContent(iVtx, iEta, value);
100       hEtaVsVtx->SetBinError(iVtx,   iEta, error);
101     }
102
103     //debug
104     for (Int_t iEta=0; iEta<=hEtaVsVtxCheck->GetNbinsY(); iEta++) {
105       Float_t value = hEtaVsVtxCheck->GetBinContent(iVtx, iEta) / nEvents;
106       if (value==0) continue;
107       Float_t error = hEtaVsVtxCheck->GetBinError(iVtx, iEta)/nEvents;
108       error = TMath::Sqrt(TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta)/
109                                        hEtaVsVtxCheck->GetBinContent(iVtx, iEta),2) +
110                           TMath::Power(nEventsError/nEvents,2));
111       hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value);
112       hEtaVsVtxCheck->SetBinError(iVtx,   iEta, error);
113     }
114   }*/
115
116   // then take the wieghted average for each eta
117   // is this the right way to do it???
118   for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++)
119   {
120     // do we have several histograms for different vertex positions?
121     Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
122     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
123     {
124       Int_t vertexBinBegin = 1;
125       Int_t vertexBinEnd = hVtx->GetNbinsX() + 1;
126
127       // the first histogram is always for the whole vertex range
128       if (vertexPos > 0)
129       {
130         vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
131         vertexBinEnd = vertexBinBegin + vertexBinWidth;
132       }
133
134       Float_t totalEvents = hVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
135       if (totalEvents == 0)
136       {
137         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
138         continue;
139       }
140
141       Float_t sum = 0;
142       Float_t sumError2 = 0;
143       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
144       {
145         if (hEtaVsVtx->GetBinContent(iVtx, iEta) != 0)
146         {
147           sum = sum + hEtaVsVtx->GetBinContent(iVtx, iEta);
148           sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
149         }
150       }
151
152       Float_t dndeta = sum / totalEvents;
153       Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
154
155       dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
156       error  = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
157
158       hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
159       hdNdEta[vertexPos]->SetBinError(iEta, error);
160     }
161   }
162 }
163
164 //____________________________________________________________________
165 void
166 dNdEtaAnalysis::SaveHistograms() {
167
168   gDirectory->mkdir(GetName());
169   gDirectory->cd(GetName());
170   
171   hEtaVsVtx  ->Write();
172   hEtaVsVtxUncorrected->Write();
173   hVtx       ->Write();
174   for (Int_t i=0; i<kVertexBinning; ++i)
175     hdNdEta[i]    ->Write();
176
177   gDirectory->cd("../");
178 }
179
180 //____________________________________________________________________
181 void dNdEtaAnalysis::DrawHistograms()
182 {
183   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 1200, 800);
184   canvas->Divide(3, 2);
185
186   canvas->cd(1);
187   if (hEtaVsVtx)
188     hEtaVsVtx->Draw("COLZ");
189
190   canvas->cd(2);
191   if (hEtaVsVtxCheck)
192     hEtaVsVtxCheck->Draw("COLZ");
193
194   canvas->cd(3);
195   if (hEtaVsVtxUncorrected)
196     hEtaVsVtxUncorrected->Draw("COLZ");
197
198   canvas->cd(4);
199   TH2F* clone = (TH2F*) hEtaVsVtxCheck->Clone("clone");
200   clone->Divide(hEtaVsVtx);
201   clone->Draw("COLZ");
202
203   canvas->cd(5);
204   if (hVtx)
205     hVtx->Draw();
206
207   canvas->cd(6);
208   if (hdNdEta[0])
209     hdNdEta[0]->Draw();
210
211     // histograms for different vertices?
212   if (kVertexBinning > 0)
213   {
214       // doesnt work, but i dont get it, giving up...
215     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
216     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
217     //printf("%d\n", yPads);
218     //canvas2->Divide(2, yPads);
219
220     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
221
222     for (Int_t i=0; i<kVertexBinning; ++i)
223     {
224       //canvas2->cd(i-1);
225       //printf("%d\n", i);
226       if (hdNdEta[i])
227       {
228         hdNdEta[i]->SetLineColor(i+1);
229         hdNdEta[i]->Draw((i == 0) ? "" : "SAME");
230         legend->AddEntry(hdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
231       }
232     }
233
234     legend->Draw();
235   }
236 }
237
238 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
239 {
240   // Merges a list of dNdEtaAnalysis objects with this one.
241   // This is needed for PROOF.
242   // Returns the number of merged objects (including this)
243
244   if (!list)
245     return 0;
246
247   if (list->IsEmpty())
248     return 1;
249
250   TIterator* iter = list->MakeIterator();
251   TObject* obj;
252
253   // sub collections
254   const Int_t nCollections = kVertexBinning + 3;
255   TList* collections[nCollections];
256   for (Int_t i=0; i<nCollections; ++i)
257     collections[i] = new TList;
258
259   Int_t count = 0;
260   while ((obj = iter->Next()))
261   {
262     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
263     if (entry == 0)
264       continue;
265
266     collections[0]->Add(entry->hEtaVsVtx);
267     collections[1]->Add(entry->hEtaVsVtxUncorrected);
268     collections[2]->Add(entry->hVtx);
269
270     for (Int_t i=0; i<kVertexBinning; ++i)
271       collections[3+i]->Add(entry->hdNdEta[i]);
272
273     ++count;
274   }
275
276   hEtaVsVtx->Merge(collections[0]);
277   hEtaVsVtxUncorrected->Merge(collections[1]);
278   hVtx->Merge(collections[2]);
279   for (Int_t i=0; i<kVertexBinning; ++i)
280     hdNdEta[i]->Merge(collections[3+i]);
281
282   for (Int_t i=0; i<nCollections; ++i)
283     delete collections[i];
284
285   return count+1;
286 }