changing from printf to AliDebug
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
CommitLineData
75ec0f41 1#include "dNdEtaAnalysis.h"
2
ceb5d1b5 3#include <TFile.h>
4#include <TH2F.h>
5#include <TH1D.h>
6#include <TMath.h>
7#include <TCanvas.h>
8
75ec0f41 9//____________________________________________________________________
b7f4a1fd 10ClassImp(dNdEtaAnalysis)
75ec0f41 11
12//____________________________________________________________________
13dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name) {
14
15 fName = TString(name);
16
ceb5d1b5 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");
75ec0f41 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();
75ec0f41 30}
31
32//____________________________________________________________________
33void
34dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t weight) {
35 hEtaVsVtx->Fill(vtx, eta, weight);
ceb5d1b5 36 hEtaVsVtxUncorrected->Fill(vtx,eta);
75ec0f41 37}
38
39//____________________________________________________________________
40void
41dNdEtaAnalysis::FillEvent(Float_t vtx) {
42 hVtx->Fill(vtx);
43}
44
45//____________________________________________________________________
46void
47dNdEtaAnalysis::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
ceb5d1b5 103 hdNdEta->SetBinContent(i_eta, dndeta);
104 hdNdEta->SetBinError(i_eta, error);
75ec0f41 105 }
106
107 }
108}
109
110
111
112//____________________________________________________________________
113void
114dNdEtaAnalysis::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
ceb5d1b5 126//____________________________________________________________________
127void 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}