PWG0/dNdEta/dNdEtaAnalysis.cxx

   1 #include "dNdEtaAnalysis.h"
   2
   3 #include <TFile.h>
   4 #include <TH2F.h>
   5 #include <TH1D.h>
   6 #include <TMath.h>
   7 #include <TCanvas.h>
   8
   9 //____________________________________________________________________
  10 ClassImp(dNdEtaAnalysis)
  11
  12 //____________________________________________________________________
  13 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name) {
  14
  15   fName = TString(name);
  16
  17   hEtaVsVtx  = new TH2F("eta_vs_vtx","",80,-20,20,120,-6,6);
  18   hEtaVsVtx->SetXTitle("vtx z [cm]");
  19   hEtaVsVtx->SetYTitle("#eta");
  20
  21   hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone("eta_vs_vtx_uncorrected"));
  22   hVtx       = hEtaVsVtx->ProjectionX("vtx");
  23   hdNdEta    = hEtaVsVtx->ProjectionY("dNdEta");
  24
  25   hdNdEta->SetYTitle("dN/d#eta");
  26
  27   hEtaVsVtx->Sumw2();
  28   hVtx->Sumw2();
  29 }
  30
  31 //____________________________________________________________________
  32 void
  33 dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) {
  34   hEtaVsVtx->Fill(vtx, eta, weight);
  35   hEtaVsVtxUncorrected->Fill(vtx,eta);
  36 }
  37
  38 //____________________________________________________________________
  39 void
  40 dNdEtaAnalysis::FillEvent(Float_t vtx) {
  41   hVtx->Fill(vtx);
  42 }
  43
  44 //____________________________________________________________________
  45 void
  46 dNdEtaAnalysis::Finish() {
  47
  48   // first normalize with n events (per vtx)
  49   for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) {
  50     Float_t nEvents      = hVtx->GetBinContent(i_vtx);
  51     Float_t nEventsError = hVtx->GetBinError(i_vtx);
  52
  53     if (nEvents==0) continue;
  54
  55     for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
  56       Float_t value = hEtaVsVtx->GetBinContent(i_vtx, i_eta)/nEvents;
  57       if (value==0) continue;
  58       Float_t error = hEtaVsVtx->GetBinError(i_vtx, i_eta)/nEvents;
  59       error = TMath::Sqrt(TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta)/
  60                                        hEtaVsVtx->GetBinContent(i_vtx, i_eta),2) +
  61                           TMath::Power(nEventsError/nEvents,2));
  62       hEtaVsVtx->SetBinContent(i_vtx, i_eta, value);
  63       hEtaVsVtx->SetBinError(i_vtx,   i_eta, error);
  64     }
  65   }
  66
  67   // then take the wieghted average for each eta
  68   // is this the right way to do it???
  69   for (Int_t i_eta=0; i_eta<=hEtaVsVtx->GetNbinsY(); i_eta++) {
  70     Float_t sum           = 0;
  71     Float_t sumError2     = 0;
  72     Int_t   nMeasurements = 0;
  73
  74     Float_t sum_xw = 0;
  75     Float_t sum_w  = 0;
  76
  77     for (Int_t i_vtx=0; i_vtx<=hVtx->GetNbinsX(); i_vtx++) {
  78       if (hVtx->GetBinContent(i_vtx)==0)             continue;
  79       if (hEtaVsVtx->GetBinContent(i_vtx, i_eta)==0) continue;
  80
  81       Float_t w = 1/TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
  82       sum_xw = sum_xw + hEtaVsVtx->GetBinContent(i_vtx, i_eta)*w;
  83       sum_w  = sum_w + w;
  84
  85       sum = sum + hEtaVsVtx->GetBinContent(i_vtx, i_eta);
  86       sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(i_vtx, i_eta),2);
  87       nMeasurements++;
  88     }
  89     Float_t dndeta = 0;
  90     Float_t error  = 0;
  91
  92     if (nMeasurements!=0) {
  93       dndeta = sum/Float_t(nMeasurements);
  94       error  = TMath::Sqrt(sumError2)/Float_t(nMeasurements);
  95
  96       dndeta = sum_xw/sum_w;
  97       error  = 1/TMath::Sqrt(sum_w);
  98
  99       dndeta = dndeta/hdNdEta->GetBinWidth(i_eta);
 100       error  = error/hdNdEta->GetBinWidth(i_eta);
 101
 102       hdNdEta->SetBinContent(i_eta, dndeta);
 103       hdNdEta->SetBinError(i_eta, error);
 104     }
 105
 106   }
 107 }
 108
 109
 110
 111 //____________________________________________________________________
 112 void
 113 dNdEtaAnalysis::SaveHistograms() {
 114
 115   gDirectory->mkdir(fName.Data());
 116   gDirectory->cd(fName.Data());
 117
 118   hEtaVsVtx  ->Write();
 119   hVtx       ->Write();
 120   hdNdEta    ->Write();
 121
 122   gDirectory->cd("../");
 123 }
 124
 125 //____________________________________________________________________
 126 void dNdEtaAnalysis::DrawHistograms()
 127 {
 128   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
 129   canvas->Divide(2, 2);
 130
 131   canvas->cd(1);
 132   if (hEtaVsVtx)
 133     hEtaVsVtx->Draw("COLZ");
 134
 135   canvas->cd(2);
 136   if (hEtaVsVtxUncorrected)
 137     hEtaVsVtxUncorrected->Draw("COLZ");
 138
 139   canvas->cd(3);
 140   if (hVtx)
 141     hVtx->Draw();
 142
 143   canvas->cd(4);
 144   if (hdNdEta)
 145     hdNdEta->Draw();
 146 }