Update
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
CommitLineData
dc740de4 1/* $Id$ */
2
75ec0f41 3#include "dNdEtaAnalysis.h"
4
ceb5d1b5 5#include <TFile.h>
45e97e28 6#include <TH3F.h>
74fd10b3 7#include <TH2F.h>
8#include <TH1F.h>
ceb5d1b5 9#include <TMath.h>
10#include <TCanvas.h>
7029240a 11#include <TCollection.h>
12#include <TIterator.h>
13#include <TList.h>
fcf2fb36 14#include <TLegend.h>
74fd10b3 15#include <TLine.h>
745d6088 16#include <TParameter.h>
fcf2fb36 17
1d107532 18#include <AliLog.h>
45e97e28 19#include "AlidNdEtaCorrection.h"
74fd10b3 20#include <AliCorrection.h>
21#include <AliPWG0Helper.h>
22#include <AliCorrectionMatrix2D.h>
23#include <AliCorrectionMatrix3D.h>
ceb5d1b5 24
75ec0f41 25//____________________________________________________________________
b7f4a1fd 26ClassImp(dNdEtaAnalysis)
75ec0f41 27
28//____________________________________________________________________
1afae8ff 29dNdEtaAnalysis::dNdEtaAnalysis() :
30 TNamed(),
31 fData(0),
7307d52c 32 fMult(0),
745d6088 33 fPtDist(0),
34 fAnalysisMode(AliPWG0Helper::kInvalid),
68fa248f 35 fTag(),
36 fvtxMin(-9.99),
37 fvtxMax(9.99)
1afae8ff 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//____________________________________________________________________
a6e0ebfe 49dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
16e24ca3 50 TNamed(name, title),
45e97e28 51 fData(0),
7307d52c 52 fMult(0),
745d6088 53 fPtDist(0),
54 fAnalysisMode(analysisMode),
68fa248f 55 fTag(),
56 fvtxMin(-9.99),
57 fvtxMax(9.99)
7029240a 58{
6bf0714d 59 // constructor
4dd2ad81 60
770a1f1d 61 fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
45e97e28 62
74fd10b3 63 // do not add this hists to the directory
64 Bool_t oldStatus = TH1::AddDirectoryStatus();
65 TH1::AddDirectory(kFALSE);
1afae8ff 66
7307d52c 67 fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5);
68
567160d6 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());
74fd10b3 71
72 fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
1afae8ff 73
45e97e28 74 for (Int_t i=1; i<kVertexBinning; ++i)
7029240a 75 {
74fd10b3 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)));
7029240a 78 }
75ec0f41 79
567160d6 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());
74fd10b3 81
82 TH1::AddDirectory(oldStatus);
75ec0f41 83}
84
85//____________________________________________________________________
16e24ca3 86dNdEtaAnalysis::~dNdEtaAnalysis()
87{
88 // destructor
89
0ab29cfa 90 if (fData)
91 {
92 delete fData;
93 fData = 0;
94 }
16e24ca3 95
7307d52c 96 if (fMult)
97 {
98 delete fMult;
99 fMult = 0;
100 }
101
16e24ca3 102 for (Int_t i=0; i<kVertexBinning; ++i)
103 {
0ab29cfa 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 }
16e24ca3 114 }
1afae8ff 115
0ab29cfa 116 if (fPtDist)
117 {
118 delete fPtDist;
119 fPtDist = 0;
120 }
16e24ca3 121}
122
123//_____________________________________________________________________________
124dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
125 TNamed(c),
45e97e28 126 fData(0),
c61a7285 127 fMult(0),
745d6088 128 fPtDist(0),
129 fAnalysisMode(AliPWG0Helper::kInvalid),
68fa248f 130 fTag(),
131 fvtxMin(-9.99),
132 fvtxMax(9.99)
16e24ca3 133{
134 //
135 // dNdEtaAnalysis copy constructor
136 //
137
138 ((dNdEtaAnalysis &) c).Copy(*this);
139}
140
141//_____________________________________________________________________________
142dNdEtaAnalysis &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//_____________________________________________________________________________
153void dNdEtaAnalysis::Copy(TObject &c) const
154{
155 //
156 // Copy function
157 //
158
159 dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
160
74fd10b3 161 target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
7307d52c 162 target.fMult = dynamic_cast<TH1F*> (fMult->Clone());
16e24ca3 163
164 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 165 {
74fd10b3 166 target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
167 target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
1afae8ff 168 }
16e24ca3 169
74fd10b3 170 target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
45e97e28 171
745d6088 172 target.fAnalysisMode = fAnalysisMode;
173 target.fTag = fTag;
174
16e24ca3 175 TNamed::Copy((TNamed &) c);
176}
177
178//____________________________________________________________________
74fd10b3 179void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
6bf0714d 180{
181 // fills a track into the histograms
182
74fd10b3 183 fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
75ec0f41 184}
185
186//____________________________________________________________________
74fd10b3 187void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
6bf0714d 188{
189 // fills an event into the histograms
190
74fd10b3 191 fData->GetEventCorrection()->FillMeas(vtx, n);
75ec0f41 192}
193
194//____________________________________________________________________
7307d52c 195void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
196{
197 // fills a triggered event into the histograms
198
199 fMult->Fill(n);
200}
201
202//____________________________________________________________________
81be4ee8 203void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag, Int_t backgroundSubtraction, TH1* combinatoricsCorrection)
fcf2fb36 204{
74fd10b3 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
3dfa46a4 208 // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
74fd10b3 209 //
210
51f6de65 211 fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
745d6088 212 Printf("\n\n%s", fTag.Data());
81be4ee8 213
214 if (combinatoricsCorrection)
215 Printf("Combinatorics subtraction active!");
3dfa46a4 216
74fd10b3 217 // set corrections to 1
218 fData->SetCorrectionToUnity();
219
220 if (correction && correctionType != AlidNdEtaCorrection::kNone)
221 {
5a6310fe 222 TH3* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
223 TH2* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
74fd10b3 224
225 if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
226 trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
6bf0714d 227
dd367a14 228 if (correctionType >= AlidNdEtaCorrection::kVertexReco)
74fd10b3 229 {
230 trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
231 eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
7307d52c 232
0fc41645 233 // set bin with multiplicity 0 to unity (correction has no meaning in this bin)
7307d52c 234 for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++)
235 eventCorr->SetBinContent(i, 1, 1);
dd367a14 236 }
74fd10b3 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 }
81be4ee8 258 case AlidNdEtaCorrection::kOnePart :
259 {
260 trackCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetTrackCorrection()->GetCorrectionHistogram());
261 eventCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->GetCorrectionHistogram());
262 break;
263 }
74fd10b3 264 default : break;
265 }
266 }
267 else
fcf2fb36 268 printf("INFO: No correction applied\n");
269
567160d6 270 TH2F* rawMeasured = (TH2F*) fData->GetEventCorrection()->GetMeasuredHistogram()->Clone("rawMeasured");
271
12bb57f1 272 fData->ResetErrorsOnCorrections();
74fd10b3 273 fData->Multiply();
81be4ee8 274
7307d52c 275 if (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
a7f69e56 282 // at the same time calculate expectation from MC (not used, just to check the value)
7307d52c 283
567160d6 284 //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
7307d52c 285 TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
286
81be4ee8 287 TH2* eAll = correction->GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram();
567160d6 288 TH2* eTrig = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
289 TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
81be4ee8 290 TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsY()+1);
567160d6 291
7307d52c 292 //new TCanvas; correctedEvents->DrawCopy("TEXT");
293
294 // start above 0 mult. bin with integration
567160d6 295 TH1* vertexDist = rawMeasured->ProjectionX("vertexdist_measured", 2, rawMeasured->GetNbinsY()+1);
7307d52c 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);
81be4ee8 302
303 if (backgroundSubtraction > 0)
304 {
305 triggeredEventsWith0Mult -= backgroundSubtraction;
306 Printf("Subtracted %d background events from 0 mult. bin", backgroundSubtraction);
307 }
7307d52c 308
51f6de65 309 TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
310 kineBias->Reset();
81be4ee8 311
a7f69e56 312 // loop over vertex bins
567160d6 313 for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
7307d52c 314 {
745d6088 315 Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
7307d52c 316 Double_t events = alpha * triggeredEventsWith0Mult;
317
567160d6 318 if (eTrigVtx_projx->GetBinContent(i) == 0)
319 continue;
320
a7f69e56 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
81be4ee8 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);
a7f69e56 342
81be4ee8 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 }
7307d52c 350
351 correctedEvents->SetBinContent(i, 1, events);
352 }
a7f69e56 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));
7307d52c 355
356 //new TCanvas; correctedEvents->DrawCopy("TEXT");
69b09e3b 357 //new TCanvas; kineBias->DrawCopy();
7307d52c 358 }
359
3dfa46a4 360 fData->PrintInfo(ptCut);
74fd10b3 361
5a6310fe 362 TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
74fd10b3 363
364 // integrate multiplicity axis out (include under/overflow bins!!!)
5a6310fe 365 TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
3dfa46a4 366
7307d52c 367 TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
fcf2fb36 368
1afae8ff 369 // create pt hist
a7f69e56 370 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
1afae8ff 371 {
372 // reset all ranges
74fd10b3 373 dataHist->GetXaxis()->SetRange(0, 0);
374 dataHist->GetYaxis()->SetRange(0, 0);
375 dataHist->GetZaxis()->SetRange(0, 0);
1afae8ff 376
377 // vtx cut
74fd10b3 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);
1afae8ff 382
74fd10b3 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;
1afae8ff 388
74fd10b3 389 TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
1afae8ff 390
74fd10b3 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 }
1afae8ff 397
74fd10b3 398 delete ptHist;
399 }
400 else
401 printf("ERROR: nEvents is 0!\n");
1afae8ff 402 }
403
404 // reset all ranges
74fd10b3 405 dataHist->GetXaxis()->SetRange(0, 0);
406 dataHist->GetYaxis()->SetRange(0, 0);
407 dataHist->GetZaxis()->SetRange(0, 0);
1afae8ff 408
51f6de65 409 // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
1afae8ff 410 Int_t ptLowBin = 1;
a7f69e56 411 if (ptCut > 0 && (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS))
74fd10b3 412 ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
69b09e3b 413
414 //new TCanvas; dataHist->DrawCopy();
1afae8ff 415
69b09e3b 416 //dataHist->Sumw2();
7307d52c 417 dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
51f6de65 418 printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
74fd10b3 419 TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
69b09e3b 420
421 //new TCanvas; vtxVsEta->Draw("COLZ");
74fd10b3 422
423 dataHist->GetZaxis()->SetRange(0, 0);
424 vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
425 vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
1afae8ff 426
847489f7 427 if (vtxVsEta == 0)
428 {
51f6de65 429 printf("ERROR: pt/multiplicity integration failed\n");
847489f7 430 return;
fcf2fb36 431 }
432
51f6de65 433 //new TCanvas(tag, tag, 800, 600); vtxVsEta->DrawCopy("COLZ");
81be4ee8 434
567160d6 435 // clear result histograms
436 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
437 {
438 fdNdEta[vertexPos]->Reset();
439 fdNdEtaPtCutOffCorrected[vertexPos]->Reset();
440 }
144ff489 441
68fa248f 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 };
b4b9cacc 446
74fd10b3 447 for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
5af55649 448 {
745d6088 449 // loop over vertex ranges
7029240a 450 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
451 {
745d6088 452 Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]);
453 Int_t vertexBinEnd = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
fcf2fb36 454
ea441adf 455 const Int_t *binBegin = 0;
68fa248f 456 const Int_t maxBins = 20;
457 const Int_t maxBins1 = 40;
1d107532 458
459 if (dataHist->GetNbinsY() != maxBins)
460 AliFatal(Form("Binning of acceptance is different from data histogram: data=%d, acceptance=%d", dataHist->GetNbinsY(), maxBins));
1c15d51a 461
462 // adjust acceptance range
69b09e3b 463 // produce with drawPlots.C: DetermineAcceptance(...)
68fa248f 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};
db7efefe 468
a7f69e56 469 if (fAnalysisMode & AliPWG0Helper::kSPD)
fcf2fb36 470 {
68fa248f 471 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};
1c15d51a 472 binBegin = binBeginSPD;
fcf2fb36 473 }
a7f69e56 474 else if (fAnalysisMode & AliPWG0Helper::kTPC)
ea441adf 475 {
68fa248f 476 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};
81be4ee8 477 binBegin = binBeginTPC;
ea441adf 478 }
a7f69e56 479 else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
ea441adf 480 {
68fa248f 481 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};
81be4ee8 482 binBegin = binBeginTPCITS;
ea441adf 483 }
484
51f6de65 485 Int_t vtxBegin = 1;
69b09e3b 486 Int_t vtxEnd = maxBins;
51f6de65 487
488 if (binBegin)
489 {
490 vtxBegin = binBegin[iEta - 1];
68fa248f 491 vtxEnd = (Int_t) vtxVsEta->GetNbinsX() + 1 - binBegin[maxBins - iEta];
51f6de65 492 }
493 else
494 Printf("WARNING: No acceptance applied!");
ea441adf 495
496 // eta range not accessible
497 if (vtxBegin == -1)
498 continue;
499
ea441adf 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));
1c15d51a 503
b9dbe907 504 //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
51f6de65 505 //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
ea441adf 506 vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
507 vertexBinEnd = TMath::Min(vertexBinEnd, vtxEnd);
51f6de65 508 //Printf(" after: %d %d", vertexBinBegin, vertexBinEnd);
7029240a 509
745d6088 510 // no data for this bin
511 if (vertexBinBegin > vertexBinEnd)
567160d6 512 {
513 //Printf("Bin empty. Skipped");
745d6088 514 continue;
567160d6 515 }
b4b9cacc 516
1c15d51a 517 Float_t totalEvents = 0;
5af55649 518 Float_t sum = 0;
519 Float_t sumError2 = 0;
1c15d51a 520 Float_t unusedTracks = 0;
521 Float_t unusedEvents = 0;
522 for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++)
567160d6 523 {
1c15d51a 524 if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd)
5af55649 525 {
1c15d51a 526 if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
527 {
528 sum += vtxVsEta->GetBinContent(iVtx, iEta);
2e88424e 529
1c15d51a 530 if (sumError2 > 10e30)
531 Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
2e88424e 532
1c15d51a 533 sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
534 }
535 totalEvents += vertexHist->GetBinContent(iVtx);
5af55649 536 }
1c15d51a 537 else
538 {
539 unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta);
69b09e3b 540 unusedEvents += vertexHist->GetBinContent(iVtx);
1c15d51a 541 }
542 }
543
544 if (totalEvents == 0)
545 {
546 printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
547 continue;
5af55649 548 }
7029240a 549
1afae8ff 550 Float_t ptCutOffCorrection = 1;
1afae8ff 551
51f6de65 552 // find pt cut off correction factor
a7f69e56 553 if ((fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) && (fAnalysisMode & AliPWG0Helper::kFieldOn))
847489f7 554 {
51f6de65 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 }
847489f7 563 }
564
69b09e3b 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);
9e952c39 566
144ff489 567 Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
745d6088 568 if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX())
144ff489 569 {
7307d52c 570 Float_t dndeta = sum / totalEvents;
571 Float_t error = TMath::Sqrt(sumError2) / totalEvents;
7029240a 572
81be4ee8 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
144ff489 582 dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
583 error = error / fdNdEta[vertexPos]->GetBinWidth(bin);
7029240a 584
144ff489 585 fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
586 fdNdEta[vertexPos]->SetBinError(bin, error);
1afae8ff 587
144ff489 588 dndeta /= ptCutOffCorrection;
589 error /= ptCutOffCorrection;
1afae8ff 590
144ff489 591 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
592 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
2e88424e 593
81be4ee8 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);
144ff489 595 }
75ec0f41 596 }
847489f7 597 }
75ec0f41 598}
599
75ec0f41 600//____________________________________________________________________
567160d6 601Float_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//____________________________________________________________________
6bf0714d 632void dNdEtaAnalysis::SaveHistograms()
633{
634 // save the histograms to a directory with the name of this class (retrieved from TNamed)
75ec0f41 635
7029240a 636 gDirectory->mkdir(GetName());
637 gDirectory->cd(GetName());
5fbd0b17 638
1afae8ff 639 if (fData)
640 {
74fd10b3 641 fData->SaveHistograms();
1afae8ff 642 }
643 else
644 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
645
7307d52c 646 if (fMult)
647 {
648 fMult->Write();
649 }
650 else
651 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n");
652
d09fb536 653 if (fPtDist)
654 fPtDist ->Write();
655 else
656 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
657
7029240a 658 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 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 }
75ec0f41 670
745d6088 671 TNamed named("fTag", fTag.Data());
672 named.Write();
673
674 TParameter<Int_t> param("fAnalysisMode", fAnalysisMode);
675 param.Write();
676
75ec0f41 677 gDirectory->cd("../");
678}
679
74fd10b3 680void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
5fbd0b17 681{
6bf0714d 682 // loads the histograms from a directory with the name of this class (retrieved from TNamed)
683
74fd10b3 684 if (!dir)
685 dir = GetName();
5fbd0b17 686
74fd10b3 687 gDirectory->cd(dir);
5fbd0b17 688
74fd10b3 689 fData->LoadHistograms();
7307d52c 690 fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName()));
5fbd0b17 691
692 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 693 {
74fd10b3 694 fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
695 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
1afae8ff 696 }
5fbd0b17 697
74fd10b3 698 fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
0ab29cfa 699
1c15d51a 700 if (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))
701 fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
745d6088 702
1c15d51a 703 if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
704 fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal();
745d6088 705
5fbd0b17 706 gDirectory->cd("../");
707}
708
ceb5d1b5 709//____________________________________________________________________
74fd10b3 710void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
ceb5d1b5 711{
6bf0714d 712 // draws the histograms
7307d52c 713
74fd10b3 714 if (!simple)
715 {
716 if (fData)
717 fData->DrawHistograms(GetName());
6bf0714d 718
74fd10b3 719 TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
720 canvas->Divide(2, 1);
ceb5d1b5 721
74fd10b3 722 canvas->cd(1);
723 if (fdNdEtaPtCutOffCorrected[0])
144ff489 724 fdNdEtaPtCutOffCorrected[0]->DrawCopy();
ceb5d1b5 725
74fd10b3 726 if (fdNdEta[0])
727 {
728 fdNdEta[0]->SetLineColor(kRed);
144ff489 729 fdNdEta[0]->DrawCopy("SAME");
74fd10b3 730 }
1afae8ff 731
74fd10b3 732 canvas->cd(2);
733 if (fPtDist)
144ff489 734 fPtDist->DrawCopy();
1afae8ff 735 }
736
fcf2fb36 737 // histograms for different vertices?
738 if (kVertexBinning > 0)
739 {
740 // doesnt work, but i dont get it, giving up...
74fd10b3 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
fcf2fb36 746 //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
747 //printf("%d\n", yPads);
748 //canvas2->Divide(2, yPads);
749
0448e811 750 TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
fcf2fb36 751
5af55649 752 for (Int_t i=0; i<kVertexBinning; ++i)
fcf2fb36 753 {
1afae8ff 754 if (fdNdEtaPtCutOffCorrected[i])
fcf2fb36 755 {
74fd10b3 756 canvas2->cd();
757
1afae8ff 758 fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
144ff489 759 fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
1afae8ff 760 legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
fcf2fb36 761 }
74fd10b3 762 if (canvas3 && fdNdEta[i])
763 {
764 canvas3->cd();
765
766 fdNdEta[i]->SetLineColor(i+1);
144ff489 767 fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
74fd10b3 768 }
fcf2fb36 769 }
770
74fd10b3 771 canvas2->cd();
fcf2fb36 772 legend->Draw();
74fd10b3 773 canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
774
775 if (canvas3)
776 {
777 canvas3->cd();
778 legend->Draw();
779 }
fcf2fb36 780 }
7307d52c 781
74fd10b3 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"));
7307d52c 786
74fd10b3 787 if (clone && clone2)
788 {
789 TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
7307d52c 790
74fd10b3 791 clone->Divide(fdNdEtaPtCutOffCorrected[0]);
792 clone->GetYaxis()->SetRangeUser(0.95, 1.05);
144ff489 793 clone->DrawCopy();
7307d52c 794
74fd10b3 795 clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
144ff489 796 clone2->DrawCopy("SAME");
74fd10b3 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 }
7029240a 804}
805
806Long64_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
7307d52c 822 const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning
7029240a 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
45e97e28 834 collections[0]->Add(entry->fData);
0fc41645 835 collections[1]->Add(entry->fMult);
836 collections[2]->Add(entry->fPtDist);
7029240a 837
838 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 839 {
7307d52c 840 collections[3+i]->Add(entry->fdNdEta[i]);
841 collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
1afae8ff 842 }
7029240a 843
844 ++count;
845 }
846
45e97e28 847 fData->Merge(collections[0]);
7307d52c 848 fMult->Merge(collections[1]);
849 fPtDist->Merge(collections[2]);
850
7029240a 851 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 852 {
7307d52c 853 fdNdEta[i]->Merge(collections[3+i]);
854 fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]);
1afae8ff 855 }
7029240a 856
857 for (Int_t i=0; i<nCollections; ++i)
858 delete collections[i];
859
860 return count+1;
ceb5d1b5 861}
567160d6 862
863