o) new way of calculating the final dndeta
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
CommitLineData
dc740de4 1/* $Id$ */
2
75ec0f41 3#include "dNdEtaAnalysis.h"
4
ceb5d1b5 5#include <TFile.h>
6#include <TH2F.h>
7#include <TH1D.h>
8#include <TMath.h>
9#include <TCanvas.h>
7029240a 10#include <TCollection.h>
11#include <TIterator.h>
12#include <TList.h>
fcf2fb36 13#include <TLegend.h>
14
15#include "dNdEtaCorrection.h"
ceb5d1b5 16
75ec0f41 17//____________________________________________________________________
b7f4a1fd 18ClassImp(dNdEtaAnalysis)
75ec0f41 19
20//____________________________________________________________________
7029240a 21dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
22TNamed(name, title)
23{
24 hEtaVsVtx = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
75ec0f41 25 hEtaVsVtx->SetXTitle("vtx z [cm]");
26 hEtaVsVtx->SetYTitle("#eta");
4dd2ad81 27
7029240a 28 hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
fcf2fb36 29 hEtaVsVtxCheck = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_check", name)));
7029240a 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 }
75ec0f41 36
37 hEtaVsVtx->Sumw2();
38 hVtx->Sumw2();
75ec0f41 39}
40
41//____________________________________________________________________
42void
fcf2fb36 43dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t c) {
ceb5d1b5 44 hEtaVsVtxUncorrected->Fill(vtx,eta);
fcf2fb36 45 hEtaVsVtxCheck->Fill(vtx, eta, c);
75ec0f41 46}
47
48//____________________________________________________________________
49void
50dNdEtaAnalysis::FillEvent(Float_t vtx) {
51 hVtx->Fill(vtx);
52}
53
54//____________________________________________________________________
fcf2fb36 55void 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)
5af55649 86/* for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) {
fcf2fb36 87 Float_t nEvents = hVtx->GetBinContent(iVtx);
88 Float_t nEventsError = hVtx->GetBinError(iVtx);
75ec0f41 89
75ec0f41 90 if (nEvents==0) continue;
fcf2fb36 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;
75ec0f41 106 if (value==0) continue;
fcf2fb36 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) +
75ec0f41 110 TMath::Power(nEventsError/nEvents,2));
fcf2fb36 111 hEtaVsVtxCheck->SetBinContent(iVtx, iEta, value);
112 hEtaVsVtxCheck->SetBinError(iVtx, iEta, error);
75ec0f41 113 }
5af55649 114 }*/
75ec0f41 115
116 // then take the wieghted average for each eta
117 // is this the right way to do it???
5af55649 118 for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++)
119 {
7029240a 120 // do we have several histograms for different vertex positions?
fcf2fb36 121 Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
7029240a 122 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
123 {
5af55649 124 Int_t vertexBinBegin = 1;
125 Int_t vertexBinEnd = hVtx->GetNbinsX() + 1;
fcf2fb36 126
127 // the first histogram is always for the whole vertex range
128 if (vertexPos > 0)
129 {
5af55649 130 vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
fcf2fb36 131 vertexBinEnd = vertexBinBegin + vertexBinWidth;
132 }
7029240a 133
5af55649 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;
7029240a 139 }
7029240a 140
5af55649 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 }
7029240a 151
5af55649 152 Float_t dndeta = sum / totalEvents;
153 Float_t error = TMath::Sqrt(sumError2) / totalEvents;
7029240a 154
5af55649 155 dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
156 error = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
7029240a 157
5af55649 158 hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
159 hdNdEta[vertexPos]->SetBinError(iEta, error);
75ec0f41 160 }
75ec0f41 161 }
162}
163
75ec0f41 164//____________________________________________________________________
165void
166dNdEtaAnalysis::SaveHistograms() {
167
7029240a 168 gDirectory->mkdir(GetName());
169 gDirectory->cd(GetName());
75ec0f41 170
171 hEtaVsVtx ->Write();
fcf2fb36 172 hEtaVsVtxUncorrected->Write();
75ec0f41 173 hVtx ->Write();
7029240a 174 for (Int_t i=0; i<kVertexBinning; ++i)
175 hdNdEta[i] ->Write();
75ec0f41 176
177 gDirectory->cd("../");
178}
179
ceb5d1b5 180//____________________________________________________________________
181void dNdEtaAnalysis::DrawHistograms()
182{
a8d3b360 183 TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 1200, 800);
184 canvas->Divide(3, 2);
ceb5d1b5 185
186 canvas->cd(1);
187 if (hEtaVsVtx)
188 hEtaVsVtx->Draw("COLZ");
189
190 canvas->cd(2);
fcf2fb36 191 if (hEtaVsVtxCheck)
192 hEtaVsVtxCheck->Draw("COLZ");
193
194 canvas->cd(3);
ceb5d1b5 195 if (hEtaVsVtxUncorrected)
196 hEtaVsVtxUncorrected->Draw("COLZ");
197
a8d3b360 198 canvas->cd(4);
199 TH2F* clone = (TH2F*) hEtaVsVtxCheck->Clone("clone");
200 clone->Divide(hEtaVsVtx);
201 clone->Draw("COLZ");
202
203 canvas->cd(5);
ceb5d1b5 204 if (hVtx)
a8d3b360 205 hVtx->Draw();
ceb5d1b5 206
a8d3b360 207 canvas->cd(6);
7029240a 208 if (hdNdEta[0])
209 hdNdEta[0]->Draw();
fcf2fb36 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
5af55649 222 for (Int_t i=0; i<kVertexBinning; ++i)
fcf2fb36 223 {
224 //canvas2->cd(i-1);
225 //printf("%d\n", i);
226 if (hdNdEta[i])
227 {
5af55649 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));
fcf2fb36 231 }
232 }
233
234 legend->Draw();
235 }
7029240a 236}
237
238Long64_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;
ceb5d1b5 286}