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(Char_t* name, Char_t* title) :
  24   TNamed(name, title),
  25   fData(0),
  26   fDataUncorrected(0),
  27   fNEvents(0),
  28   fVtx(0)
  29 {
  30   // constructor
  31
  32   fData  = new TH3F(Form("%s_analysis", name),"dNdEtaAnalysis",80,-20,20,120,-6,6,100, 0, 10);
  33   fData->SetXTitle("vtx z [cm]");
  34   fData->SetYTitle("#eta");
  35   fData->SetZTitle("p_{T}");
  36
  37   fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
  38   fVtx       = dynamic_cast<TH1D*> (fData->Project3D("x"));
  39
  40   fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
  41   for (Int_t i=1; i<kVertexBinning; ++i)
  42   {
  43     fdNdEta[i]    = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
  44     fdNdEta[i]->SetYTitle("dN/d#eta");
  45   }
  46
  47   fData->Sumw2();
  48   fVtx->Sumw2();
  49 }
  50
  51 //____________________________________________________________________
  52 dNdEtaAnalysis::~dNdEtaAnalysis()
  53 {
  54   // destructor
  55
  56   delete fData;
  57   fData = 0;
  58
  59   delete fDataUncorrected;
  60   fDataUncorrected = 0;
  61
  62   delete fVtx;
  63   fVtx = 0;
  64
  65   for (Int_t i=0; i<kVertexBinning; ++i)
  66   {
  67     delete fdNdEta[i];
  68     fdNdEta[i] = 0;
  69   }
  70 }
  71
  72 //_____________________________________________________________________________
  73 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
  74   TNamed(c),
  75   fData(0),
  76   fDataUncorrected(0),
  77   fNEvents(0),
  78   fVtx(0)
  79 {
  80   //
  81   // dNdEtaAnalysis copy constructor
  82   //
  83
  84   ((dNdEtaAnalysis &) c).Copy(*this);
  85 }
  86
  87 //_____________________________________________________________________________
  88 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
  89 {
  90   //
  91   // Assignment operator
  92   //
  93
  94   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
  95   return *this;
  96 }
  97
  98 //_____________________________________________________________________________
  99 void dNdEtaAnalysis::Copy(TObject &c) const
 100 {
 101   //
 102   // Copy function
 103   //
 104
 105   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
 106
 107   target.fData = dynamic_cast<TH3F*> (fData->Clone());
 108   target.fDataUncorrected = dynamic_cast<TH3F*> (fDataUncorrected->Clone());
 109   target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
 110
 111   for (Int_t i=0; i<kVertexBinning; ++i)
 112     target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
 113
 114   target.fNEvents = fNEvents;
 115
 116   TNamed::Copy((TNamed &) c);
 117 }
 118
 119 //____________________________________________________________________
 120 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight)
 121 {
 122   // fills a track into the histograms
 123
 124   fDataUncorrected->Fill(vtx, eta, pt);
 125   fData->Fill(vtx, eta, pt, weight);
 126 }
 127
 128 //____________________________________________________________________
 129 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
 130 {
 131   // fills an event into the histograms
 132
 133   fVtx->Fill(vtx, weight);
 134
 135   fNEvents += weight;
 136 }
 137
 138 //____________________________________________________________________
 139 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
 140 {
 141   // correct with correction values if available
 142
 143   // TODO what do we do with the error?
 144   if (!correction)
 145     printf("INFO: No correction applied\n");
 146
 147   // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
 148
 149   // integrate over pt (with pt cut)
 150   fData->GetZaxis()->SetRange(fData->GetZaxis()->FindBin(ptCut), fData->GetZaxis()->GetNbins());
 151   TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2"));
 152   if (vtxVsEta == 0)
 153   {
 154     printf("ERROR: pt integration failed\n");
 155     return;
 156   }
 157
 158   new TCanvas;
 159   vtxVsEta->Draw("COLZ");
 160
 161   for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
 162   {
 163     // do we have several histograms for different vertex positions?
 164     Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
 165     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
 166     {
 167       Int_t vertexBinBegin = 1;
 168       Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
 169
 170       // the first histogram is always for the whole vertex range
 171       if (vertexPos > 0)
 172       {
 173         vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
 174         vertexBinEnd = vertexBinBegin + vertexBinWidth;
 175       }
 176
 177       Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
 178       if (totalEvents == 0)
 179       {
 180         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
 181         continue;
 182       }
 183
 184       Float_t sum = 0;
 185       Float_t sumError2 = 0;
 186       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
 187       {
 188         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
 189         {
 190           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
 191           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
 192         }
 193       }
 194
 195       Float_t ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
 196       //ptCutOffCorrection = 1;
 197       if (ptCutOffCorrection <= 0)
 198       {
 199         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
 200         continue;
 201       }
 202
 203       Float_t dndeta = sum / totalEvents / ptCutOffCorrection;
 204       Float_t error  = TMath::Sqrt(sumError2) / totalEvents / ptCutOffCorrection;
 205
 206       dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
 207       error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
 208
 209       fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
 210       fdNdEta[vertexPos]->SetBinError(iEta, error);
 211     }
 212   }
 213 }
 214
 215 //____________________________________________________________________
 216 void dNdEtaAnalysis::SaveHistograms()
 217 {
 218   // save the histograms to a directory with the name of this class (retrieved from TNamed)
 219
 220   gDirectory->mkdir(GetName());
 221   gDirectory->cd(GetName());
 222
 223   fData  ->Write();
 224   AliPWG0Helper::CreateProjections(fData);
 225   fDataUncorrected->Write();
 226   AliPWG0Helper::CreateProjections(fDataUncorrected);
 227
 228   fVtx       ->Write();
 229   for (Int_t i=0; i<kVertexBinning; ++i)
 230     fdNdEta[i]    ->Write();
 231
 232   gDirectory->cd("../");
 233 }
 234
 235 void dNdEtaAnalysis::LoadHistograms()
 236 {
 237   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
 238
 239   gDirectory->cd(GetName());
 240
 241   fData = dynamic_cast<TH3F*> (gDirectory->Get(fData->GetName()));
 242   fDataUncorrected = dynamic_cast<TH3F*> (gDirectory->Get(fDataUncorrected->GetName()));
 243
 244   fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
 245
 246   for (Int_t i=0; i<kVertexBinning; ++i)
 247     fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
 248
 249   gDirectory->cd("../");
 250 }
 251
 252 //____________________________________________________________________
 253 void dNdEtaAnalysis::DrawHistograms()
 254 {
 255   // draws the histograms
 256
 257   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
 258   canvas->Divide(2, 2);
 259
 260   canvas->cd(1);
 261   if (fData)
 262     fData->Draw("COLZ");
 263
 264   canvas->cd(2);
 265   if (fDataUncorrected)
 266     fDataUncorrected->Draw("COLZ");
 267
 268   canvas->cd(3);
 269   if (fVtx)
 270     fVtx->Draw();
 271
 272   canvas->cd(4);
 273   if (fdNdEta[0])
 274     fdNdEta[0]->Draw();
 275
 276     // histograms for different vertices?
 277   if (kVertexBinning > 0)
 278   {
 279       // doesnt work, but i dont get it, giving up...
 280     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
 281     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
 282     //printf("%d\n", yPads);
 283     //canvas2->Divide(2, yPads);
 284
 285     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
 286
 287     for (Int_t i=0; i<kVertexBinning; ++i)
 288     {
 289       //canvas2->cd(i-1);
 290       //printf("%d\n", i);
 291       if (fdNdEta[i])
 292       {
 293         fdNdEta[i]->SetLineColor(i+1);
 294         fdNdEta[i]->Draw((i == 0) ? "" : "SAME");
 295         legend->AddEntry(fdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
 296       }
 297     }
 298
 299     legend->Draw();
 300   }
 301 }
 302
 303 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
 304 {
 305   // Merges a list of dNdEtaAnalysis objects with this one.
 306   // This is needed for PROOF.
 307   // Returns the number of merged objects (including this)
 308
 309   if (!list)
 310     return 0;
 311
 312   if (list->IsEmpty())
 313     return 1;
 314
 315   TIterator* iter = list->MakeIterator();
 316   TObject* obj;
 317
 318   // sub collections
 319   const Int_t nCollections = kVertexBinning + 3;
 320   TList* collections[nCollections];
 321   for (Int_t i=0; i<nCollections; ++i)
 322     collections[i] = new TList;
 323
 324   Int_t count = 0;
 325   while ((obj = iter->Next()))
 326   {
 327     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
 328     if (entry == 0)
 329       continue;
 330
 331     collections[0]->Add(entry->fData);
 332     collections[1]->Add(entry->fDataUncorrected);
 333     collections[2]->Add(entry->fVtx);
 334
 335     for (Int_t i=0; i<kVertexBinning; ++i)
 336       collections[3+i]->Add(entry->fdNdEta[i]);
 337
 338     ++count;
 339   }
 340
 341   fData->Merge(collections[0]);
 342   fDataUncorrected->Merge(collections[1]);
 343   fVtx->Merge(collections[2]);
 344   for (Int_t i=0; i<kVertexBinning; ++i)
 345     fdNdEta[i]->Merge(collections[3+i]);
 346
 347   for (Int_t i=0; i<nCollections; ++i)
 348     delete collections[i];
 349
 350   return count+1;
 351 }