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 }