Removed ; after ClassImp(dNdEtaAnalysis)
[u/mrichter/AliRoot.git] / 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