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