2c629b0c7d1e1beadd6dfd487ddeb4e5a4ec139d
[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 "AlidNdEtaCorrection.h"
19 #include <AliCorrection.h>
20 #include <AliPWG0Helper.h>
21 #include <AliCorrectionMatrix2D.h>
22 #include <AliCorrectionMatrix3D.h>
23
24 //____________________________________________________________________
25 ClassImp(dNdEtaAnalysis)
26
27 //____________________________________________________________________
28 dNdEtaAnalysis::dNdEtaAnalysis() :
29   TNamed(),
30   fData(0),
31   fMult(0),
32   fPtDist(0),
33   fAnalysisMode(AliPWG0Helper::kInvalid),
34   fTag()
35 {
36   // default constructor
37
38   for (Int_t i=0; i<kVertexBinning; ++i)
39   {
40     fdNdEta[i] = 0;
41     fdNdEtaPtCutOffCorrected[i] = 0;
42   }
43 }
44
45 //____________________________________________________________________
46 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, AliPWG0Helper::AnalysisMode analysisMode) :
47   TNamed(name, title),
48   fData(0),
49   fMult(0),
50   fPtDist(0),
51   fAnalysisMode(analysisMode),
52   fTag()
53 {
54   // constructor
55
56   // TODO this binning has to be the same than in AliCorrection, somehow passed?!
57   Float_t binLimitsPt[] = {0.0, 0.05, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
58
59   fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysisMode);
60
61   // do not add this hists to the directory
62   Bool_t oldStatus = TH1::AddDirectoryStatus();
63   TH1::AddDirectory(kFALSE);
64
65   fMult = new TH1F("TriggeredMultiplicity", "Triggered Events;raw multiplicity;entries", 1000, -0.5, 999.5);
66
67   // TODO get binning from AliCorrection
68   fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 20, -2, 2);
69
70   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
71
72   for (Int_t i=1; i<kVertexBinning; ++i)
73   {
74     fdNdEta[i]    = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
75     fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
76   }
77
78   fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", 28, binLimitsPt);
79
80   TH1::AddDirectory(oldStatus);
81 }
82
83 //____________________________________________________________________
84 dNdEtaAnalysis::~dNdEtaAnalysis()
85 {
86   // destructor
87
88   if (fData)
89   {
90     delete fData;
91     fData = 0;
92   }
93
94   if (fMult)
95   {
96     delete fMult;
97     fMult = 0;
98   }
99
100   for (Int_t i=0; i<kVertexBinning; ++i)
101   {
102     if (fdNdEta[i])
103     {
104       delete fdNdEta[i];
105       fdNdEta[i] = 0;
106     }
107     if (fdNdEtaPtCutOffCorrected[i])
108     {
109       delete fdNdEtaPtCutOffCorrected[i];
110       fdNdEtaPtCutOffCorrected[i] = 0;
111     }
112   }
113
114   if (fPtDist)
115   {
116     delete fPtDist;
117     fPtDist = 0;
118   }
119 }
120
121 //_____________________________________________________________________________
122 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
123   TNamed(c),
124   fData(0),
125   fMult(0),
126   fPtDist(0),
127   fAnalysisMode(AliPWG0Helper::kInvalid),
128   fTag()
129 {
130   //
131   // dNdEtaAnalysis copy constructor
132   //
133
134   ((dNdEtaAnalysis &) c).Copy(*this);
135 }
136
137 //_____________________________________________________________________________
138 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
139 {
140   //
141   // Assignment operator
142   //
143
144   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
145   return *this;
146 }
147
148 //_____________________________________________________________________________
149 void dNdEtaAnalysis::Copy(TObject &c) const
150 {
151   //
152   // Copy function
153   //
154
155   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
156
157   target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
158   target.fMult = dynamic_cast<TH1F*> (fMult->Clone());
159
160   for (Int_t i=0; i<kVertexBinning; ++i)
161   {
162     target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
163     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
164   }
165
166   target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
167
168   target.fAnalysisMode = fAnalysisMode;
169   target.fTag = fTag;
170
171   TNamed::Copy((TNamed &) c);
172 }
173
174 //____________________________________________________________________
175 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
176 {
177   // fills a track into the histograms
178
179   fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
180 }
181
182 //____________________________________________________________________
183 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
184 {
185   // fills an event into the histograms
186
187   fData->GetEventCorrection()->FillMeas(vtx, n);
188 }
189
190 //____________________________________________________________________
191 void dNdEtaAnalysis::FillTriggeredEvent(Float_t n)
192 {
193   // fills a triggered event into the histograms
194
195   fMult->Fill(n);
196 }
197
198 //____________________________________________________________________
199 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType, const char* tag)
200 {
201   //
202   // correct with the given correction values and calculate dNdEta and pT distribution
203   // the corrections that are applied can be steered by the flag correctionType
204   // the measured result is not used up to a multiplicity of multCut (the bin at multCut is the first that is used!)
205   //
206
207   fTag.Form("Correcting dN/deta spectrum >>> %s <<<. Correction type: %d, pt cut: %.2f.", tag, (Int_t) correctionType, ptCut);
208   Printf("\n\n%s", fTag.Data());
209
210   // set corrections to 1
211   fData->SetCorrectionToUnity();
212
213   if (correction && correctionType != AlidNdEtaCorrection::kNone)
214   {
215     TH3F* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
216     TH2F* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
217
218     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
219       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
220
221     if (correctionType >= AlidNdEtaCorrection::kVertexReco)
222     {
223       trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
224       eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
225
226       // set bin with multiplicity 0 to unity (correction has no meaning in this bin)
227       for (Int_t i=0; i<=eventCorr->GetNbinsX()+1; i++)
228         eventCorr->SetBinContent(i, 1, 1);
229     }
230
231     switch (correctionType)
232     {
233       case AlidNdEtaCorrection::kINEL :
234       {
235         trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram());
236         eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram());
237         break;
238       }
239       case AlidNdEtaCorrection::kNSD :
240       {
241         trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram());
242         eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram());
243         break;
244       }
245       case AlidNdEtaCorrection::kND :
246       {
247         trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram());
248         eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
249         break;
250       }
251       default : break;
252     }
253   }
254   else
255     printf("INFO: No correction applied\n");
256
257   fData->Multiply();
258
259   if (correctionType >= AlidNdEtaCorrection::kVertexReco)
260   {
261     // There are no events with vertex that have 0 multiplicity, therefore
262     //   populate bin with 0 multiplicity with the following idea:
263     //     alpha = triggered events with vertex at a given vertex position / all triggered events with vertex
264     //     triggered events without vertex and 0 multiplicity at a given vertex position = alpha * all triggered events with 0 multiplicity
265     //   afterwards we still correct for the trigger efficiency
266
267     TH2* measuredEvents = fData->GetEventCorrection()->GetMeasuredHistogram();
268     TH2* correctedEvents = fData->GetEventCorrection()->GetGeneratedHistogram();
269
270     //new TCanvas; correctedEvents->DrawCopy("TEXT");
271
272     // start above 0 mult. bin with integration
273     TH1* vertexDist = measuredEvents->ProjectionX("vertexdist_measured", 2, measuredEvents->GetNbinsY()+1);
274     //new TCanvas; vertexDist->DrawCopy();
275
276     Int_t allEventsWithVertex = (Int_t) vertexDist->Integral(0, vertexDist->GetNbinsX()+1); // include under/overflow!
277     Int_t triggeredEventsWith0Mult = (Int_t) fMult->GetBinContent(1);
278
279     Printf("%d triggered events with 0 mult. -- %d triggered events with vertex", triggeredEventsWith0Mult, allEventsWithVertex);
280
281     for (Int_t i = 1; i <= measuredEvents->GetNbinsX(); i++)
282     {
283       Double_t alpha = (Double_t) vertexDist->GetBinContent(i) / allEventsWithVertex;
284       Double_t events = alpha * triggeredEventsWith0Mult;
285
286       // multiply with trigger correction if set above
287       events *= fData->GetEventCorrection()->GetCorrectionHistogram()->GetBinContent(i, 1);
288
289       Printf("Bin %d, alpha is %.2f, number of events with 0 mult.: %.2f", i, alpha * 100., events);
290
291       correctedEvents->SetBinContent(i, 1, events);
292     }
293
294     //new TCanvas; correctedEvents->DrawCopy("TEXT");
295   }
296
297   fData->PrintInfo(ptCut);
298
299   TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
300
301   // integrate multiplicity axis out (include under/overflow bins!!!)
302   TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
303
304   TH1D* vertexHist = (TH1D*) tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
305
306   // create pt hist
307   {
308     // reset all ranges
309     dataHist->GetXaxis()->SetRange(0, 0);
310     dataHist->GetYaxis()->SetRange(0, 0);
311     dataHist->GetZaxis()->SetRange(0, 0);
312
313     // vtx cut
314     Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5);
315     Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5);
316     dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
317     Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
318
319     if (nEvents > 0)
320     {
321       // eta cut
322       dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8));
323       Float_t etaWidth = 1.6;
324
325       TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
326
327       for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
328       {
329         Float_t binSize = fPtDist->GetBinWidth(i);
330         fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
331         fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
332       }
333
334       delete ptHist;
335     }
336     else
337       printf("ERROR: nEvents is 0!\n");
338   }
339
340   // reset all ranges
341   dataHist->GetXaxis()->SetRange(0, 0);
342   dataHist->GetYaxis()->SetRange(0, 0);
343   dataHist->GetZaxis()->SetRange(0, 0);
344
345   // integrate over pt (with pt cut)
346   Int_t ptLowBin = 1;
347   if (ptCut > 0)
348     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
349
350   dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
351   printf("pt range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins()+1);
352   TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
353
354   dataHist->GetZaxis()->SetRange(0, 0);
355   vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
356   vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
357
358   if (vtxVsEta == 0)
359   {
360     printf("ERROR: pt integration failed\n");
361     return;
362   }
363
364   const Float_t vertexRangeBegin[kVertexBinning] = { -9.99, -9.99, 0 };
365   const Float_t vertexRangeEnd[kVertexBinning]   = { 9.99,  0,     9.99 };
366
367   // TODO acceptance range binning-independent; only for SPD
368   const Int_t binBegin[20] = { 11, 10, 9, 7, 5, 4, 3, 3, 3, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1 };
369
370   for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
371   {
372     // loop over vertex ranges
373     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
374     {
375       Int_t vertexBinBegin = vertexHist->GetXaxis()->FindBin(vertexRangeBegin[vertexPos]);
376       Int_t vertexBinEnd   = vertexHist->GetXaxis()->FindBin(vertexRangeEnd[vertexPos]);
377
378       // adjust acceptance range
379
380       if (fAnalysisMode == AliPWG0Helper::kSPD)
381       {
382         //printf("Eta bin: %d Vertex range: %d; before: %d %d ", iEta, vertexPos, vertexBinBegin, vertexBinEnd);
383         vertexBinBegin = TMath::Max(vertexBinBegin, binBegin[iEta - 1]);
384         vertexBinEnd =   TMath::Min(vertexBinEnd, 15 - binBegin[20 - iEta]);
385         //Printf("after:  %d %d", vertexBinBegin, vertexBinEnd);
386       }
387
388       // no data for this bin
389       if (vertexBinBegin > vertexBinEnd)
390         continue;
391
392       Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
393       if (totalEvents == 0)
394       {
395         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
396         continue;
397       }
398
399       Float_t sum = 0;
400       Float_t sumError2 = 0;
401       for (Int_t iVtx = vertexBinBegin; iVtx <= vertexBinEnd; iVtx++)
402       {      
403         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
404         {
405           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
406
407           if (sumError2 > 10e30)
408             Printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
409
410           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
411         }
412       }
413
414       Float_t ptCutOffCorrection = 1;
415       if (correction && ptCut > 0)
416         ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
417
418       if (ptCutOffCorrection <= 0)
419       {
420         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
421         continue;
422       }
423
424       //printf("Eta: %d Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f\n", iEta, vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
425
426       Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
427       if (bin > 0 && bin <= fdNdEta[vertexPos]->GetNbinsX())
428       {
429         Float_t dndeta = sum / totalEvents;
430         Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
431
432         dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
433         error  = error / fdNdEta[vertexPos]->GetBinWidth(bin);
434
435         fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
436         fdNdEta[vertexPos]->SetBinError(bin, error);
437
438         dndeta /= ptCutOffCorrection;
439         error  /= ptCutOffCorrection;
440
441         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
442         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
443
444         //Printf("Bin %d has dN/deta = %f", bin, dndeta);
445       }
446     }
447   }
448 }
449
450 //____________________________________________________________________
451 void dNdEtaAnalysis::SaveHistograms()
452 {
453   // save the histograms to a directory with the name of this class (retrieved from TNamed)
454
455   gDirectory->mkdir(GetName());
456   gDirectory->cd(GetName());
457
458   if (fData)
459   {
460     fData->SaveHistograms();
461   }
462   else
463     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
464
465   if (fMult)
466   {
467     fMult->Write();
468   }
469   else
470     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fMult is 0\n");
471
472   if (fPtDist)
473     fPtDist       ->Write();
474   else
475     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
476
477   for (Int_t i=0; i<kVertexBinning; ++i)
478   {
479     if (fdNdEta[i])
480       fdNdEta[i]->Write();
481     else
482       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
483
484     if (fdNdEtaPtCutOffCorrected[i])
485       fdNdEtaPtCutOffCorrected[i]->Write();
486     else
487       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
488   }
489
490   TNamed named("fTag", fTag.Data());
491   named.Write();
492
493   TParameter<Int_t> param("fAnalysisMode", fAnalysisMode);
494   param.Write();
495
496   gDirectory->cd("../");
497 }
498
499 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
500 {
501   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
502
503   if (!dir)
504     dir = GetName();
505
506   gDirectory->cd(dir);
507
508   fData->LoadHistograms();
509   fMult = dynamic_cast<TH1F*> (gDirectory->Get(fMult->GetName()));
510
511   for (Int_t i=0; i<kVertexBinning; ++i)
512   {
513     fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
514     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
515   }
516
517   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
518
519   if (dynamic_cast<TNamed*> (gFile->Get("fTag")))
520     fTag = (dynamic_cast<TNamed*> (gFile->Get("fTag")))->GetTitle();
521
522   if (dynamic_cast<TParameter<Int_t>*> (gFile->Get("fAnalysisMode")))
523     fAnalysisMode = (AliPWG0Helper::AnalysisMode) (dynamic_cast<TParameter<Int_t>*> (gFile->Get("fAnalysisMode")))->GetVal();
524
525   gDirectory->cd("../");
526 }
527
528 //____________________________________________________________________
529 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
530 {
531   // draws the histograms
532
533   if (!simple)
534   {
535     if (fData)
536       fData->DrawHistograms(GetName());
537
538     TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
539     canvas->Divide(2, 1);
540
541     canvas->cd(1);
542     if (fdNdEtaPtCutOffCorrected[0])
543       fdNdEtaPtCutOffCorrected[0]->DrawCopy();
544
545     if (fdNdEta[0])
546     {
547       fdNdEta[0]->SetLineColor(kRed);
548       fdNdEta[0]->DrawCopy("SAME");
549     }
550
551     canvas->cd(2);
552     if (fPtDist)
553       fPtDist->DrawCopy();
554   }
555
556     // histograms for different vertices?
557   if (kVertexBinning > 0)
558   {
559       // doesnt work, but i dont get it, giving up...
560     TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
561     TCanvas* canvas3 = 0;
562     if (!simple)
563       canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
564
565     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
566     //printf("%d\n", yPads);
567     //canvas2->Divide(2, yPads);
568
569     TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
570
571     for (Int_t i=0; i<kVertexBinning; ++i)
572     {
573       if (fdNdEtaPtCutOffCorrected[i])
574       {
575         canvas2->cd();
576
577         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
578         fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
579         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
580       }
581       if (canvas3 && fdNdEta[i])
582       {
583         canvas3->cd();
584
585         fdNdEta[i]->SetLineColor(i+1);
586         fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
587       }
588     }
589
590     canvas2->cd();
591     legend->Draw();
592     canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
593
594     if (canvas3)
595     {
596       canvas3->cd();
597       legend->Draw();
598     }
599   }
600
601   if (kVertexBinning == 3)
602   {
603      TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
604      TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
605
606      if (clone && clone2)
607      {
608         TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
609
610         clone->Divide(fdNdEtaPtCutOffCorrected[0]);
611         clone->GetYaxis()->SetRangeUser(0.95, 1.05);
612         clone->DrawCopy();
613
614         clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
615         clone2->DrawCopy("SAME");
616
617         TLine* line = new TLine(-1, 1, 1, 1);
618         line->Draw();
619
620         canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
621      }
622    }
623 }
624
625 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
626 {
627   // Merges a list of dNdEtaAnalysis objects with this one.
628   // This is needed for PROOF.
629   // Returns the number of merged objects (including this)
630
631   if (!list)
632     return 0;
633
634   if (list->IsEmpty())
635     return 1;
636
637   TIterator* iter = list->MakeIterator();
638   TObject* obj;
639
640   // sub collections
641   const Int_t nCollections = 2 * kVertexBinning + 3; // 3 standalone hists, 3 arrays of size kVertexBinning
642   TList* collections[nCollections];
643   for (Int_t i=0; i<nCollections; ++i)
644     collections[i] = new TList;
645
646   Int_t count = 0;
647   while ((obj = iter->Next()))
648   {
649     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
650     if (entry == 0)
651       continue;
652
653     collections[0]->Add(entry->fData);
654     collections[1]->Add(entry->fMult);
655     collections[2]->Add(entry->fPtDist);
656
657     for (Int_t i=0; i<kVertexBinning; ++i)
658     {
659       collections[3+i]->Add(entry->fdNdEta[i]);
660       collections[3+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
661     }
662
663     ++count;
664   }
665
666   fData->Merge(collections[0]);
667   fMult->Merge(collections[1]);
668   fPtDist->Merge(collections[2]);
669
670   for (Int_t i=0; i<kVertexBinning; ++i)
671   {
672     fdNdEta[i]->Merge(collections[3+i]);
673     fdNdEtaPtCutOffCorrected[i]->Merge(collections[3+kVertexBinning+i]);
674   }
675
676   for (Int_t i=0; i<nCollections; ++i)
677     delete collections[i];
678
679   return count+1;
680 }