PWG0/dNdEta/dNdEtaAnalysis.cxx

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