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