]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/dNdEta/dNdEtaAnalysis.cxx
o) more debug histograms
[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)
86 for (Int_t iVtx=0; iVtx<=hVtx->GetNbinsX(); iVtx++) {
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 }
114 }
115
116 // then take the wieghted average for each eta
117 // is this the right way to do it???
fcf2fb36 118 for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++) {
75ec0f41 119 Float_t sum = 0;
120 Float_t sumError2 = 0;
121 Int_t nMeasurements = 0;
122
fcf2fb36 123 Float_t sumXw = 0;
124 Float_t sumW = 0;
75ec0f41 125
7029240a 126 // do we have several histograms for different vertex positions?
fcf2fb36 127 Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
7029240a 128 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
129 {
fcf2fb36 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 }
7029240a 139
fcf2fb36 140 for (Int_t iVtx=vertexBinBegin; iVtx<=vertexBinEnd; iVtx++) {
141 if (hVtx->GetBinContent(iVtx)==0) continue;
a8d3b360 142 if (hEtaVsVtxCheck->GetBinContent(iVtx, iEta)==0) continue;
7029240a 143
fcf2fb36 144 Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
a8d3b360 145 sumXw = sumXw + hEtaVsVtxCheck->GetBinContent(iVtx, iEta)*w;
fcf2fb36 146 sumW = sumW + w;
7029240a 147
a8d3b360 148 sum = sum + hEtaVsVtxCheck->GetBinContent(iVtx, iEta);
149 sumError2 = sumError2 + TMath::Power(hEtaVsVtxCheck->GetBinError(iVtx, iEta),2);
7029240a 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
fcf2fb36 159 dndeta = sumXw/sumW;
160 error = 1/TMath::Sqrt(sumW);
7029240a 161
fcf2fb36 162 dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
163 error = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
7029240a 164
fcf2fb36 165 hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
166 hdNdEta[vertexPos]->SetBinError(iEta, error);
7029240a 167 }
75ec0f41 168 }
75ec0f41 169 }
170}
171
75ec0f41 172//____________________________________________________________________
173void
174dNdEtaAnalysis::SaveHistograms() {
175
7029240a 176 gDirectory->mkdir(GetName());
177 gDirectory->cd(GetName());
75ec0f41 178
179 hEtaVsVtx ->Write();
fcf2fb36 180 hEtaVsVtxUncorrected->Write();
75ec0f41 181 hVtx ->Write();
7029240a 182 for (Int_t i=0; i<kVertexBinning; ++i)
183 hdNdEta[i] ->Write();
75ec0f41 184
185 gDirectory->cd("../");
186}
187
ceb5d1b5 188//____________________________________________________________________
189void dNdEtaAnalysis::DrawHistograms()
190{
a8d3b360 191 TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 1200, 800);
192 canvas->Divide(3, 2);
ceb5d1b5 193
194 canvas->cd(1);
195 if (hEtaVsVtx)
196 hEtaVsVtx->Draw("COLZ");
197
198 canvas->cd(2);
fcf2fb36 199 if (hEtaVsVtxCheck)
200 hEtaVsVtxCheck->Draw("COLZ");
201
202 canvas->cd(3);
ceb5d1b5 203 if (hEtaVsVtxUncorrected)
204 hEtaVsVtxUncorrected->Draw("COLZ");
205
a8d3b360 206 canvas->cd(4);
207 TH2F* clone = (TH2F*) hEtaVsVtxCheck->Clone("clone");
208 clone->Divide(hEtaVsVtx);
209 clone->Draw("COLZ");
210
211 canvas->cd(5);
ceb5d1b5 212 if (hVtx)
a8d3b360 213 hVtx->Draw();
ceb5d1b5 214
a8d3b360 215 canvas->cd(6);
7029240a 216 if (hdNdEta[0])
217 hdNdEta[0]->Draw();
fcf2fb36 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 }
7029240a 244}
245
246Long64_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;
ceb5d1b5 294}