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