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