PWG0/dNdEta/dNdEtaAnalysis.cxx

   1 /* $Id$ */
   2
   3 #include "dNdEtaAnalysis.h"
   4
   5 #include <TFile.h>
   6 #include <TH3F.h>
   7 #include <TH2D.h>
   8 #include <TH1D.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
  16 #include "AlidNdEtaCorrection.h"
  17 #include "AliPWG0Helper.h"
  18
  19 //____________________________________________________________________
  20 ClassImp(dNdEtaAnalysis)
  21
  22 //____________________________________________________________________
  23 dNdEtaAnalysis::dNdEtaAnalysis() :
  24   TNamed(),
  25   fData(0),
  26   fDataUncorrected(0),
  27   fVtx(0),
  28   fPtDist(0)
  29 {
  30   // default constructor
  31
  32   for (Int_t i=0; i<kVertexBinning; ++i)
  33   {
  34     fdNdEta[i] = 0;
  35     fdNdEtaPtCutOffCorrected[i] = 0;
  36   }
  37 }
  38
  39 //____________________________________________________________________
  40 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
  41   TNamed(name, title),
  42   fData(0),
  43   fDataUncorrected(0),
  44   fVtx(0),
  45   fPtDist(0)
  46 {
  47   // constructor
  48
  49   fData  = new TH3F(Form("%s_analysis", name),"Input data",80,-20,20,40,-2,2,100, 0, 10);
  50   fData->SetXTitle("vtx z [cm]");
  51   fData->SetYTitle("#eta");
  52   fData->SetZTitle("p_{T}");
  53   fData->Sumw2();
  54
  55   fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
  56   fDataUncorrected->SetTitle(Form("%s uncorrected", fData->GetTitle()));
  57   fVtx       = dynamic_cast<TH1D*> (fData->Project3D("x"));
  58   fVtx->SetName(Form("%s_vtx", name));
  59   fVtx->SetTitle("Vertex distribution");
  60   fVtx->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
  61   //fVtx->Sumw2();
  62
  63   fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
  64   fdNdEta[0]->SetName(Form("%s_dNdEta", name));
  65   fdNdEta[0]->SetTitle("dN_{ch}/d#eta");
  66   fdNdEta[0]->GetXaxis()->SetTitle(fData->GetYaxis()->GetTitle());
  67   fdNdEta[0]->SetYTitle("dN_{ch}/d#eta");
  68
  69   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
  70
  71   for (Int_t i=1; i<kVertexBinning; ++i)
  72   {
  73     fdNdEta[i]    = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
  74     fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
  75   }
  76
  77   fPtDist = dynamic_cast<TH1D*> (fData->Project3D("z"));
  78   fPtDist->SetName(Form("%s_pt", name));
  79   fPtDist->SetTitle("p_{T} [GeV/c]");
  80   fPtDist->GetXaxis()->SetTitle(fData->GetZaxis()->GetTitle());
  81   fPtDist->SetYTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
  82 }
  83
  84 //____________________________________________________________________
  85 dNdEtaAnalysis::~dNdEtaAnalysis()
  86 {
  87   // destructor
  88
  89   if (fData)
  90   {
  91     delete fData;
  92     fData = 0;
  93   }
  94
  95   if (fDataUncorrected)
  96   {
  97     delete fDataUncorrected;
  98     fDataUncorrected = 0;
  99   }
 100
 101   if (fVtx)
 102   {
 103     delete fVtx;
 104     fVtx = 0;
 105   }
 106
 107   for (Int_t i=0; i<kVertexBinning; ++i)
 108   {
 109     if (fdNdEta[i])
 110     {
 111       delete fdNdEta[i];
 112       fdNdEta[i] = 0;
 113     }
 114     if (fdNdEtaPtCutOffCorrected[i])
 115     {
 116       delete fdNdEtaPtCutOffCorrected[i];
 117       fdNdEtaPtCutOffCorrected[i] = 0;
 118     }
 119   }
 120
 121   if (fPtDist)
 122   {
 123     delete fPtDist;
 124     fPtDist = 0;
 125   }
 126 }
 127
 128 //_____________________________________________________________________________
 129 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
 130   TNamed(c),
 131   fData(0),
 132   fDataUncorrected(0),
 133   fVtx(0),
 134   fPtDist(0)
 135 {
 136   //
 137   // dNdEtaAnalysis copy constructor
 138   //
 139
 140   ((dNdEtaAnalysis &) c).Copy(*this);
 141 }
 142
 143 //_____________________________________________________________________________
 144 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
 145 {
 146   //
 147   // Assignment operator
 148   //
 149
 150   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
 151   return *this;
 152 }
 153
 154 //_____________________________________________________________________________
 155 void dNdEtaAnalysis::Copy(TObject &c) const
 156 {
 157   //
 158   // Copy function
 159   //
 160
 161   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
 162
 163   target.fData = dynamic_cast<TH3F*> (fData->Clone());
 164   target.fDataUncorrected = dynamic_cast<TH3F*> (fDataUncorrected->Clone());
 165   target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
 166
 167   for (Int_t i=0; i<kVertexBinning; ++i)
 168   {
 169     target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
 170     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[i]->Clone());
 171   }
 172
 173   target.fPtDist = dynamic_cast<TH1D*> (fPtDist->Clone());
 174
 175   TNamed::Copy((TNamed &) c);
 176 }
 177
 178 //____________________________________________________________________
 179 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight)
 180 {
 181   // fills a track into the histograms
 182
 183   fDataUncorrected->Fill(vtx, eta, pt);
 184   fData->Fill(vtx, eta, pt, weight);
 185 }
 186
 187 //____________________________________________________________________
 188 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
 189 {
 190   // fills an event into the histograms
 191
 192   fVtx->Fill(vtx, weight);
 193 }
 194
 195 //____________________________________________________________________
 196 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
 197 {
 198   // correct with correction values if available
 199
 200   // TODO what do we do with the error?
 201   if (!correction)
 202     printf("INFO: No correction applied\n");
 203
 204   // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
 205
 206   // create pt hist
 207   {
 208     // reset all ranges
 209     fData->GetXaxis()->SetRange(0, 0);
 210     fData->GetYaxis()->SetRange(0, 0);
 211     fData->GetZaxis()->SetRange(0, 0);
 212
 213     // vtx cut
 214     Int_t vertexBinBegin = fData->GetXaxis()->FindBin(-5);
 215     Int_t vertexBinEnd = fData->GetXaxis()->FindBin(5);
 216     fData->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
 217     Float_t nEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd);
 218
 219     // eta cut
 220     fData->GetYaxis()->SetRange(fData->GetYaxis()->FindBin(-0.8), fData->GetYaxis()->FindBin(0.8));
 221     Float_t etaWidth = 1.6;
 222
 223     TH1D* ptHist = dynamic_cast<TH1D*> (fData->Project3D("ze"));
 224
 225     for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
 226     {
 227       Float_t binSize = fPtDist->GetBinWidth(i);
 228       fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
 229       fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
 230     }
 231
 232     delete ptHist;
 233   }
 234
 235   // reset all ranges
 236   fData->GetXaxis()->SetRange(0, 0);
 237   fData->GetYaxis()->SetRange(0, 0);
 238   fData->GetZaxis()->SetRange(0, 0);
 239
 240   // integrate over pt (with pt cut)
 241   Int_t ptLowBin = 1;
 242   if (ptCut > 0)
 243     ptLowBin = fData->GetZaxis()->FindBin(ptCut);
 244
 245   fData->GetZaxis()->SetRange(ptLowBin, fData->GetZaxis()->GetNbins());
 246   TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2e"));
 247   vtxVsEta->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
 248   vtxVsEta->GetYaxis()->SetTitle(fData->GetYaxis()->GetTitle());
 249
 250   if (vtxVsEta == 0)
 251   {
 252     printf("ERROR: pt integration failed\n");
 253     return;
 254   }
 255
 256   //new TCanvas;
 257   //vtxVsEta->Draw("COLZ");
 258
 259   for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
 260   {
 261     // do we have several histograms for different vertex positions?
 262     Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
 263     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
 264     {
 265       Int_t vertexBinBegin = 1;
 266       Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
 267
 268       // the first histogram is always for the whole vertex range
 269       if (vertexPos > 0)
 270       {
 271         vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
 272         vertexBinEnd = vertexBinBegin + vertexBinWidth;
 273       }
 274
 275       Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
 276       if (totalEvents == 0)
 277       {
 278         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
 279         continue;
 280       }
 281
 282       Float_t sum = 0;
 283       Float_t sumError2 = 0;
 284       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
 285       {
 286         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
 287         {
 288           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
 289           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
 290         }
 291       }
 292
 293       Float_t ptCutOffCorrection = 1;
 294       if (correction)
 295         ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
 296
 297       if (ptCutOffCorrection <= 0)
 298       {
 299         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
 300         continue;
 301       }
 302
 303       Float_t dndeta = sum / totalEvents;
 304       Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
 305
 306       dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
 307       error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
 308
 309       fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
 310       fdNdEta[vertexPos]->SetBinError(iEta, error);
 311
 312       dndeta /= ptCutOffCorrection;
 313       error  /= ptCutOffCorrection;
 314
 315       fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(iEta, dndeta);
 316       fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(iEta, error);
 317     }
 318   }
 319 }
 320
 321 //____________________________________________________________________
 322 void dNdEtaAnalysis::SaveHistograms()
 323 {
 324   // save the histograms to a directory with the name of this class (retrieved from TNamed)
 325
 326   gDirectory->mkdir(GetName());
 327   gDirectory->cd(GetName());
 328
 329   if (fData)
 330   {
 331     fData->Write();
 332     AliPWG0Helper::CreateProjections(fData);
 333   }
 334   else
 335     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
 336
 337   if (fDataUncorrected)
 338   {
 339     fDataUncorrected->Write();
 340     AliPWG0Helper::CreateProjections(fDataUncorrected);
 341   }
 342   else
 343     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fDataUncorrected is 0\n");
 344
 345   if (fVtx)
 346     fVtx       ->Write();
 347   else
 348     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fVtx is 0\n");
 349
 350   if (fPtDist)
 351     fPtDist       ->Write();
 352   else
 353     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
 354
 355   for (Int_t i=0; i<kVertexBinning; ++i)
 356   {
 357     if (fdNdEta[i])
 358       fdNdEta[i]->Write();
 359     else
 360       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
 361
 362     if (fdNdEtaPtCutOffCorrected[i])
 363       fdNdEtaPtCutOffCorrected[i]->Write();
 364     else
 365       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
 366   }
 367
 368   gDirectory->cd("../");
 369 }
 370
 371 void dNdEtaAnalysis::LoadHistograms()
 372 {
 373   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
 374
 375   gDirectory->cd(GetName());
 376
 377   fData = dynamic_cast<TH3F*> (gDirectory->Get(fData->GetName()));
 378   fDataUncorrected = dynamic_cast<TH3F*> (gDirectory->Get(fDataUncorrected->GetName()));
 379
 380   fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
 381
 382   for (Int_t i=0; i<kVertexBinning; ++i)
 383   {
 384     fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
 385     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
 386   }
 387
 388   fPtDist = dynamic_cast<TH1D*> (gDirectory->Get(fPtDist->GetName()));
 389
 390   gDirectory->cd("../");
 391 }
 392
 393 //____________________________________________________________________
 394 void dNdEtaAnalysis::DrawHistograms()
 395 {
 396   // draws the histograms
 397
 398   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
 399   canvas->Divide(2, 2);
 400
 401   canvas->cd(1);
 402   if (fData)
 403     fData->Draw("COLZ");
 404
 405   /*canvas->cd(2);
 406   if (fDataUncorrected)
 407     fDataUncorrected->Draw("COLZ");*/
 408
 409   canvas->cd(2);
 410   if (fVtx)
 411     fVtx->Draw();
 412
 413   canvas->cd(3);
 414   if (fdNdEtaPtCutOffCorrected[0])
 415     fdNdEtaPtCutOffCorrected[0]->Draw();
 416
 417   if (fdNdEta[0])
 418   {
 419     fdNdEta[0]->SetLineColor(kRed);
 420     fdNdEta[0]->Draw("SAME");
 421   }
 422
 423   canvas->cd(4);
 424   if (fPtDist)
 425     fPtDist->Draw();
 426
 427     // histograms for different vertices?
 428   if (kVertexBinning > 0)
 429   {
 430       // doesnt work, but i dont get it, giving up...
 431     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
 432     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
 433     //printf("%d\n", yPads);
 434     //canvas2->Divide(2, yPads);
 435
 436     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
 437
 438     for (Int_t i=0; i<kVertexBinning; ++i)
 439     {
 440       //canvas2->cd(i-1);
 441       //printf("%d\n", i);
 442       if (fdNdEtaPtCutOffCorrected[i])
 443       {
 444         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
 445         fdNdEtaPtCutOffCorrected[i]->Draw((i == 0) ? "" : "SAME");
 446         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
 447       }
 448     }
 449
 450     legend->Draw();
 451   }
 452 }
 453
 454 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
 455 {
 456   // Merges a list of dNdEtaAnalysis objects with this one.
 457   // This is needed for PROOF.
 458   // Returns the number of merged objects (including this)
 459
 460   if (!list)
 461     return 0;
 462
 463   if (list->IsEmpty())
 464     return 1;
 465
 466   TIterator* iter = list->MakeIterator();
 467   TObject* obj;
 468
 469   // sub collections
 470   const Int_t nCollections = 2 * kVertexBinning + 4; // 4 standalone hists, two arrays of size kVertexBinning
 471   TList* collections[nCollections];
 472   for (Int_t i=0; i<nCollections; ++i)
 473     collections[i] = new TList;
 474
 475   Int_t count = 0;
 476   while ((obj = iter->Next()))
 477   {
 478     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
 479     if (entry == 0)
 480       continue;
 481
 482     collections[0]->Add(entry->fData);
 483     if (fDataUncorrected)
 484       collections[1]->Add(entry->fDataUncorrected);
 485     collections[2]->Add(entry->fVtx);
 486     collections[3]->Add(entry->fPtDist);
 487
 488     for (Int_t i=0; i<kVertexBinning; ++i)
 489     {
 490       collections[4+i]->Add(entry->fdNdEta[i]);
 491       collections[4+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
 492     }
 493
 494     ++count;
 495   }
 496
 497   fData->Merge(collections[0]);
 498   if (fDataUncorrected)
 499     fDataUncorrected->Merge(collections[1]);
 500   fVtx->Merge(collections[2]);
 501   fPtDist->Merge(collections[3]);
 502   for (Int_t i=0; i<kVertexBinning; ++i)
 503   {
 504     fdNdEta[i]->Merge(collections[4+i]);
 505     fdNdEta[i]->Merge(collections[4+kVertexBinning+i]);
 506   }
 507
 508   for (Int_t i=0; i<nCollections; ++i)
 509     delete collections[i];
 510
 511   return count+1;
 512 }