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
  14 //____________________________________________________________________
  15 ClassImp(dNdEtaAnalysis)
  16
  17 //____________________________________________________________________
  18 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
  19 TNamed(name, title)
  20 {
  21   hEtaVsVtx  = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
  22   hEtaVsVtx->SetXTitle("vtx z [cm]");
  23   hEtaVsVtx->SetYTitle("#eta");
  24
  25   hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
  26   hVtx       = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
  27   for (Int_t i=0; i<kVertexBinning; ++i)
  28   {
  29     hdNdEta[i]    = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
  30     hdNdEta[i]->SetYTitle("dN/d#eta");
  31   }
  32
  33   hEtaVsVtx->Sumw2();
  34   hVtx->Sumw2();
  35 }
  36
  37 //____________________________________________________________________
  38 void
  39 dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) {
  40   hEtaVsVtx->Fill(vtx, eta, weight);
  41   hEtaVsVtxUncorrected->Fill(vtx,eta);
  42 }
  43
  44 //____________________________________________________________________
  45 void
  46 dNdEtaAnalysis::FillEvent(Float_t vtx) {
  47   hVtx->Fill(vtx);
  48 }
  49
  50 //____________________________________________________________________
  51 void
  52 dNdEtaAnalysis::Finish() {
  53
  54   // first normalize with n events (per vtx)
  55   for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) {
  56     Float_t nEvents      = hVtx->GetBinContent(i_vtx);
  57     Float_t nEventsError = hVtx->GetBinError(i_vtx);
  58
  59     if (nEvents==0) continue;
  60
  61     for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
  62       Float_t value = hEtaVsVtx->GetBinContent(i_vtx, i_eta)/nEvents;
  63       if (value==0) continue;
  64       Float_t error = hEtaVsVtx->GetBinError(i_vtx, i_eta)/nEvents;
  65       error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta)/
  66                                        hEtaVsVtx->GetBinContent(i_vtx, i_eta),2) +
  67                           TMath::Power(nEventsError/nEvents,2));
  68       hEtaVsVtx->SetBinContent(i_vtx, i_eta, value);
  69       hEtaVsVtx->SetBinError(i_vtx,   i_eta, error);
  70     }
  71   }
  72
  73   // then take the wieghted average for each eta
  74   // is this the right way to do it???
  75   for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
  76     Float_t sum           = 0;
  77     Float_t sumError2     = 0;
  78     Int_t   nMeasurements = 0;
  79
  80     Float_t sum_xw = 0;
  81     Float_t sum_w  = 0;
  82
  83     // do we have several histograms for different vertex positions?
  84     Int_t vertexBinWidth = hVtx->GetNbinsX() / kVertexBinning;
  85     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
  86     {
  87       Int_t vertexBinBegin = vertexBinWidth * vertexPos;
  88       Int_t vertexBinEnd = vertexBinBegin + vertexBinWidth;
  89
  90       for (Int_t i_vtx=vertexBinBegin; i_vtx<=vertexBinEnd; i_vtx++) {
  91         if (hVtx->GetBinContent(i_vtx)==0)             continue;
  92         if (hEtaVsVtx->GetBinContent(i_vtx, i_eta)==0) continue;
  93
  94         Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
  95         sum_xw = sum_xw + hEtaVsVtx->GetBinContent(i_vtx, i_eta)*w;
  96         sum_w  = sum_w + w;
  97
  98         sum = sum + hEtaVsVtx->GetBinContent(i_vtx, i_eta);
  99         sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
 100         nMeasurements++;
 101       }
 102       Float_t dndeta = 0;
 103       Float_t error  = 0;
 104
 105       if (nMeasurements!=0) {
 106         dndeta = sum/Float_t(nMeasurements);
 107         error  = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
 108
 109         dndeta = sum_xw/sum_w;
 110         error  = 1/TMath::Sqrt(sum_w);
 111
 112         dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(i_eta);
 113         error  = error/hdNdEta[vertexPos]->GetBinWidth(i_eta);
 114
 115         hdNdEta[vertexPos]->SetBinContent(i_eta, dndeta);
 116         hdNdEta[vertexPos]->SetBinError(i_eta, error);
 117       }
 118     }
 119
 120   }
 121 }
 122
 123 //____________________________________________________________________
 124 void
 125 dNdEtaAnalysis::SaveHistograms() {
 126
 127   gDirectory->mkdir(GetName());
 128   gDirectory->cd(GetName());
 129
 130   hEtaVsVtx  ->Write();
 131   hVtx       ->Write();
 132   for (Int_t i=0; i<kVertexBinning; ++i)
 133     hdNdEta[i]    ->Write();
 134
 135   gDirectory->cd("../");
 136 }
 137
 138 //____________________________________________________________________
 139 void dNdEtaAnalysis::DrawHistograms()
 140 {
 141   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
 142   canvas->Divide(2, 2);
 143
 144   canvas->cd(1);
 145   if (hEtaVsVtx)
 146     hEtaVsVtx->Draw("COLZ");
 147
 148   canvas->cd(2);
 149   if (hEtaVsVtxUncorrected)
 150     hEtaVsVtxUncorrected->Draw("COLZ");
 151
 152   canvas->cd(3);
 153   if (hVtx)
 154     hVtx->Draw();
 155
 156   canvas->cd(4);
 157   if (hdNdEta[0])
 158     hdNdEta[0]->Draw();
 159 }
 160
 161 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
 162 {
 163   // Merges a list of dNdEtaAnalysis objects with this one.
 164   // This is needed for PROOF.
 165   // Returns the number of merged objects (including this)
 166
 167   if (!list)
 168     return 0;
 169
 170   if (list->IsEmpty())
 171     return 1;
 172
 173   TIterator* iter = list->MakeIterator();
 174   TObject* obj;
 175
 176   // sub collections
 177   const Int_t nCollections = kVertexBinning + 3;
 178   TList* collections[nCollections];
 179   for (Int_t i=0; i<nCollections; ++i)
 180     collections[i] = new TList;
 181
 182   Int_t count = 0;
 183   while ((obj = iter->Next()))
 184   {
 185     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
 186     if (entry == 0)
 187       continue;
 188
 189     collections[0]->Add(entry->hEtaVsVtx);
 190     collections[1]->Add(entry->hEtaVsVtxUncorrected);
 191     collections[2]->Add(entry->hVtx);
 192
 193     for (Int_t i=0; i<kVertexBinning; ++i)
 194       collections[3+i]->Add(entry->hdNdEta[i]);
 195
 196     ++count;
 197   }
 198
 199   hEtaVsVtx->Merge(collections[0]);
 200   hEtaVsVtxUncorrected->Merge(collections[1]);
 201   hVtx->Merge(collections[2]);
 202   for (Int_t i=0; i<kVertexBinning; ++i)
 203     hdNdEta[i]->Merge(collections[3+i]);
 204
 205   for (Int_t i=0; i<nCollections; ++i)
 206     delete collections[i];
 207
 208   return count+1;
 209 }