updating makefile to include evgen headers to compile systematics selector
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
CommitLineData
dc740de4 1/* $Id$ */
2
75ec0f41 3#include "dNdEtaAnalysis.h"
4
ceb5d1b5 5#include <TFile.h>
45e97e28 6#include <TH3F.h>
847489f7 7#include <TH2D.h>
ceb5d1b5 8#include <TH1D.h>
9#include <TMath.h>
10#include <TCanvas.h>
7029240a 11#include <TCollection.h>
12#include <TIterator.h>
13#include <TList.h>
fcf2fb36 14#include <TLegend.h>
15
45e97e28 16#include "AlidNdEtaCorrection.h"
847489f7 17#include "AliPWG0Helper.h"
ceb5d1b5 18
75ec0f41 19//____________________________________________________________________
b7f4a1fd 20ClassImp(dNdEtaAnalysis)
75ec0f41 21
1afae8ff 22//____________________________________________________________________
23dNdEtaAnalysis::dNdEtaAnalysis() :
24 TNamed(),
25 fData(0),
26 fDataUncorrected(0),
27 fVtx(0),
28 fPtDist(0)
29{
30 // default constructor
31
32 for (Int_t i=0; i<kVertexBinning; ++i)
33 {
34 fdNdEta[i] = 0;
35 fdNdEtaPtCutOffCorrected[i] = 0;
36 }
37}
38
75ec0f41 39//____________________________________________________________________
7029240a 40dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title) :
16e24ca3 41 TNamed(name, title),
45e97e28 42 fData(0),
43 fDataUncorrected(0),
1afae8ff 44 fVtx(0),
45 fPtDist(0)
7029240a 46{
6bf0714d 47 // constructor
4dd2ad81 48
5c495d37 49 fData = new TH3F(Form("%s_analysis", name),"Input data",80,-20,20,40,-2,2,100, 0, 10);
45e97e28 50 fData->SetXTitle("vtx z [cm]");
51 fData->SetYTitle("#eta");
52 fData->SetZTitle("p_{T}");
0ab29cfa 53 fData->Sumw2();
6bf0714d 54
45e97e28 55 fDataUncorrected = dynamic_cast<TH3F*> (fData->Clone(Form("%s_analysis_uncorrected", name)));
1afae8ff 56 fDataUncorrected->SetTitle(Form("%s uncorrected", fData->GetTitle()));
45e97e28 57 fVtx = dynamic_cast<TH1D*> (fData->Project3D("x"));
1afae8ff 58 fVtx->SetName(Form("%s_vtx", name));
59 fVtx->SetTitle("Vertex distribution");
60 fVtx->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
0ab29cfa 61 //fVtx->Sumw2();
45e97e28 62
63 fdNdEta[0] = dynamic_cast<TH1D*> (fData->Project3D("y"));
1afae8ff 64 fdNdEta[0]->SetName(Form("%s_dNdEta", name));
25db2d85 65 fdNdEta[0]->SetTitle("dN_{ch}/d#eta");
1afae8ff 66 fdNdEta[0]->GetXaxis()->SetTitle(fData->GetYaxis()->GetTitle());
25db2d85 67 fdNdEta[0]->SetYTitle("dN_{ch}/d#eta");
1afae8ff 68
69 fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
70
45e97e28 71 for (Int_t i=1; i<kVertexBinning; ++i)
7029240a 72 {
45e97e28 73 fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
1afae8ff 74 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
7029240a 75 }
75ec0f41 76
1afae8ff 77 fPtDist = dynamic_cast<TH1D*> (fData->Project3D("z"));
78 fPtDist->SetName(Form("%s_pt", name));
d09fb536 79 fPtDist->SetTitle("p_{T} [GeV/c]");
1afae8ff 80 fPtDist->GetXaxis()->SetTitle(fData->GetZaxis()->GetTitle());
d09fb536 81 fPtDist->SetYTitle("#frac{dN}{d#eta dp_{T}} [c/GeV]");
75ec0f41 82}
83
16e24ca3 84//____________________________________________________________________
85dNdEtaAnalysis::~dNdEtaAnalysis()
86{
87 // destructor
88
0ab29cfa 89 if (fData)
90 {
91 delete fData;
92 fData = 0;
93 }
16e24ca3 94
0ab29cfa 95 if (fDataUncorrected)
96 {
97 delete fDataUncorrected;
98 fDataUncorrected = 0;
99 }
16e24ca3 100
0ab29cfa 101 if (fVtx)
102 {
103 delete fVtx;
104 fVtx = 0;
105 }
16e24ca3 106
107 for (Int_t i=0; i<kVertexBinning; ++i)
108 {
0ab29cfa 109 if (fdNdEta[i])
110 {
111 delete fdNdEta[i];
112 fdNdEta[i] = 0;
113 }
114 if (fdNdEtaPtCutOffCorrected[i])
115 {
116 delete fdNdEtaPtCutOffCorrected[i];
117 fdNdEtaPtCutOffCorrected[i] = 0;
118 }
16e24ca3 119 }
1afae8ff 120
0ab29cfa 121 if (fPtDist)
122 {
123 delete fPtDist;
124 fPtDist = 0;
125 }
16e24ca3 126}
127
128//_____________________________________________________________________________
129dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
130 TNamed(c),
45e97e28 131 fData(0),
132 fDataUncorrected(0),
1afae8ff 133 fVtx(0),
134 fPtDist(0)
16e24ca3 135{
136 //
137 // dNdEtaAnalysis copy constructor
138 //
139
140 ((dNdEtaAnalysis &) c).Copy(*this);
141}
142
143//_____________________________________________________________________________
144dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
145{
146 //
147 // Assignment operator
148 //
149
150 if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
151 return *this;
152}
153
154//_____________________________________________________________________________
155void dNdEtaAnalysis::Copy(TObject &c) const
156{
157 //
158 // Copy function
159 //
160
161 dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
162
45e97e28 163 target.fData = dynamic_cast<TH3F*> (fData->Clone());
164 target.fDataUncorrected = dynamic_cast<TH3F*> (fDataUncorrected->Clone());
16e24ca3 165 target.fVtx = dynamic_cast<TH1D*> (fVtx->Clone());
166
167 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 168 {
16e24ca3 169 target.fdNdEta[i] = dynamic_cast<TH1D*> (fdNdEta[i]->Clone());
1afae8ff 170 target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (fdNdEtaPtCutOffCorrected[i]->Clone());
171 }
16e24ca3 172
1afae8ff 173 target.fPtDist = dynamic_cast<TH1D*> (fPtDist->Clone());
45e97e28 174
16e24ca3 175 TNamed::Copy((TNamed &) c);
176}
177
75ec0f41 178//____________________________________________________________________
45e97e28 179void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt, Float_t weight)
6bf0714d 180{
181 // fills a track into the histograms
182
45e97e28 183 fDataUncorrected->Fill(vtx, eta, pt);
184 fData->Fill(vtx, eta, pt, weight);
75ec0f41 185}
186
187//____________________________________________________________________
45e97e28 188void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t weight)
6bf0714d 189{
190 // fills an event into the histograms
191
847489f7 192 fVtx->Fill(vtx, weight);
75ec0f41 193}
194
195//____________________________________________________________________
847489f7 196void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut)
fcf2fb36 197{
198 // correct with correction values if available
6bf0714d 199
fcf2fb36 200 // TODO what do we do with the error?
201 if (!correction)
202 printf("INFO: No correction applied\n");
203
847489f7 204 // In fData we have the track2particle and vertex reconstruction efficiency correction already applied
fcf2fb36 205
1afae8ff 206 // create pt hist
207 {
208 // reset all ranges
209 fData->GetXaxis()->SetRange(0, 0);
210 fData->GetYaxis()->SetRange(0, 0);
211 fData->GetZaxis()->SetRange(0, 0);
212
213 // vtx cut
214 Int_t vertexBinBegin = fData->GetXaxis()->FindBin(-5);
215 Int_t vertexBinEnd = fData->GetXaxis()->FindBin(5);
216 fData->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
217 Float_t nEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd);
218
219 // eta cut
220 fData->GetYaxis()->SetRange(fData->GetYaxis()->FindBin(-0.8), fData->GetYaxis()->FindBin(0.8));
221 Float_t etaWidth = 1.6;
222
223 TH1D* ptHist = dynamic_cast<TH1D*> (fData->Project3D("ze"));
224
225 for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
226 {
227 Float_t binSize = fPtDist->GetBinWidth(i);
228 fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
229 fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
230 }
231
232 delete ptHist;
233 }
234
235 // reset all ranges
236 fData->GetXaxis()->SetRange(0, 0);
237 fData->GetYaxis()->SetRange(0, 0);
238 fData->GetZaxis()->SetRange(0, 0);
239
847489f7 240 // integrate over pt (with pt cut)
1afae8ff 241 Int_t ptLowBin = 1;
242 if (ptCut > 0)
243 ptLowBin = fData->GetZaxis()->FindBin(ptCut);
244
245 fData->GetZaxis()->SetRange(ptLowBin, fData->GetZaxis()->GetNbins());
246 TH2D* vtxVsEta = dynamic_cast<TH2D*> (fData->Project3D("yx2e"));
247 vtxVsEta->GetXaxis()->SetTitle(fData->GetXaxis()->GetTitle());
248 vtxVsEta->GetYaxis()->SetTitle(fData->GetYaxis()->GetTitle());
249
847489f7 250 if (vtxVsEta == 0)
251 {
252 printf("ERROR: pt integration failed\n");
253 return;
fcf2fb36 254 }
255
d09fb536 256 //new TCanvas;
257 //vtxVsEta->Draw("COLZ");
847489f7 258
259 for (Int_t iEta=0; iEta<=vtxVsEta->GetNbinsY(); iEta++)
5af55649 260 {
7029240a 261 // do we have several histograms for different vertex positions?
6bf0714d 262 Int_t vertexBinWidth = fVtx->GetNbinsX() / (kVertexBinning-1);
7029240a 263 for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
264 {
5af55649 265 Int_t vertexBinBegin = 1;
6bf0714d 266 Int_t vertexBinEnd = fVtx->GetNbinsX() + 1;
fcf2fb36 267
268 // the first histogram is always for the whole vertex range
269 if (vertexPos > 0)
270 {
5af55649 271 vertexBinBegin = 1 + vertexBinWidth * (vertexPos-1);
fcf2fb36 272 vertexBinEnd = vertexBinBegin + vertexBinWidth;
273 }
7029240a 274
6bf0714d 275 Float_t totalEvents = fVtx->Integral(vertexBinBegin, vertexBinEnd - 1);
5af55649 276 if (totalEvents == 0)
277 {
278 printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
279 continue;
7029240a 280 }
7029240a 281
5af55649 282 Float_t sum = 0;
283 Float_t sumError2 = 0;
284 for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
285 {
847489f7 286 if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
5af55649 287 {
847489f7 288 sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
289 sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
5af55649 290 }
291 }
7029240a 292
1afae8ff 293 Float_t ptCutOffCorrection = 1;
294 if (correction)
295 ptCutOffCorrection = correction->GetMeasuredFraction(ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
296
847489f7 297 if (ptCutOffCorrection <= 0)
298 {
299 printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
300 continue;
301 }
302
1afae8ff 303 Float_t dndeta = sum / totalEvents;
304 Float_t error = TMath::Sqrt(sumError2) / totalEvents;
7029240a 305
6bf0714d 306 dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
307 error = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
7029240a 308
6bf0714d 309 fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
310 fdNdEta[vertexPos]->SetBinError(iEta, error);
1afae8ff 311
312 dndeta /= ptCutOffCorrection;
313 error /= ptCutOffCorrection;
314
315 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(iEta, dndeta);
316 fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(iEta, error);
75ec0f41 317 }
847489f7 318 }
75ec0f41 319}
320
75ec0f41 321//____________________________________________________________________
6bf0714d 322void dNdEtaAnalysis::SaveHistograms()
323{
324 // save the histograms to a directory with the name of this class (retrieved from TNamed)
75ec0f41 325
7029240a 326 gDirectory->mkdir(GetName());
327 gDirectory->cd(GetName());
5fbd0b17 328
1afae8ff 329 if (fData)
330 {
331 fData->Write();
332 AliPWG0Helper::CreateProjections(fData);
333 }
334 else
335 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
336
337 if (fDataUncorrected)
338 {
339 fDataUncorrected->Write();
340 AliPWG0Helper::CreateProjections(fDataUncorrected);
341 }
342 else
343 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fDataUncorrected is 0\n");
344
d09fb536 345 if (fVtx)
1afae8ff 346 fVtx ->Write();
347 else
348 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fVtx is 0\n");
847489f7 349
d09fb536 350 if (fPtDist)
351 fPtDist ->Write();
352 else
353 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
354
7029240a 355 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 356 {
357 if (fdNdEta[i])
358 fdNdEta[i]->Write();
359 else
360 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
361
362 if (fdNdEtaPtCutOffCorrected[i])
363 fdNdEtaPtCutOffCorrected[i]->Write();
364 else
365 printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
366 }
75ec0f41 367
368 gDirectory->cd("../");
369}
370
5fbd0b17 371void dNdEtaAnalysis::LoadHistograms()
372{
6bf0714d 373 // loads the histograms from a directory with the name of this class (retrieved from TNamed)
374
5fbd0b17 375 gDirectory->cd(GetName());
376
45e97e28 377 fData = dynamic_cast<TH3F*> (gDirectory->Get(fData->GetName()));
378 fDataUncorrected = dynamic_cast<TH3F*> (gDirectory->Get(fDataUncorrected->GetName()));
5fbd0b17 379
6bf0714d 380 fVtx = dynamic_cast<TH1D*> (gDirectory->Get(fVtx->GetName()));
5fbd0b17 381
382 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 383 {
6bf0714d 384 fdNdEta[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEta[i]->GetName()));
1afae8ff 385 fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1D*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
386 }
5fbd0b17 387
0ab29cfa 388 fPtDist = dynamic_cast<TH1D*> (gDirectory->Get(fPtDist->GetName()));
389
5fbd0b17 390 gDirectory->cd("../");
391}
392
ceb5d1b5 393//____________________________________________________________________
394void dNdEtaAnalysis::DrawHistograms()
395{
6bf0714d 396 // draws the histograms
397
5fbd0b17 398 TCanvas* canvas = new TCanvas("dNdEtaAnalysis", "dNdEtaAnalysis", 800, 800);
399 canvas->Divide(2, 2);
ceb5d1b5 400
401 canvas->cd(1);
45e97e28 402 if (fData)
403 fData->Draw("COLZ");
ceb5d1b5 404
1afae8ff 405 /*canvas->cd(2);
45e97e28 406 if (fDataUncorrected)
1afae8ff 407 fDataUncorrected->Draw("COLZ");*/
ceb5d1b5 408
1afae8ff 409 canvas->cd(2);
6bf0714d 410 if (fVtx)
411 fVtx->Draw();
ceb5d1b5 412
1afae8ff 413 canvas->cd(3);
414 if (fdNdEtaPtCutOffCorrected[0])
415 fdNdEtaPtCutOffCorrected[0]->Draw();
416
6bf0714d 417 if (fdNdEta[0])
1afae8ff 418 {
419 fdNdEta[0]->SetLineColor(kRed);
420 fdNdEta[0]->Draw("SAME");
421 }
422
423 canvas->cd(4);
424 if (fPtDist)
425 fPtDist->Draw();
fcf2fb36 426
427 // histograms for different vertices?
428 if (kVertexBinning > 0)
429 {
430 // doesnt work, but i dont get it, giving up...
431 /*TCanvas* canvas2 =*/ new TCanvas("dNdEtaAnalysisVtx", "dNdEtaAnalysisVtx", 450, 450);
432 //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
433 //printf("%d\n", yPads);
434 //canvas2->Divide(2, yPads);
435
436 TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
437
5af55649 438 for (Int_t i=0; i<kVertexBinning; ++i)
fcf2fb36 439 {
440 //canvas2->cd(i-1);
441 //printf("%d\n", i);
1afae8ff 442 if (fdNdEtaPtCutOffCorrected[i])
fcf2fb36 443 {
1afae8ff 444 fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
445 fdNdEtaPtCutOffCorrected[i]->Draw((i == 0) ? "" : "SAME");
446 legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
fcf2fb36 447 }
448 }
449
450 legend->Draw();
451 }
7029240a 452}
453
454Long64_t dNdEtaAnalysis::Merge(TCollection* list)
455{
456 // Merges a list of dNdEtaAnalysis objects with this one.
457 // This is needed for PROOF.
458 // Returns the number of merged objects (including this)
459
460 if (!list)
461 return 0;
462
463 if (list->IsEmpty())
464 return 1;
465
466 TIterator* iter = list->MakeIterator();
467 TObject* obj;
468
469 // sub collections
1afae8ff 470 const Int_t nCollections = 2 * kVertexBinning + 4; // 4 standalone hists, two arrays of size kVertexBinning
7029240a 471 TList* collections[nCollections];
472 for (Int_t i=0; i<nCollections; ++i)
473 collections[i] = new TList;
474
475 Int_t count = 0;
476 while ((obj = iter->Next()))
477 {
478 dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
479 if (entry == 0)
480 continue;
481
45e97e28 482 collections[0]->Add(entry->fData);
0ab29cfa 483 if (fDataUncorrected)
484 collections[1]->Add(entry->fDataUncorrected);
6bf0714d 485 collections[2]->Add(entry->fVtx);
1afae8ff 486 collections[3]->Add(entry->fPtDist);
7029240a 487
488 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 489 {
490 collections[4+i]->Add(entry->fdNdEta[i]);
491 collections[4+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
492 }
7029240a 493
494 ++count;
495 }
496
45e97e28 497 fData->Merge(collections[0]);
0ab29cfa 498 if (fDataUncorrected)
499 fDataUncorrected->Merge(collections[1]);
6bf0714d 500 fVtx->Merge(collections[2]);
1afae8ff 501 fPtDist->Merge(collections[3]);
7029240a 502 for (Int_t i=0; i<kVertexBinning; ++i)
1afae8ff 503 {
504 fdNdEta[i]->Merge(collections[4+i]);
505 fdNdEta[i]->Merge(collections[4+kVertexBinning+i]);
506 }
7029240a 507
508 for (Int_t i=0; i<nCollections; ++i)
509 delete collections[i];
510
511 return count+1;
ceb5d1b5 512}