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