first checkin of Claus' code with a few modifications
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
CommitLineData
75ec0f41 1#include "dNdEtaAnalysis.h"
2
3//____________________________________________________________________
4ClassImp(dNdEtaAnalysis);
5
6//____________________________________________________________________
7dNdEtaAnalysis::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//____________________________________________________________________
27void
28dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) {
29 hEtaVsVtx->Fill(vtx, eta, weight);
30}
31
32//____________________________________________________________________
33void
34dNdEtaAnalysis::FillEvent(Float_t vtx) {
35 hVtx->Fill(vtx);
36}
37
38//____________________________________________________________________
39void
40dNdEtaAnalysis::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//____________________________________________________________________
106void
107dNdEtaAnalysis::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