ecb6f7bf67d97a949731072eeedee875978c25b8
[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     Float_t sum           = 0;
120     Float_t sumError2     = 0;
121     Int_t   nMeasurements = 0;    
122
123     Float_t sumXw = 0;
124     Float_t sumW  = 0;
125     
126     // do we have several histograms for different vertex positions?
127     Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
128     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
129     {
130       Int_t vertexBinBegin = 0;
131       Int_t vertexBinEnd = hVtx->GetNbinsX();
132
133       // the first histogram is always for the whole vertex range
134       if (vertexPos > 0)
135       {
136         vertexBinBegin = vertexBinWidth * (vertexPos-1);
137         vertexBinEnd = vertexBinBegin + vertexBinWidth;
138       }
139
140       for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) {
141         if (hVtx->GetBinContent(iVtx)==0)             continue;
142         if (hEtaVsVtxCheck->GetBinContent(iVtx, iEta)==0) continue;
143
144         Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
145         sumXw = sumXw + hEtaVsVtxCheck->GetBinContent(iVtx, iEta)*w;
146         sumW  = sumW + w;
147
148         sum = sum + hEtaVsVtxCheck->GetBinContent(iVtx, iEta);
149         sumError2 = sumError2 + TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta),2);      
150         nMeasurements++;
151       }
152       Float_t dndeta = 0;
153       Float_t error  = 0;
154
155       if (nMeasurements!=0) {
156         dndeta = sum/Float_t(nMeasurements);
157         error  = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
158
159         dndeta = sumXw/sumW;
160         error  = 1/TMath::Sqrt(sumW);
161
162         dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
163         error  = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
164
165         hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
166         hdNdEta[vertexPos]->SetBinError(iEta, error);
167       }
168     }
169   }
170 }
171
172 //____________________________________________________________________
173 void
174 dNdEtaAnalysis::SaveHistograms() {
175
176   gDirectory->mkdir(GetName());
177   gDirectory->cd(GetName());
178   
179   hEtaVsVtx  ->Write();
180   hEtaVsVtxUncorrected->Write();
181   hVtx       ->Write();
182   for (Int_t i=0; i<kVertexBinning; ++i)
183     hdNdEta[i]    ->Write();
184
185   gDirectory->cd("../");
186 }
187
188 //____________________________________________________________________
189 void dNdEtaAnalysis::DrawHistograms()
190 {
191   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 1200, 800);
192   canvas->Divide(3, 2);
193
194   canvas->cd(1);
195   if (hEtaVsVtx)
196     hEtaVsVtx->Draw("COLZ");
197
198   canvas->cd(2);
199   if (hEtaVsVtxCheck)
200     hEtaVsVtxCheck->Draw("COLZ");
201
202   canvas->cd(3);
203   if (hEtaVsVtxUncorrected)
204     hEtaVsVtxUncorrected->Draw("COLZ");
205
206   canvas->cd(4);
207   TH2F* clone = (TH2F*) hEtaVsVtxCheck->Clone("clone");
208   clone->Divide(hEtaVsVtx);
209   clone->Draw("COLZ");
210
211   canvas->cd(5);
212   if (hVtx)
213     hVtx->Draw();
214
215   canvas->cd(6);
216   if (hdNdEta[0])
217     hdNdEta[0]->Draw();
218
219     // histograms for different vertices?
220   if (kVertexBinning > 0)
221   {
222       // doesnt work, but i dont get it, giving up...
223     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
224     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
225     //printf("%d\n", yPads);
226     //canvas2->Divide(2, yPads);
227
228     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
229
230     for (Int_t i=1; i<kVertexBinning; ++i)
231     {
232       //canvas2->cd(i-1);
233       //printf("%d\n", i);
234       if (hdNdEta[i])
235       {
236         hdNdEta[i]->SetLineColor(i);
237         hdNdEta[i]->Draw((i == 1) ? "" : "SAME");
238         legend->AddEntry(hdNdEta[i], Form("Vtx Bin %d", i-1));
239       }
240     }
241
242     legend->Draw();
243   }
244 }
245
246 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
247 {
248   // Merges a list of dNdEtaAnalysis objects with this one.
249   // This is needed for PROOF.
250   // Returns the number of merged objects (including this)
251
252   if (!list)
253     return 0;
254
255   if (list->IsEmpty())
256     return 1;
257
258   TIterator* iter = list->MakeIterator();
259   TObject* obj;
260
261   // sub collections
262   const Int_t nCollections = kVertexBinning + 3;
263   TList* collections[nCollections];
264   for (Int_t i=0; i<nCollections; ++i)
265     collections[i] = new TList;
266
267   Int_t count = 0;
268   while ((obj = iter->Next()))
269   {
270     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
271     if (entry == 0)
272       continue;
273
274     collections[0]->Add(entry->hEtaVsVtx);
275     collections[1]->Add(entry->hEtaVsVtxUncorrected);
276     collections[2]->Add(entry->hVtx);
277
278     for (Int_t i=0; i<kVertexBinning; ++i)
279       collections[3+i]->Add(entry->hdNdEta[i]);
280
281     ++count;
282   }
283
284   hEtaVsVtx->Merge(collections[0]);
285   hEtaVsVtxUncorrected->Merge(collections[1]);
286   hVtx->Merge(collections[2]);
287   for (Int_t i=0; i<kVertexBinning; ++i)
288     hdNdEta[i]->Merge(collections[3+i]);
289
290   for (Int_t i=0; i<nCollections; ++i)
291     delete collections[i];
292
293   return count+1;
294 }