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