updated dndeta analysis
[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 {
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
46 //____________________________________________________________________
47 dNdEtaAnalysis::dNdEtaAnalysis(const Char_t* name, const Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
48   TNamed(name, title),
49   fData(0),
50   fMult(0),
51   fPtDist(0),
52   fAnalysisMode(analysisMode),
53   fTag()
54 {
55   // constructor
56
57   fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
58
59   // do not add this hists to the directory
60   Bool_t oldStatus = TH1::AddDirectoryStatus();
61   TH1::AddDirectory(kFALSE);
62
63   fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5);
64
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());
67
68   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
69
70   for (Int_t i=1; i<kVertexBinning; ++i)
71   {
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)));
74   }
75
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());
77
78   TH1::AddDirectory(oldStatus);
79 }
80
81 //____________________________________________________________________
82 dNdEtaAnalysis::~dNdEtaAnalysis()
83 {
84   // destructor
85
86   if (fData)
87   {
88     delete fData;
89     fData = 0;
90   }
91
92   if (fMult)
93   {
94     delete fMult;
95     fMult = 0;
96   }
97
98   for (Int_t i=0; i<kVertexBinning; ++i)
99   {
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     }
110   }
111
112   if (fPtDist)
113   {
114     delete fPtDist;
115     fPtDist = 0;
116   }
117 }
118
119 //_____________________________________________________________________________
120 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
121   TNamed(c),
122   fData(0),
123   fMult(0),
124   fPtDist(0),
125   fAnalysisMode(AliPWG0Helper::kInvalid),
126   fTag()
127 {
128   //
129   // dNdEtaAnalysis copy constructor
130   //
131
132   ((dNdEtaAnalysis &) c).Copy(*this);
133 }
134
135 //_____________________________________________________________________________
136 dNdEtaAnalysis &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 //_____________________________________________________________________________
147 void dNdEtaAnalysis::Copy(TObject &c) const
148 {
149   //
150   // Copy function
151   //
152
153   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
154
155   target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
156   target.fMult = dynamic_cast<TH1F*> (fMult->Clone());
157
158   for (Int_t i=0; i<kVertexBinning; ++i)
159   {
160     target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
161     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
162   }
163
164   target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
165
166   target.fAnalysisMode = fAnalysisMode;
167   target.fTag = fTag;
168
169   TNamed::Copy((TNamed &) c);
170 }
171
172 //____________________________________________________________________
173 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
174 {
175   // fills a track into the histograms
176
177   fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
178 }
179
180 //____________________________________________________________________
181 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
182 {
183   // fills an event into the histograms
184
185   fData->GetEventCorrection()->FillMeas(vtx, n);
186 }
187
188 //____________________________________________________________________
189 void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
190 {
191   // fills a triggered event into the histograms
192
193   fMult->Fill(n);
194 }
195
196 //____________________________________________________________________
197 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag, Int_t backgroundSubtraction, TH1* combinatoricsCorrection)
198 {
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
202   // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
203   //
204
205   fTag.Form("Correcting dN/deta spectrum (data: %d) >>> %s <<<. Correction type: %d, pt cut: %.2f.", (Int_t) fAnalysisMode, tag, (Int_t) correctionType, ptCut);
206   Printf("\n\n%s", fTag.Data());
207   
208   if (combinatoricsCorrection)
209     Printf("Combinatorics subtraction active!");
210
211   // set corrections to 1
212   fData->SetCorrectionToUnity();
213
214   if (correction && correctionType != AlidNdEtaCorrection::kNone)
215   {
216     TH3* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
217     TH2* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
218
219     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
220       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
221
222     if (correctionType >= AlidNdEtaCorrection::kVertexReco)
223     {
224       trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
225       eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
226
227       // set bin with multiplicity 0 to unity (correction has no meaning in this bin)
228       for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++)
229         eventCorr->SetBinContent(i, 1, 1);
230     }
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       }
252       case AlidNdEtaCorrection::kOnePart :
253       {
254         trackCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetTrackCorrection()->GetCorrectionHistogram());
255         eventCorr->Multiply(correction->GetTriggerBiasCorrectionOnePart()->GetEventCorrection()->GetCorrectionHistogram());
256         break;
257       }
258       default : break;
259     }
260   }
261   else
262     printf("INFO: No correction applied\n");
263
264   TH2F* rawMeasured = (TH2F*) fData->GetEventCorrection()->GetMeasuredHistogram()->Clone("rawMeasured");
265
266   fData->ResetErrorsOnCorrections();
267   fData->Multiply();
268   
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
276     // at the same time calculate expectation from MC (not used, just to check the value)
277
278     //TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
279     TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
280
281     TH2* eAll  =    correction->GetCorrection(correctionType)->GetEventCorrection()->GetGeneratedHistogram();
282     TH2* eTrig =    correction->GetVertexRecoCorrection()->GetEventCorrection()->GetGeneratedHistogram();
283     TH2* eTrigVtx = correction->GetVertexRecoCorrection()->GetEventCorrection()->GetMeasuredHistogram();
284     TH1* eTrigVtx_projx = eTrigVtx->ProjectionX("eTrigVtx_projx", 2, eTrigVtx->GetNbinsY()+1);
285
286     //new TCanvas; correctedEvents->DrawCopy("TEXT");
287
288     // start above 0 mult. bin with integration
289     TH1* vertexDist = rawMeasured->ProjectionX("vertexdist_measured", 2, rawMeasured->GetNbinsY()+1);
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);
296     
297     if (backgroundSubtraction > 0)
298     {
299       triggeredEventsWith0Mult -= backgroundSubtraction;
300       Printf("Subtracted %d background events from 0 mult. bin", backgroundSubtraction);
301     }
302
303     TH1* kineBias = (TH1*) vertexDist->Clone("kineBias");
304     kineBias->Reset();
305     
306     // loop over vertex bins
307     for (Int_t i = 1; i <= rawMeasured->GetNbinsX(); i++)
308     {
309       Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
310       Double_t events = alpha * triggeredEventsWith0Mult;
311
312       if (eTrigVtx_projx->GetBinContent(i) == 0)
313         continue;
314
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
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);
336
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       }
344
345       correctedEvents->SetBinContent(i, 1, events);
346     }
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));
349
350     //new TCanvas; correctedEvents->DrawCopy("TEXT");
351     //new TCanvas; kineBias->DrawCopy();
352   }
353
354   fData->PrintInfo(ptCut);
355
356   TH3* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
357
358   // integrate multiplicity axis out (include under/overflow bins!!!)
359   TH2* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
360
361   TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
362
363   // create pt hist
364   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
365   {
366     // reset all ranges
367     dataHist->GetXaxis()->SetRange(0, 0);
368     dataHist->GetYaxis()->SetRange(0, 0);
369     dataHist->GetZaxis()->SetRange(0, 0);
370
371     // vtx cut
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);
376
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;
382
383       TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
384
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       }
391
392       delete ptHist;
393     }
394     else
395       printf("ERROR: nEvents is 0!\n");
396   }
397
398   // reset all ranges
399   dataHist->GetXaxis()->SetRange(0, 0);
400   dataHist->GetYaxis()->SetRange(0, 0);
401   dataHist->GetZaxis()->SetRange(0, 0);
402
403   // integrate over pt (with pt cut) (TPC, TPCITS) or multiplicity (SPD)
404   Int_t ptLowBin = 1;
405   if (ptCut > 0 && (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS))
406     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
407     
408   //new TCanvas; dataHist->DrawCopy();
409
410   //dataHist->Sumw2();
411   dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
412   printf("pt/multiplicity range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
413   TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
414   
415   //new TCanvas; vtxVsEta->Draw("COLZ");
416
417   dataHist->GetZaxis()->SetRange(0, 0);
418   vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
419   vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
420
421   if (vtxVsEta == 0)
422   {
423     printf("ERROR: pt/multiplicity integration failed\n");
424     return;
425   }
426
427   //new TCanvas(tag, tag, 800, 600);  vtxVsEta->DrawCopy("COLZ");
428   
429   // clear result histograms
430   for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
431   {
432     fdNdEta[vertexPos]->Reset();
433     fdNdEtaPtCutOffCorrected[vertexPos]->Reset();
434   }
435
436   const Float_t vertexRangeBegin[kVertexBinning] = { -9.99,  -9.99,  0.01 };
437   const Float_t vertexRangeEnd[kVertexBinning]   = {  9.99,  -0.01,  9.99 };
438
439   for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
440   {
441     // loop over vertex ranges
442     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
443     {
444       Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]);
445       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
446
447       const Int_t *binBegin = 0;
448       const Int_t maxBins = 40;
449       
450       if (dataHist->GetNbinsY() != maxBins)
451         AliFatal(Form("Binning of acceptance is different from data histogram: data=%d, acceptance=%d", dataHist->GetNbinsY(), maxBins));
452
453       // adjust acceptance range
454       // produce with drawPlots.C: DetermineAcceptance(...)
455       if (fAnalysisMode & AliPWG0Helper::kSPD)
456       {
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         
460         binBegin = binBeginSPD;
461       }
462       else if (fAnalysisMode & AliPWG0Helper::kTPC)
463       {
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};
465
466         binBegin = binBeginTPC;
467       }
468       else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
469       {
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;
473       }
474
475       Int_t vtxBegin = 1;
476       Int_t vtxEnd   = maxBins;
477
478       if (binBegin)
479       {
480         vtxBegin = binBegin[iEta - 1];
481         vtxEnd = 18 + 1 - binBegin[maxBins - iEta];
482       }
483       else
484         Printf("WARNING: No acceptance applied!");
485       
486       // eta range not accessible
487       if (vtxBegin == -1)
488         continue;
489       
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));
493
494       //Float_t eta = vtxVsEta->GetYaxis()->GetBinCenter(iEta);
495       //printf("Eta bin: %d (%f) Vertex range: %d; before: %d %d (range) %d %d (acceptance)", iEta, eta, vertexPos, vertexBinBegin, vertexBinEnd, vtxBegin, vtxEnd);
496       vertexBinBegin = TMath::Max(vertexBinBegin, vtxBegin);
497       vertexBinEnd =   TMath::Min(vertexBinEnd, vtxEnd);
498       //Printf(" after:  %d %d", vertexBinBegin, vertexBinEnd);
499
500       // no data for this bin
501       if (vertexBinBegin > vertexBinEnd)
502       {
503         //Printf("Bin empty. Skipped");
504         continue;
505       }
506
507       Float_t totalEvents = 0;
508       Float_t sum = 0;
509       Float_t sumError2 = 0;
510       Float_t unusedTracks = 0;
511       Float_t unusedEvents = 0;
512       for (Int_t iVtx = 1; iVtx <= vtxVsEta->GetNbinsX(); iVtx++)
513       {
514         if (iVtx >= vertexBinBegin && iVtx <= vertexBinEnd)
515         {
516           if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
517           {
518             sum += vtxVsEta->GetBinContent(iVtx, iEta);
519
520             if (sumError2 > 10e30)
521               Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
522
523             sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
524           }
525           totalEvents += vertexHist->GetBinContent(iVtx);
526         }
527         else
528         {
529           unusedTracks += vtxVsEta->GetBinContent(iVtx, iEta);
530           unusedEvents += vertexHist->GetBinContent(iVtx);
531         }
532       }
533
534       if (totalEvents == 0)
535       {
536         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
537         continue;
538       }
539
540       Float_t ptCutOffCorrection = 1;
541
542       // find pt cut off correction factor
543       if ((fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) && (fAnalysisMode & AliPWG0Helper::kFieldOn))
544       {
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         }
553       }
554
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);
556
557       Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
558       if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX())
559       {
560         Float_t dndeta = sum / totalEvents;
561         Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
562
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         
572         dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
573         error  = error / fdNdEta[vertexPos]->GetBinWidth(bin);
574
575         fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
576         fdNdEta[vertexPos]->SetBinError(bin, error);
577
578         dndeta /= ptCutOffCorrection;
579         error  /= ptCutOffCorrection;
580
581         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
582         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
583
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);
585       }
586     }
587   }
588 }
589
590 //____________________________________________________________________
591 Float_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
621 //____________________________________________________________________
622 void dNdEtaAnalysis::SaveHistograms()
623 {
624   // save the histograms to a directory with the name of this class (retrieved from TNamed)
625
626   gDirectory->mkdir(GetName());
627   gDirectory->cd(GetName());
628
629   if (fData)
630   {
631     fData->SaveHistograms();
632   }
633   else
634     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
635
636   if (fMult)
637   {
638     fMult->Write();
639   }
640   else
641     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n");
642
643   if (fPtDist)
644     fPtDist       ->Write();
645   else
646     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
647
648   for (Int_t i=0; i<kVertexBinning; ++i)
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   }
660
661   TNamed named("fTag", fTag.Data());
662   named.Write();
663
664   TParameter<Int_t> param("fAnalysisMode", fAnalysisMode);
665   param.Write();
666
667   gDirectory->cd("../");
668 }
669
670 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
671 {
672   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
673
674   if (!dir)
675     dir = GetName();
676
677   gDirectory->cd(dir);
678
679   fData->LoadHistograms();
680   fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName()));
681
682   for (Int_t i=0; i<kVertexBinning; ++i)
683   {
684     fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
685     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
686   }
687
688   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
689
690   if (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))
691     fTag = (dynamic_cast<TNamed*> (gDirectory->Get("fTag")))->GetTitle();
692
693   if (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))
694     fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gDirectory->Get("fAnalysisMode")))->GetVal();
695
696   gDirectory->cd("../");
697 }
698
699 //____________________________________________________________________
700 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
701 {
702   // draws the histograms
703
704   if (!simple)
705   {
706     if (fData)
707       fData->DrawHistograms(GetName());
708
709     TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
710     canvas->Divide(2, 1);
711
712     canvas->cd(1);
713     if (fdNdEtaPtCutOffCorrected[0])
714       fdNdEtaPtCutOffCorrected[0]->DrawCopy();
715
716     if (fdNdEta[0])
717     {
718       fdNdEta[0]->SetLineColor(kRed);
719       fdNdEta[0]->DrawCopy("SAME");
720     }
721
722     canvas->cd(2);
723     if (fPtDist)
724       fPtDist->DrawCopy();
725   }
726
727     // histograms for different vertices?
728   if (kVertexBinning > 0)
729   {
730       // doesnt work, but i dont get it, giving up...
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
736     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
737     //printf("%d\n", yPads);
738     //canvas2->Divide(2, yPads);
739
740     TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
741
742     for (Int_t i=0; i<kVertexBinning; ++i)
743     {
744       if (fdNdEtaPtCutOffCorrected[i])
745       {
746         canvas2->cd();
747
748         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
749         fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
750         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
751       }
752       if (canvas3 && fdNdEta[i])
753       {
754         canvas3->cd();
755
756         fdNdEta[i]->SetLineColor(i+1);
757         fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
758       }
759     }
760
761     canvas2->cd();
762     legend->Draw();
763     canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
764
765     if (canvas3)
766     {
767       canvas3->cd();
768       legend->Draw();
769     }
770   }
771
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"));
776
777      if (clone && clone2)
778      {
779         TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
780
781         clone->Divide(fdNdEtaPtCutOffCorrected[0]);
782         clone->GetYaxis()->SetRangeUser(0.95, 1.05);
783         clone->DrawCopy();
784
785         clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
786         clone2->DrawCopy("SAME");
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    }
794 }
795
796 Long64_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
812   const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning
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
824     collections[0]->Add(entry->fData);
825     collections[1]->Add(entry->fMult);
826     collections[2]->Add(entry->fPtDist);
827
828     for (Int_t i=0; i<kVertexBinning; ++i)
829     {
830       collections[3+i]->Add(entry->fdNdEta[i]);
831       collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
832     }
833
834     ++count;
835   }
836
837   fData->Merge(collections[0]);
838   fMult->Merge(collections[1]);
839   fPtDist->Merge(collections[2]);
840
841   for (Int_t i=0; i<kVertexBinning; ++i)
842   {
843     fdNdEta[i]->Merge(collections[3+i]);
844     fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]);
845   }
846
847   for (Int_t i=0; i<nCollections; ++i)
848     delete collections[i];
849
850   return count+1;
851 }
852
853