Update
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
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 (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 = dynamic_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       if (fAnalysisMode & AliPWG0Helper::kSPD)
470       {
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};
472         binBegin = binBeginSPD;
473       }
474       else if (fAnalysisMode & AliPWG0Helper::kTPC)
475       {
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};
477         binBegin = binBeginTPC;
478       }
479       else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
480       {
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};
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 (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))
701     fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
702
703   if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
704     fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_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