PWG0/dNdEta/dNdEtaAnalysis.cxx

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