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