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