o) splitted analysis and correction of dneta
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
1 /* $Id$ */
2
3 #include "dNdEtaAnalysis.h"
4
5 #include <TFile.h>
6 #include <TH2F.h>
7 #include <TH1D.h>
8 #include <TMath.h>
9 #include <TCanvas.h>
10 #include <TCollection.h>
11 #include <TIterator.h>
12 #include <TList.h>
13 #include <TLegend.h>
14
15 #include "dNdEtaCorrection.h"
16
17 //____________________________________________________________________
18 ClassImp(dNdEtaAnalysis)
19
20 //____________________________________________________________________
21 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
22 TNamed(name, title)
23 {
24   hEtaVsVtx  = new TH2F(Form("%s_eta_vs_vtx", name),"",80,-20,20,120,-6,6);
25   hEtaVsVtx->SetXTitle("vtx z [cm]");
26   hEtaVsVtx->SetYTitle("#eta");
27
28   hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (hEtaVsVtx->Clone(Form("%s_eta_vs_vtx_uncorrected", name)));
29   hVtx       = hEtaVsVtx->ProjectionX(Form("%s_vtx", name));
30   for (Int_t i=0; i<kVertexBinning; ++i)
31   {
32     hdNdEta[i]    = hEtaVsVtx->ProjectionY(Form("%s_dNdEta_%d", name, i));
33     hdNdEta[i]->SetYTitle("dN/d#eta");
34   }
35
36   hEtaVsVtx->Sumw2();
37   hVtx->Sumw2();
38 }
39
40 //____________________________________________________________________
41 void
42 dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta) {
43   hEtaVsVtxUncorrected->Fill(vtx,eta);
44 }
45
46 //____________________________________________________________________
47 void
48 dNdEtaAnalysis::FillEvent(Float_t vtx) {
49   hVtx->Fill(vtx);
50 }
51
52 //____________________________________________________________________
53 void dNdEtaAnalysis::Finish(dNdEtaCorrection* correction)
54 {
55   // correct with correction values if available
56   // TODO what do we do with the error?
57   if (!correction)
58     printf("INFO: No correction applied\n");
59
60   // this can be replaced by TH2F::Divide if we agree that the binning will be always the same
61   for (Int_t iVtx=0; iVtx<=hEtaVsVtxUncorrected->GetNbinsX(); iVtx++)
62   {
63     for (Int_t iEta=0; iEta<=hEtaVsVtxUncorrected->GetNbinsY(); iEta++)
64     {
65       Float_t correctionValue = 1;
66       if (correction)
67         correctionValue = correction->GetCorrection(hEtaVsVtxUncorrected->GetXaxis()->GetBinCenter(iVtx), hEtaVsVtxUncorrected->GetYaxis()->GetBinCenter(iEta));
68
69       Float_t value = hEtaVsVtxUncorrected->GetBinContent(iVtx, iEta);
70       Float_t error = hEtaVsVtxUncorrected->GetBinError(iVtx, iEta);
71
72       Float_t correctedValue = value * correctionValue;
73       Float_t correctedError = error * correctionValue;
74
75       if (correctedValue != 0)
76       {
77         hEtaVsVtx->SetBinContent(iVtx, iEta, correctedValue);
78         hEtaVsVtx->SetBinError(iVtx, iEta, correctedError);
79       }
80     }
81   }
82
83   for (Int_t iEta=0; iEta<=hEtaVsVtx->GetNbinsY(); iEta++)
84   {
85     // do we have several histograms for different vertex positions?
86     Int_t vertexBinWidth = hVtx->GetNbinsX() / (kVertexBinning-1);
87     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
88     {
89       Int_t vertexBinBegin = 1;
90       Int_t vertexBinEnd = hVtx->GetNbinsX() + 1;
91
92       // the first histogram is always for the whole vertex range
93       if (vertexPos > 0)
94       {
95         vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
96         vertexBinEnd = vertexBinBegin + vertexBinWidth;
97       }
98
99       Float_t totalEvents = hVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
100       if (totalEvents == 0)
101       {
102         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
103         continue;
104       }
105
106       Float_t sum = 0;
107       Float_t sumError2 = 0;
108       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
109       {
110         if (hEtaVsVtx->GetBinContent(iVtx, iEta) != 0)
111         {
112           sum = sum + hEtaVsVtx->GetBinContent(iVtx, iEta);
113           sumError2 = sumError2 + TMath::Power(hEtaVsVtx->GetBinError(iVtx, iEta),2);
114         }
115       }
116
117       Float_t dndeta = sum / totalEvents;
118       Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
119
120       dndeta = dndeta/hdNdEta[vertexPos]->GetBinWidth(iEta);
121       error  = error/hdNdEta[vertexPos]->GetBinWidth(iEta);
122
123       hdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
124       hdNdEta[vertexPos]->SetBinError(iEta, error);
125     }
126   }
127 }
128
129 //____________________________________________________________________
130 void
131 dNdEtaAnalysis::SaveHistograms() {
132
133   gDirectory->mkdir(GetName());
134   gDirectory->cd(GetName());
135
136   hEtaVsVtx  ->Write();
137   hEtaVsVtxUncorrected->Write();
138   hVtx       ->Write();
139   for (Int_t i=0; i<kVertexBinning; ++i)
140     hdNdEta[i]    ->Write();
141
142   gDirectory->cd("../");
143 }
144
145 void dNdEtaAnalysis::LoadHistograms()
146 {
147   gDirectory->cd(GetName());
148
149   hEtaVsVtx = dynamic_cast<TH2F*> (gDirectory->Get(hEtaVsVtx->GetName()));
150   hEtaVsVtxUncorrected = dynamic_cast<TH2F*> (gDirectory->Get(hEtaVsVtxUncorrected->GetName()));
151
152   hVtx = dynamic_cast<TH1D*> (gDirectory->Get(hVtx->GetName()));
153
154   for (Int_t i=0; i<kVertexBinning; ++i)
155     hdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(hdNdEta[i]->GetName()));
156
157   gDirectory->cd("../");
158 }
159
160 //____________________________________________________________________
161 void dNdEtaAnalysis::DrawHistograms()
162 {
163   TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
164   canvas->Divide(2, 2);
165
166   canvas->cd(1);
167   if (hEtaVsVtx)
168     hEtaVsVtx->Draw("COLZ");
169
170   canvas->cd(2);
171   if (hEtaVsVtxUncorrected)
172     hEtaVsVtxUncorrected->Draw("COLZ");
173
174   canvas->cd(3);
175   if (hVtx)
176     hVtx->Draw();
177
178   canvas->cd(4);
179   if (hdNdEta[0])
180     hdNdEta[0]->Draw();
181
182     // histograms for different vertices?
183   if (kVertexBinning > 0)
184   {
185       // doesnt work, but i dont get it, giving up...
186     /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
187     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
188     //printf("%d\n", yPads);
189     //canvas2->Divide(2, yPads);
190
191     TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
192
193     for (Int_t i=0; i<kVertexBinning; ++i)
194     {
195       //canvas2->cd(i-1);
196       //printf("%d\n", i);
197       if (hdNdEta[i])
198       {
199         hdNdEta[i]->SetLineColor(i+1);
200         hdNdEta[i]->Draw((i == 0) ? "" : "SAME");
201         legend->AddEntry(hdNdEta[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
202       }
203     }
204
205     legend->Draw();
206   }
207 }
208
209 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
210 {
211   // Merges a list of dNdEtaAnalysis objects with this one.
212   // This is needed for PROOF.
213   // Returns the number of merged objects (including this)
214
215   if (!list)
216     return 0;
217
218   if (list->IsEmpty())
219     return 1;
220
221   TIterator* iter = list->MakeIterator();
222   TObject* obj;
223
224   // sub collections
225   const Int_t nCollections = kVertexBinning + 3;
226   TList* collections[nCollections];
227   for (Int_t i=0; i<nCollections; ++i)
228     collections[i] = new TList;
229
230   Int_t count = 0;
231   while ((obj = iter->Next()))
232   {
233     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
234     if (entry == 0)
235       continue;
236
237     collections[0]->Add(entry->hEtaVsVtx);
238     collections[1]->Add(entry->hEtaVsVtxUncorrected);
239     collections[2]->Add(entry->hVtx);
240
241     for (Int_t i=0; i<kVertexBinning; ++i)
242       collections[3+i]->Add(entry->hdNdEta[i]);
243
244     ++count;
245   }
246
247   hEtaVsVtx->Merge(collections[0]);
248   hEtaVsVtxUncorrected->Merge(collections[1]);
249   hVtx->Merge(collections[2]);
250   for (Int_t i=0; i<kVertexBinning; ++i)
251     hdNdEta[i]->Merge(collections[3+i]);
252
253   for (Int_t i=0; i<nCollections; ++i)
254     delete collections[i];
255
256   return count+1;
257 }