]>
Commit | Line | Data |
---|---|---|
1 | /* $Id$ */ | |
2 | ||
3 | #include "dNdEtaAnalysis.h" | |
4 | ||
5 | #include <TFile.h> | |
6 | #include <TH3F.h> | |
7 | #include <TH2F.h> | |
8 | #include <TH1F.h> | |
9 | #include <TMath.h> | |
10 | #include <TCanvas.h> | |
11 | #include <TCollection.h> | |
12 | #include <TIterator.h> | |
13 | #include <TList.h> | |
14 | #include <TLegend.h> | |
15 | #include <TLine.h> | |
16 | #include <TParameter.h> | |
17 | ||
18 | #include <AliLog.h> | |
19 | #include "AlidNdEtaCorrection.h" | |
20 | #include <AliCorrection.h> | |
21 | #include <AliPWG0Helper.h> | |
22 | #include <AliCorrectionMatrix2D.h> | |
23 | #include <AliCorrectionMatrix3D.h> | |
24 | ||
25 | //____________________________________________________________________ | |
26 | ClassImp(dNdEtaAnalysis) | |
27 | ||
28 | //____________________________________________________________________ | |
29 | dNdEtaAnalysis::dNdEtaAnalysis() : | |
30 | TNamed(), | |
31 | fData(0), | |
32 | fMult(0), | |
33 | fPtDist(0), | |
34 | fAnalysisMode(AliPWG0Helper::kInvalid), | |
35 | fTag(), | |
36 | fvtxMin(-9.99), | |
37 | fvtxMax(9.99) | |
38 | { | |
39 | // default constructor | |
40 | ||
41 | for (Int_t i=0; i<kVertexBinning; ++i) | |
42 | { | |
43 | fdNdEta[i] = 0; | |
44 | fdNdEtaPtCutOffCorrected[i] = 0; | |
45 | } | |
46 | } | |
47 | ||
48 | //____________________________________________________________________ | |
49 | dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) : | |
50 | TNamed(name, title), | |
51 | fData(0), | |
52 | fMult(0), | |
53 | fPtDist(0), | |
54 | fAnalysisMode(analysisMode), | |
55 | fTag(), | |
56 | fvtxMin(-9.99), | |
57 | fvtxMax(9.99) | |
58 | { | |
59 | // constructor | |
60 | ||
61 | fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode); | |
62 | ||
63 | // do not add this hists to the directory | |
64 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
65 | TH1::AddDirectory(kFALSE); | |
66 | ||
67 | fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5); | |
68 | ||
69 | TH1* histForBinning = fData->GetTrackCorrection()->GetGeneratedHistogram(); | |
70 | fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", histForBinning->GetNbinsY(), histForBinning->GetYaxis()->GetXbins()->GetArray()); | |
71 | ||
72 | fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName()))); | |
73 | ||
74 | for (Int_t i=1; i<kVertexBinning; ++i) | |
75 | { | |
76 | fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i))); | |
77 | fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i))); | |
78 | } | |
79 | ||
80 | fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", histForBinning->GetNbinsZ(), histForBinning->GetZaxis()->GetXbins()->GetArray()); | |
81 | ||
82 | TH1::AddDirectory(oldStatus); | |
83 | } | |
84 | ||
85 | //____________________________________________________________________ | |
86 | dNdEtaAnalysis::~dNdEtaAnalysis() | |
87 | { | |
88 | // destructor | |
89 | ||
90 | if (fData) | |
91 | { | |
92 | delete fData; | |
93 | fData = 0; | |
94 | } | |
95 | ||
96 | if (fMult) | |
97 | { | |
98 | delete fMult; | |
99 | fMult = 0; | |
100 | } | |
101 | ||
102 | for (Int_t i=0; i<kVertexBinning; ++i) | |
103 | { | |
104 | if (fdNdEta[i]) | |
105 | { | |
106 | delete fdNdEta[i]; | |
107 | fdNdEta[i] = 0; | |
108 | } | |
109 | if (fdNdEtaPtCutOffCorrected[i]) | |
110 | { | |
111 | delete fdNdEtaPtCutOffCorrected[i]; | |
112 | fdNdEtaPtCutOffCorrected[i] = 0; | |
113 | } | |
114 | } | |
115 | ||
116 | if (fPtDist) | |
117 | { | |
118 | delete fPtDist; | |
119 | fPtDist = 0; | |
120 | } | |
121 | } | |
122 | ||
123 | //_____________________________________________________________________________ | |
124 | dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) : | |
125 | TNamed(c), | |
126 | fData(0), | |
127 | fMult(0), | |
128 | fPtDist(0), | |
129 | fAnalysisMode(AliPWG0Helper::kInvalid), | |
130 | fTag(), | |
131 | fvtxMin(-9.99), | |
132 | fvtxMax(9.99) | |
133 | { | |
134 | // | |
135 | // dNdEtaAnalysis copy constructor | |
136 | // | |
137 | ||
138 | ((dNdEtaAnalysis &) c).Copy(*this); | |
139 | } | |
140 | ||
141 | //_____________________________________________________________________________ | |
142 | dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c) | |
143 | { | |
144 | // | |
145 | // Assignment operator | |
146 | // | |
147 | ||
148 | if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this); | |
149 | return *this; | |
150 | } | |
151 | ||
152 | //_____________________________________________________________________________ | |
153 | void dNdEtaAnalysis::Copy(TObject &c) const | |
154 | { | |
155 | // | |
156 | // Copy function | |
157 | // | |
158 | ||
159 | dNdEtaAnalysis& target = (dNdEtaAnalysis &) c; | |
160 | ||
161 | target.fData = dynamic_cast<AliCorrection*> (fData->Clone()); | |
162 | target.fMult = dynamic_cast<TH1F*> (fMult->Clone()); | |
163 | ||
164 | for (Int_t i=0; i<kVertexBinning; ++i) | |
165 | { | |
166 | target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone()); | |
167 | target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone()); | |
168 | } | |
169 | ||
170 | target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone()); | |
171 | ||
172 | target.fAnalysisMode = fAnalysisMode; | |
173 | target.fTag = fTag; | |
174 | ||
175 | TNamed::Copy((TNamed &) c); | |
176 | } | |
177 | ||
178 | //____________________________________________________________________ | |
179 | void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt) | |
180 | { | |
181 | // fills a track into the histograms | |
182 | ||
183 | fData->GetTrackCorrection()->FillMeas(vtx, eta, pt); | |
184 | } | |
185 | ||
186 | //____________________________________________________________________ | |
187 | void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n) | |
188 | { | |
189 | // fills an event into the histograms | |
190 | ||
191 | fData->GetEventCorrection()->FillMeas(vtx, n); | |
192 | } | |
193 | ||
194 | //____________________________________________________________________ | |
195 | void dNdEtaAnalysis::FillTriggeredEvent(Float_t n) | |
196 | { | |
197 | // fills a triggered event into the histograms | |
198 | ||
199 | fMult->Fill(n); | |
200 | } | |
201 | ||
202 | //____________________________________________________________________ | |
203 | void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag, Int_t backgroundSubtraction, TH1* combinatoricsCorrection) | |
204 | { | |
205 | // | |
206 | // correct with the given correction values and calculate dNdEta and pT distribution | |
207 | // the corrections that are applied can be steered by the flag correctionType | |
208 | // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!) | |
209 | // | |
210 | ||
211 | fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut); | |
212 | Printf("\n\n%s", fTag.Data()); | |
213 | ||
214 | if (combinatoricsCorrection) | |
215 | Printf("Combinatorics subtraction active!"); | |
216 | ||
217 | // set corrections to 1 | |
218 | fData->SetCorrectionToUnity(); | |
219 | ||
220 | if (correction && correctionType != AlidNdEtaCorrection::kNone) | |
221 | { | |
222 | TH3* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram(); | |
223 | TH2* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram(); | |
224 | ||
225 | if (correctionType >= AlidNdEtaCorrection::kTrack2Particle) | |
226 | trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram()); | |
227 | ||
228 | if (correctionType >= AlidNdEtaCorrection::kVertexReco) | |
229 | { | |
230 | trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram()); | |
231 | eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram()); | |
232 | ||
233 | // set bin with multiplicity 0 to unity (correction has no meaning in this bin) | |
234 | for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++) | |
235 | eventCorr->SetBinContent(i, 1, 1); | |
236 | } | |
237 | ||
238 | switch (correctionType) | |
239 | { | |
240 | case AlidNdEtaCorrection::kINEL : | |
241 | { | |
242 | trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram()); | |
243 | eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram()); | |
244 | break; | |
245 | } | |
246 | case AlidNdEtaCorrection::kNSD : | |
247 | { | |
248 | trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram()); | |
249 | eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram()); | |
250 | break; | |
251 | } | |
252 | case AlidNdEtaCorrection::kND : | |
253 | { | |
254 | trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram()); | |
255 | eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram()); | |
256 | break; | |
257 | } | |
258 | case AlidNdEtaCorrection::kOnePart : | |
259 | { | |
260 | trackCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetTrackCorrection()->GetCorrectionHistogram()); | |
261 | eventCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->GetCorrectionHistogram()); | |
262 | break; | |
263 | } | |
264 | default : break; | |
265 | } | |
266 | } | |
267 | else | |
268 | printf("INFO: No correction applied\n"); | |
269 | ||
270 | TH2F* rawMeasured = (TH2F*) fData->GetEventCorrection()->GetMeasuredHistogram()->Clone("rawMeasured"); | |
271 | ||
272 | fData->ResetErrorsOnCorrections(); | |
273 | fData->Multiply(); | |
274 | ||
275 | if (correction && correctionType >= AlidNdEtaCorrection::kVertexReco) | |
276 | { | |
277 | // There are no events with vertex that have 0 multiplicity, therefore | |
278 | // populate bin with 0 multiplicity with the following idea: | |
279 | // alpha = triggered events with vertex at a given vertex position / all triggered events with vertex | |
280 | // triggered events without vertex and 0 multiplicity at a given vertex position = alpha * all triggered events with 0 multiplicity | |
281 | // afterwards we still correct for the trigger efficiency | |
282 | // at the same time calculate expectation from MC (not used, just to check the value) | |
283 | ||
284 | //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram(); | |
285 | TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram(); | |
286 | ||
287 | TH2* eAll = correction->GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram(); | |
288 | TH2* eTrig = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram(); | |
289 | TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram(); | |
290 | TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsY()+1); | |
291 | ||
292 | //new TCanvas; correctedEvents->DrawCopy("TEXT"); | |
293 | ||
294 | // start above 0 mult. bin with integration | |
295 | TH1* vertexDist = rawMeasured->ProjectionX("vertexdist_measured", 2, rawMeasured->GetNbinsY()+1); | |
296 | //new TCanvas; vertexDist->DrawCopy(); | |
297 | ||
298 | Int_t allEventsWithVertex = (Int_t) vertexDist->Integral(0, vertexDist->GetNbinsX()+1); // include under/overflow! | |
299 | Int_t triggeredEventsWith0Mult = (Int_t) fMult->GetBinContent(1); | |
300 | ||
301 | Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex); | |
302 | ||
303 | if (backgroundSubtraction > 0) | |
304 | { | |
305 | triggeredEventsWith0Mult -= backgroundSubtraction; | |
306 | Printf("Subtracted %d background events from 0 mult. bin", backgroundSubtraction); | |
307 | } | |
308 | ||
309 | TH1* kineBias = (TH1*) vertexDist->Clone("kineBias"); | |
310 | kineBias->Reset(); | |
311 | ||
312 | // loop over vertex bins | |
313 | for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++) | |
314 | { | |
315 | Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex; | |
316 | Double_t events = alpha * triggeredEventsWith0Mult; | |
317 | ||
318 | if (eTrigVtx_projx->GetBinContent(i) == 0) | |
319 | continue; | |
320 | ||
321 | // calculate how many events we would have got with a pure MC-based correction | |
322 | // in the given bin: measured events with vertex (mult > 0) * triggered events with mult 0 (mc) / events with vertex and mult > 0 (mc) * trigger correction for bin 0 | |
323 | Printf("+++ 0-Bin Correction for bin %d +++", i); | |
324 | Printf(" Events: %f", vertexDist->GetBinContent(i)); | |
325 | Printf(" Ratio triggered N==0 / triggered vertex N>0: %f", eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i)); | |
326 | Printf(" Ratio all N==0 / triggered vertex N>0: %f", eAll->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i)); | |
327 | Printf(" Correction factor: %f", fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1)); | |
328 | ||
329 | //Double_t mcEvents = vertexDist->GetBinContent(i) * eTrig->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i) * fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1); | |
330 | Double_t mcEvents = vertexDist->GetBinContent(i) * eAll->GetBinContent(i, 1) / eTrigVtx_projx->GetBinContent(i); | |
331 | if (backgroundSubtraction == -1) | |
332 | { | |
333 | Printf("Using MC value for 0-bin correction!"); | |
334 | events = mcEvents; | |
335 | } | |
336 | else | |
337 | { | |
338 | Double_t fZ = eTrigVtx_projx->Integral(0, eTrigVtx_projx->GetNbinsX()+1) / eTrigVtx_projx->GetBinContent(i) * | |
339 | eTrig->GetBinContent(i, 1) / eTrig->Integral(0, eTrig->GetNbinsX()+1, 1, 1); | |
340 | ||
341 | kineBias->SetBinContent(i, fZ); | |
342 | ||
343 | events *= fZ; | |
344 | ||
345 | // multiply with trigger correction if set above | |
346 | events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1); | |
347 | ||
348 | Printf(" Bin %d, alpha is %.2f%%, fZ is %.3f, number of events with 0 mult.: %.2f (MC comparison: %.2f)", i, alpha * 100., fZ, events, mcEvents); | |
349 | } | |
350 | ||
351 | correctedEvents->SetBinContent(i, 1, events); | |
352 | } | |
353 | ||
354 | Printf("In |vtx-z| < 10 cm: %d events have been added", (Int_t) correctedEvents->Integral(vertexDist->FindBin(-9.9), vertexDist->FindBin(9.9), 1, 1)); | |
355 | ||
356 | //new TCanvas; correctedEvents->DrawCopy("TEXT"); | |
357 | //new TCanvas; kineBias->DrawCopy(); | |
358 | } | |
359 | ||
360 | fData->PrintInfo(ptCut); | |
361 | ||
362 | TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram(); | |
363 | ||
364 | // integrate multiplicity axis out (include under/overflow bins!!!) | |
365 | TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram(); | |
366 | ||
367 | TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e"); | |
368 | ||
369 | // create pt hist | |
370 | if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) | |
371 | { | |
372 | // reset all ranges | |
373 | dataHist->GetXaxis()->SetRange(0, 0); | |
374 | dataHist->GetYaxis()->SetRange(0, 0); | |
375 | dataHist->GetZaxis()->SetRange(0, 0); | |
376 | ||
377 | // vtx cut | |
378 | Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5); | |
379 | Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5); | |
380 | dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd); | |
381 | Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd); | |
382 | ||
383 | if (nEvents > 0) | |
384 | { | |
385 | // eta cut | |
386 | dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8)); | |
387 | Float_t etaWidth = 1.6; | |
388 | ||
389 | TH1D* ptHist = static_cast<TH1D*> (dataHist->Project3D("ze")); | |
390 | ||
391 | for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i) | |
392 | { | |
393 | Float_t binSize = fPtDist->GetBinWidth(i); | |
394 | fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth); | |
395 | fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth); | |
396 | } | |
397 | ||
398 | delete ptHist; | |
399 | } | |
400 | else | |
401 | printf("ERROR: nEvents is 0!\n"); | |
402 | } | |
403 | ||
404 | // reset all ranges | |
405 | dataHist->GetXaxis()->SetRange(0, 0); | |
406 | dataHist->GetYaxis()->SetRange(0, 0); | |
407 | dataHist->GetZaxis()->SetRange(0, 0); | |
408 | ||
409 | // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD) | |
410 | Int_t ptLowBin = 1; | |
411 | if (ptCut > 0 && (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)) | |
412 | ptLowBin = dataHist->GetZaxis()->FindBin(ptCut); | |
413 | ||
414 | //new TCanvas; dataHist->DrawCopy(); | |
415 | ||
416 | //dataHist->Sumw2(); | |
417 | dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1); | |
418 | printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1); | |
419 | TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e")); | |
420 | ||
421 | //new TCanvas; vtxVsEta->Draw("COLZ"); | |
422 | ||
423 | dataHist->GetZaxis()->SetRange(0, 0); | |
424 | vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle()); | |
425 | vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle()); | |
426 | ||
427 | if (vtxVsEta == 0) | |
428 | { | |
429 | printf("ERROR: pt/multiplicity integration failed\n"); | |
430 | return; | |
431 | } | |
432 | ||
433 | //new TCanvas(tag, tag, 800, 600); vtxVsEta->DrawCopy("COLZ"); | |
434 | ||
435 | // clear result histograms | |
436 | for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos) | |
437 | { | |
438 | fdNdEta[vertexPos]->Reset(); | |
439 | fdNdEtaPtCutOffCorrected[vertexPos]->Reset(); | |
440 | } | |
441 | ||
442 | //const Float_t vertexRangeBegin[kVertexBinning] = { -9.99, -9.99, 0.01 }; | |
443 | //const Float_t vertexRangeEnd[kVertexBinning] = { 9.99, -0.01, 9.99 }; | |
444 | const Float_t vertexRangeBegin[kVertexBinning] = { fvtxMin, fvtxMin, 0.01 }; | |
445 | const Float_t vertexRangeEnd[kVertexBinning] = { fvtxMax, -0.01, fvtxMax }; | |
446 | ||
447 | for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++) | |
448 | { | |
449 | // loop over vertex ranges | |
450 | for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos) | |
451 | { | |
452 | Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]); | |
453 | Int_t vertexBinEnd = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]); | |
454 | ||
455 | const Int_t *binBegin = 0; | |
456 | const Int_t maxBins = 20; | |
457 | const Int_t maxBins1 = 40; | |
458 | ||
459 | if (dataHist->GetNbinsY() != maxBins) | |
460 | AliFatal(Form("Binning of acceptance is different from data histogram: data=%d, acceptance=%d", dataHist->GetNbinsY(), maxBins)); | |
461 | ||
462 | // adjust acceptance range | |
463 | // produce with drawPlots.C: DetermineAcceptance(...) | |
464 | ||
465 | //const Int_t binBeginSPD[maxBins] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; | |
466 | //const Int_t binBeginTPC[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1}; | |
467 | //const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, -1, 14, 10, 8, 7, 6, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1, -1}; | |
468 | ||
469 | const Int_t binBeginSPD[maxBins1] = {19, 18, 17, 15, 14, 12, 10, 9, 8, 7, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; | |
470 | const Int_t binBeginTPC[maxBins1] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1}; | |
471 | const Int_t binBeginTPCITS[maxBins] = {-1, -1, -1, -1, -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, -1, -1, -1, -1, -1, -1}; | |
472 | if (fAnalysisMode & AliPWG0Helper::kSPD) | |
473 | { | |
474 | binBegin = binBeginSPD; | |
475 | } | |
476 | else if (fAnalysisMode & AliPWG0Helper::kTPC) | |
477 | { | |
478 | binBegin = binBeginTPC; | |
479 | } | |
480 | else if (fAnalysisMode & AliPWG0Helper::kTPCITS) | |
481 | { | |
482 | binBegin = binBeginTPCITS; | |
483 | } | |
484 | ||
485 | Int_t vtxBegin = 1; | |
486 | Int_t vtxEnd = maxBins; | |
487 | ||
488 | if (binBegin) | |
489 | { | |
490 | vtxBegin = binBegin[iEta - 1]; | |
491 | vtxEnd = (Int_t) vtxVsEta->GetNbinsX() + 1 - binBegin[maxBins - iEta]; | |
492 | } | |
493 | else | |
494 | Printf("WARNING: No acceptance applied!"); | |
495 | ||
496 | // eta range not accessible | |
497 | if (vtxBegin == -1) | |
498 | continue; | |
499 | ||
500 | //Printf("%d %d | %d %d", vtxBegin, vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)), vtxEnd, vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta))); | |
501 | //vtxBegin = vertexHist->GetXaxis()->FindBin(GetVtxMin(eta)); | |
502 | //vtxEnd = vertexHist->GetXaxis()->FindBin(-GetVtxMin(-eta)); | |
503 | ||
504 | //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta); | |
505 | //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd); | |
506 | vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin); | |
507 | vertexBinEnd = TMath::Min(vertexBinEnd, vtxEnd); | |
508 | //Printf(" after: %d %d", vertexBinBegin, vertexBinEnd); | |
509 | ||
510 | // no data for this bin | |
511 | if (vertexBinBegin > vertexBinEnd) | |
512 | { | |
513 | //Printf("Bin empty. Skipped"); | |
514 | continue; | |
515 | } | |
516 | ||
517 | Float_t totalEvents = 0; | |
518 | Float_t sum = 0; | |
519 | Float_t sumError2 = 0; | |
520 | Float_t unusedTracks = 0; | |
521 | Float_t unusedEvents = 0; | |
522 | for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++) | |
523 | { | |
524 | if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd) | |
525 | { | |
526 | if (vtxVsEta->GetBinContent(iVtx, iEta) != 0) | |
527 | { | |
528 | sum += vtxVsEta->GetBinContent(iVtx, iEta); | |
529 | ||
530 | if (sumError2 > 10e30) | |
531 | Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... "); | |
532 | ||
533 | sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2); | |
534 | } | |
535 | totalEvents += vertexHist->GetBinContent(iVtx); | |
536 | } | |
537 | else | |
538 | { | |
539 | unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta); | |
540 | unusedEvents += vertexHist->GetBinContent(iVtx); | |
541 | } | |
542 | } | |
543 | ||
544 | if (totalEvents == 0) | |
545 | { | |
546 | printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd); | |
547 | continue; | |
548 | } | |
549 | ||
550 | Float_t ptCutOffCorrection = 1; | |
551 | ||
552 | // find pt cut off correction factor | |
553 | if ((fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) && (fAnalysisMode & AliPWG0Helper::kFieldOn)) | |
554 | { | |
555 | if (correction && ptCut > 0) | |
556 | ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd); | |
557 | ||
558 | if (ptCutOffCorrection <= 0) | |
559 | { | |
560 | printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd); | |
561 | continue; | |
562 | } | |
563 | } | |
564 | ||
565 | //printf("Eta: %d (%f) Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f \n", iEta, vtxVsEta->GetYaxis()->GetBinCenter(iEta), vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection); | |
566 | ||
567 | Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta)); | |
568 | if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX()) | |
569 | { | |
570 | Float_t dndeta = sum / totalEvents; | |
571 | Float_t error = TMath::Sqrt(sumError2) / totalEvents; | |
572 | ||
573 | // correct for additional combinatorics | |
574 | Float_t combCorr = 0; | |
575 | if (combinatoricsCorrection) | |
576 | { | |
577 | combCorr = combinatoricsCorrection->GetBinContent(combinatoricsCorrection->GetXaxis()->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta))); | |
578 | dndeta *= combCorr; | |
579 | error *= combCorr; | |
580 | } | |
581 | ||
582 | dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin); | |
583 | error = error / fdNdEta[vertexPos]->GetBinWidth(bin); | |
584 | ||
585 | fdNdEta[vertexPos]->SetBinContent(bin, dndeta); | |
586 | fdNdEta[vertexPos]->SetBinError(bin, error); | |
587 | ||
588 | dndeta /= ptCutOffCorrection; | |
589 | error /= ptCutOffCorrection; | |
590 | ||
591 | fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta); | |
592 | fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error); | |
593 | ||
594 | //Printf("Bin %d has dN/deta = %f +- %f; %.2f tracks %.2f events (outside acceptance: %.2f tracks, %.2f events) (comb. corr: %f)", bin, dndeta, error, sum, totalEvents, unusedTracks, unusedEvents, combCorr); | |
595 | } | |
596 | } | |
597 | } | |
598 | } | |
599 | ||
600 | //____________________________________________________________________ | |
601 | Float_t dNdEtaAnalysis::GetVtxMin(Float_t eta) | |
602 | { | |
603 | // helper function for the SPD acceptance | |
604 | // the function returns the beginning of the acceptance window in z vertex position as function of eta | |
605 | // to get the maximum: -GetVtxMin(-eta) | |
606 | ||
607 | Float_t a[2] = { -15, 0 }; | |
608 | Float_t b[2] = { 0, -1.4 }; | |
609 | Float_t c[2] = { 15, -2.2 }; | |
610 | ||
611 | Float_t meanAB[2]; | |
612 | meanAB[0] = (b[0] + a[0]) / 2; | |
613 | meanAB[1] = (b[1] + a[1]) / 2; | |
614 | ||
615 | Float_t meanBC[2]; | |
616 | meanBC[0] = (c[0] + b[0]) / 2; | |
617 | meanBC[1] = (c[1] + b[1]) / 2; | |
618 | ||
619 | Float_t mAB = (b[1] - a[1]) / (b[0] - a[0]); | |
620 | Float_t mBC = (c[1] - b[1]) / (c[0] - b[0]); | |
621 | ||
622 | Float_t bAB = meanAB[1] - mAB * meanAB[0]; | |
623 | Float_t bBC = meanBC[1] - mBC * meanBC[0]; | |
624 | ||
625 | if (eta > b[1]) | |
626 | return 1.0 / mAB * eta - bAB / mAB; | |
627 | ||
628 | return 1.0 / mBC * eta - bBC / mBC; | |
629 | } | |
630 | ||
631 | //____________________________________________________________________ | |
632 | void dNdEtaAnalysis::SaveHistograms() | |
633 | { | |
634 | // save the histograms to a directory with the name of this class (retrieved from TNamed) | |
635 | ||
636 | gDirectory->mkdir(GetName()); | |
637 | gDirectory->cd(GetName()); | |
638 | ||
639 | if (fData) | |
640 | { | |
641 | fData->SaveHistograms(); | |
642 | } | |
643 | else | |
644 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n"); | |
645 | ||
646 | if (fMult) | |
647 | { | |
648 | fMult->Write(); | |
649 | } | |
650 | else | |
651 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n"); | |
652 | ||
653 | if (fPtDist) | |
654 | fPtDist ->Write(); | |
655 | else | |
656 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n"); | |
657 | ||
658 | for (Int_t i=0; i<kVertexBinning; ++i) | |
659 | { | |
660 | if (fdNdEta[i]) | |
661 | fdNdEta[i]->Write(); | |
662 | else | |
663 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i); | |
664 | ||
665 | if (fdNdEtaPtCutOffCorrected[i]) | |
666 | fdNdEtaPtCutOffCorrected[i]->Write(); | |
667 | else | |
668 | printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i); | |
669 | } | |
670 | ||
671 | TNamed named("fTag", fTag.Data()); | |
672 | named.Write(); | |
673 | ||
674 | TParameter<Int_t> param("fAnalysisMode", fAnalysisMode); | |
675 | param.Write(); | |
676 | ||
677 | gDirectory->cd("../"); | |
678 | } | |
679 | ||
680 | void dNdEtaAnalysis::LoadHistograms(const Char_t* dir) | |
681 | { | |
682 | // loads the histograms from a directory with the name of this class (retrieved from TNamed) | |
683 | ||
684 | if (!dir) | |
685 | dir = GetName(); | |
686 | ||
687 | gDirectory->cd(dir); | |
688 | ||
689 | fData->LoadHistograms(); | |
690 | fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName())); | |
691 | ||
692 | for (Int_t i=0; i<kVertexBinning; ++i) | |
693 | { | |
694 | fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName())); | |
695 | fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName())); | |
696 | } | |
697 | ||
698 | fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName())); | |
699 | ||
700 | if (static_cast<TNamed*> (gDirectory->Get("fTag"))) | |
701 | fTag = (static_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle(); | |
702 | ||
703 | if (static_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode"))) | |
704 | fAnalysisMode = (AliPWG0Helper::AnalysisMode) (static_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal(); | |
705 | ||
706 | gDirectory->cd("../"); | |
707 | } | |
708 | ||
709 | //____________________________________________________________________ | |
710 | void dNdEtaAnalysis::DrawHistograms(Bool_t simple) | |
711 | { | |
712 | // draws the histograms | |
713 | ||
714 | if (!simple) | |
715 | { | |
716 | if (fData) | |
717 | fData->DrawHistograms(GetName()); | |
718 | ||
719 | TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400); | |
720 | canvas->Divide(2, 1); | |
721 | ||
722 | canvas->cd(1); | |
723 | if (fdNdEtaPtCutOffCorrected[0]) | |
724 | fdNdEtaPtCutOffCorrected[0]->DrawCopy(); | |
725 | ||
726 | if (fdNdEta[0]) | |
727 | { | |
728 | fdNdEta[0]->SetLineColor(kRed); | |
729 | fdNdEta[0]->DrawCopy("SAME"); | |
730 | } | |
731 | ||
732 | canvas->cd(2); | |
733 | if (fPtDist) | |
734 | fPtDist->DrawCopy(); | |
735 | } | |
736 | ||
737 | // histograms for different vertices? | |
738 | if (kVertexBinning > 0) | |
739 | { | |
740 | // doesnt work, but i dont get it, giving up... | |
741 | TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450); | |
742 | TCanvas* canvas3 = 0; | |
743 | if (!simple) | |
744 | canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450); | |
745 | ||
746 | //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2); | |
747 | //printf("%d\n", yPads); | |
748 | //canvas2->Divide(2, yPads); | |
749 | ||
750 | TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4); | |
751 | ||
752 | for (Int_t i=0; i<kVertexBinning; ++i) | |
753 | { | |
754 | if (fdNdEtaPtCutOffCorrected[i]) | |
755 | { | |
756 | canvas2->cd(); | |
757 | ||
758 | fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1); | |
759 | fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME"); | |
760 | legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1)); | |
761 | } | |
762 | if (canvas3 && fdNdEta[i]) | |
763 | { | |
764 | canvas3->cd(); | |
765 | ||
766 | fdNdEta[i]->SetLineColor(i+1); | |
767 | fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME"); | |
768 | } | |
769 | } | |
770 | ||
771 | canvas2->cd(); | |
772 | legend->Draw(); | |
773 | canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName())); | |
774 | ||
775 | if (canvas3) | |
776 | { | |
777 | canvas3->cd(); | |
778 | legend->Draw(); | |
779 | } | |
780 | } | |
781 | ||
782 | if (kVertexBinning == 3) | |
783 | { | |
784 | TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone")); | |
785 | TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2")); | |
786 | ||
787 | if (clone && clone2) | |
788 | { | |
789 | TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450); | |
790 | ||
791 | clone->Divide(fdNdEtaPtCutOffCorrected[0]); | |
792 | clone->GetYaxis()->SetRangeUser(0.95, 1.05); | |
793 | clone->DrawCopy(); | |
794 | ||
795 | clone2->Divide(fdNdEtaPtCutOffCorrected[0]); | |
796 | clone2->DrawCopy("SAME"); | |
797 | ||
798 | TLine* line = new TLine(-1, 1, 1, 1); | |
799 | line->Draw(); | |
800 | ||
801 | canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName())); | |
802 | } | |
803 | } | |
804 | } | |
805 | ||
806 | Long64_t dNdEtaAnalysis::Merge(TCollection* list) | |
807 | { | |
808 | // Merges a list of dNdEtaAnalysis objects with this one. | |
809 | // This is needed for PROOF. | |
810 | // Returns the number of merged objects (including this) | |
811 | ||
812 | if (!list) | |
813 | return 0; | |
814 | ||
815 | if (list->IsEmpty()) | |
816 | return 1; | |
817 | ||
818 | TIterator* iter = list->MakeIterator(); | |
819 | TObject* obj; | |
820 | ||
821 | // sub collections | |
822 | const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning | |
823 | TList* collections[nCollections]; | |
824 | for (Int_t i=0; i<nCollections; ++i) | |
825 | collections[i] = new TList; | |
826 | ||
827 | Int_t count = 0; | |
828 | while ((obj = iter->Next())) | |
829 | { | |
830 | dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj); | |
831 | if (entry == 0) | |
832 | continue; | |
833 | ||
834 | collections[0]->Add(entry->fData); | |
835 | collections[1]->Add(entry->fMult); | |
836 | collections[2]->Add(entry->fPtDist); | |
837 | ||
838 | for (Int_t i=0; i<kVertexBinning; ++i) | |
839 | { | |
840 | collections[3+i]->Add(entry->fdNdEta[i]); | |
841 | collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]); | |
842 | } | |
843 | ||
844 | ++count; | |
845 | } | |
846 | ||
847 | fData->Merge(collections[0]); | |
848 | fMult->Merge(collections[1]); | |
849 | fPtDist->Merge(collections[2]); | |
850 | ||
851 | for (Int_t i=0; i<kVertexBinning; ++i) | |
852 | { | |
853 | fdNdEta[i]->Merge(collections[3+i]); | |
854 | fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]); | |
855 | } | |
856 | ||
857 | for (Int_t i=0; i<nCollections; ++i) | |
858 | delete collections[i]; | |
859 | ||
860 | return count+1; | |
861 | } | |
862 | ||
863 |