b4341dd7440c9f1f5b08c61afb12fc1eb46d2d1b
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / dNdEtaAnalysis.cxx
1 /* $Id$ */
2
3 #include "dNdEtaAnalysis.h"
4
5 #include <TFile.h>
6 #include <TH3F.h>
7 #include <TH2F.h>
8 #include <TH1F.h>
9 #include <TMath.h>
10 #include <TCanvas.h>
11 #include <TCollection.h>
12 #include <TIterator.h>
13 #include <TList.h>
14 #include <TLegend.h>
15 #include <TLine.h>
16
17 #include "AlidNdEtaCorrection.h"
18 #include <AliCorrection.h>
19 #include <AliPWG0Helper.h>
20 #include <AliCorrectionMatrix2D.h>
21 #include <AliCorrectionMatrix3D.h>
22
23 //____________________________________________________________________
24 ClassImp(dNdEtaAnalysis)
25
26 //____________________________________________________________________
27 dNdEtaAnalysis::dNdEtaAnalysis() :
28   TNamed(),
29   fData(0),
30   fPtDist(0)
31 {
32   // default constructor
33
34   for (Int_t i=0; i<kVertexBinning; ++i)
35   {
36     fdNdEta[i] = 0;
37     fdNdEtaPtCutOffCorrected[i] = 0;
38   }
39 }
40
41 //____________________________________________________________________
42 dNdEtaAnalysis::dNdEtaAnalysis(Char_t* name, Char_t* title, const char* analysis) :
43   TNamed(name, title),
44   fData(0),
45   fPtDist(0)
46 {
47   // constructor
48
49   // TODO this binning has to be the same than in AliCorrection, somehow passed?!
50   Float_t binLimitsPt[] = {0.0, 0.05, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4, 0.425, 0.45, 0.475, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
51
52   fData = new AliCorrection("Analysis", Form("%s Analysis", title), analysis);
53
54   // do not add this hists to the directory
55   Bool_t oldStatus = TH1::AddDirectoryStatus();
56   TH1::AddDirectory(kFALSE);
57
58   fdNdEta[0] = new TH1F("dNdEta", "dN_{ch}/d#eta;#eta;dN_{ch}/d#eta", 40, -2, 2);
59
60   fdNdEtaPtCutOffCorrected[0] = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_corrected", fdNdEta[0]->GetName())));
61
62   for (Int_t i=1; i<kVertexBinning; ++i)
63   {
64     fdNdEta[i]    = dynamic_cast<TH1F*> (fdNdEta[0]->Clone(Form("%s_%d", fdNdEta[0]->GetName(), i)));
65     fdNdEtaPtCutOffCorrected[i]    = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[0]->Clone(Form("%s_%d", fdNdEtaPtCutOffCorrected[0]->GetName(), i)));
66   }
67
68   fPtDist = new TH1F("Pt", "p_{T} distribution;p_{T} [GeV/c];#frac{dN}{d#eta dp_{T}} [c/GeV]", 28, binLimitsPt);
69
70   TH1::AddDirectory(oldStatus);
71 }
72
73 //____________________________________________________________________
74 dNdEtaAnalysis::~dNdEtaAnalysis()
75 {
76   // destructor
77
78   if (fData)
79   {
80     delete fData;
81     fData = 0;
82   }
83
84   for (Int_t i=0; i<kVertexBinning; ++i)
85   {
86     if (fdNdEta[i])
87     {
88       delete fdNdEta[i];
89       fdNdEta[i] = 0;
90     }
91     if (fdNdEtaPtCutOffCorrected[i])
92     {
93       delete fdNdEtaPtCutOffCorrected[i];
94       fdNdEtaPtCutOffCorrected[i] = 0;
95     }
96   }
97
98   if (fPtDist)
99   {
100     delete fPtDist;
101     fPtDist = 0;
102   }
103 }
104
105 //_____________________________________________________________________________
106 dNdEtaAnalysis::dNdEtaAnalysis(const dNdEtaAnalysis &c) :
107   TNamed(c),
108   fData(0),
109   fPtDist(0)
110 {
111   //
112   // dNdEtaAnalysis copy constructor
113   //
114
115   ((dNdEtaAnalysis &) c).Copy(*this);
116 }
117
118 //_____________________________________________________________________________
119 dNdEtaAnalysis &dNdEtaAnalysis::operator=(const dNdEtaAnalysis &c)
120 {
121   //
122   // Assignment operator
123   //
124
125   if (this != &c) ((dNdEtaAnalysis &) c).Copy(*this);
126   return *this;
127 }
128
129 //_____________________________________________________________________________
130 void dNdEtaAnalysis::Copy(TObject &c) const
131 {
132   //
133   // Copy function
134   //
135
136   dNdEtaAnalysis& target = (dNdEtaAnalysis &) c;
137
138   target.fData = dynamic_cast<AliCorrection*> (fData->Clone());
139
140   for (Int_t i=0; i<kVertexBinning; ++i)
141   {
142     target.fdNdEta[i] = dynamic_cast<TH1F*> (fdNdEta[i]->Clone());
143     target.fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (fdNdEtaPtCutOffCorrected[i]->Clone());
144   }
145
146   target.fPtDist = dynamic_cast<TH1F*> (fPtDist->Clone());
147
148   TNamed::Copy((TNamed &) c);
149 }
150
151 //____________________________________________________________________
152 void dNdEtaAnalysis::FillTrack(Float_t vtx, Float_t eta, Float_t pt)
153 {
154   // fills a track into the histograms
155
156   fData->GetTrackCorrection()->FillMeas(vtx, eta, pt);
157 }
158
159 //____________________________________________________________________
160 void dNdEtaAnalysis::FillEvent(Float_t vtx, Float_t n)
161 {
162   // fills an event into the histograms
163
164   fData->GetEventCorrection()->FillMeas(vtx, n);
165 }
166
167 //____________________________________________________________________
168 void dNdEtaAnalysis::Finish(AlidNdEtaCorrection* correction, Float_t ptCut, AlidNdEtaCorrection::CorrectionType correctionType)
169 {
170   //
171   // correct with the given correction values and calculate dNdEta and pT distribution
172   // the corrections that are applied can be steered by the flag correctionType
173   //
174
175   // TODO put tag somewhere which corrections have been applied
176
177   // set corrections to 1
178   fData->SetCorrectionToUnity();
179
180   if (correction && correctionType != AlidNdEtaCorrection::kNone)
181   {
182     TH3F* trackCorr = fData->GetTrackCorrection()->GetCorrectionHistogram();
183     TH2F* eventCorr = fData->GetEventCorrection()->GetCorrectionHistogram();
184
185     if (correctionType >= AlidNdEtaCorrection::kTrack2Particle)
186       trackCorr->Multiply(correction->GetTrack2ParticleCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
187
188     /*if (correctionType >= AlidNdEtaCorrection::kVertexReco)
189     {
190       trackCorr->Multiply(correction->GetVertexRecoCorrection()->GetTrackCorrection()->GetCorrectionHistogram());
191       eventCorr->Multiply(correction->GetVertexRecoCorrection()->GetEventCorrection()->GetCorrectionHistogram());
192     }*/
193
194     switch (correctionType)
195     {
196       case AlidNdEtaCorrection::kINEL :
197       {
198         trackCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetTrackCorrection()->GetCorrectionHistogram());
199         eventCorr->Multiply(correction->GetTriggerBiasCorrectionINEL()->GetEventCorrection()->GetCorrectionHistogram());
200         break;
201       }
202       case AlidNdEtaCorrection::kNSD :
203       {
204         trackCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetTrackCorrection()->GetCorrectionHistogram());
205         eventCorr->Multiply(correction->GetTriggerBiasCorrectionNSD()->GetEventCorrection()->GetCorrectionHistogram());
206         break;
207       }
208       case AlidNdEtaCorrection::kND :
209       {
210         trackCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetTrackCorrection()->GetCorrectionHistogram());
211         eventCorr->Multiply(correction->GetTriggerBiasCorrectionND()->GetEventCorrection()->GetCorrectionHistogram());
212         break;
213       }
214       default : break;
215     }
216   }
217   else
218     printf("INFO: No correction applied\n");
219
220   fData->Multiply();
221
222   TH3F* dataHist = fData->GetTrackCorrection()->GetGeneratedHistogram();
223
224   // integrate multiplicity axis out (include under/overflow bins!!!)
225   TH2F* tmp = fData->GetEventCorrection()->GetGeneratedHistogram();
226   TH1D* vertexHist = tmp->ProjectionX("_px", 0, tmp->GetNbinsY() + 1, "e");
227
228   // create pt hist
229   {
230     // reset all ranges
231     dataHist->GetXaxis()->SetRange(0, 0);
232     dataHist->GetYaxis()->SetRange(0, 0);
233     dataHist->GetZaxis()->SetRange(0, 0);
234
235     // vtx cut
236     Int_t vertexBinBegin = dataHist->GetXaxis()->FindBin(-5);
237     Int_t vertexBinEnd = dataHist->GetXaxis()->FindBin(5);
238     dataHist->GetXaxis()->SetRange(vertexBinBegin, vertexBinEnd);
239     Float_t nEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd);
240
241     if (nEvents > 0)
242     {
243       // eta cut
244       dataHist->GetYaxis()->SetRange(dataHist->GetYaxis()->FindBin(-0.8), dataHist->GetYaxis()->FindBin(0.8));
245       Float_t etaWidth = 1.6;
246
247       TH1D* ptHist = dynamic_cast<TH1D*> (dataHist->Project3D("ze"));
248
249       for (Int_t i=1; i<=fPtDist->GetNbinsX(); ++i)
250       {
251         Float_t binSize = fPtDist->GetBinWidth(i);
252         fPtDist->SetBinContent(i, ptHist->GetBinContent(i) / binSize / nEvents / etaWidth);
253         fPtDist->SetBinError(i, ptHist->GetBinError(i) / binSize / nEvents / etaWidth);
254       }
255
256       delete ptHist;
257     }
258     else
259       printf("ERROR: nEvents is 0!\n");
260   }
261
262   // reset all ranges
263   dataHist->GetXaxis()->SetRange(0, 0);
264   dataHist->GetYaxis()->SetRange(0, 0);
265   dataHist->GetZaxis()->SetRange(0, 0);
266
267   // integrate over pt (with pt cut)
268   Int_t ptLowBin = 1;
269   if (ptCut > 0)
270     ptLowBin = dataHist->GetZaxis()->FindBin(ptCut);
271
272   dataHist->GetZaxis()->SetRange(ptLowBin, dataHist->GetZaxis()->GetNbins());
273   printf("range %d %d\n", ptLowBin, dataHist->GetZaxis()->GetNbins());
274   TH2D* vtxVsEta = dynamic_cast<TH2D*> (dataHist->Project3D("yx2e"));
275
276   dataHist->GetZaxis()->SetRange(0, 0);
277   vtxVsEta->GetXaxis()->SetTitle(dataHist->GetXaxis()->GetTitle());
278   vtxVsEta->GetYaxis()->SetTitle(dataHist->GetYaxis()->GetTitle());
279
280   if (vtxVsEta == 0)
281   {
282     printf("ERROR: pt integration failed\n");
283     return;
284   }
285
286   const Float_t vertexRange = 9.99;
287
288   for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
289   {
290     // do we have several histograms for different vertex positions?
291     Int_t vertexBinGlobalBegin = vertexHist->GetXaxis()->FindBin(-vertexRange);
292     Int_t vertexBinWidth = (vertexHist->GetXaxis()->FindBin(vertexRange) - vertexBinGlobalBegin + 1) / (kVertexBinning-1);
293     //printf("vertexBinGlobalBegin = %d, vertexBinWidth = %d\n", vertexBinGlobalBegin, vertexBinWidth);
294     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
295     {
296
297       Int_t vertexBinBegin = vertexBinGlobalBegin;
298       Int_t vertexBinEnd   = vertexBinGlobalBegin + vertexBinWidth * (kVertexBinning-1);
299
300       // the first histogram is always for the whole vertex range
301       if (vertexPos > 0)
302       {
303         vertexBinBegin = vertexBinGlobalBegin + vertexBinWidth * (vertexPos-1);
304         vertexBinEnd = vertexBinBegin + vertexBinWidth;
305       }
306
307       //printf("vertexBinBegin = %d, vertexBinEnd = %d\n", vertexBinBegin, vertexBinEnd);
308
309       Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
310       if (totalEvents == 0)
311       {
312         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
313         continue;
314       }
315
316       Float_t sum = 0;
317       Float_t sumError2 = 0;
318       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
319       {      
320         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
321         {
322           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
323
324           if (sumError2 > 10e30)
325             printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
326
327           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
328         }
329       }
330
331       Float_t ptCutOffCorrection = 1;
332       if (correction && ptCut > 0)
333         ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
334
335       if (ptCutOffCorrection <= 0)
336       {
337         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
338         continue;
339       }
340
341       //printf("Eta: %d Vertex Range: %d %d, Event Count %f, Track Sum: %f, Track Sum corrected: %f\n", iEta, vertexBinBegin, vertexBinEnd, totalEvents, sum, sum / ptCutOffCorrection);
342
343       Float_t dndeta = sum / totalEvents;
344       Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
345
346       dndeta = dndeta/fdNdEta[vertexPos]->GetBinWidth(iEta);
347       error  = error/fdNdEta[vertexPos]->GetBinWidth(iEta);
348
349       fdNdEta[vertexPos]->SetBinContent(iEta, dndeta);
350       fdNdEta[vertexPos]->SetBinError(iEta, error);
351
352       dndeta /= ptCutOffCorrection;
353       error  /= ptCutOffCorrection;
354
355       fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(iEta, dndeta);
356       fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(iEta, error);
357
358     }
359   }
360 }
361
362 //____________________________________________________________________
363 void dNdEtaAnalysis::SaveHistograms()
364 {
365   // save the histograms to a directory with the name of this class (retrieved from TNamed)
366
367   gDirectory->mkdir(GetName());
368   gDirectory->cd(GetName());
369
370   if (fData)
371   {
372     fData->SaveHistograms();
373   }
374   else
375     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
376
377   if (fPtDist)
378     fPtDist       ->Write();
379   else
380     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
381
382   for (Int_t i=0; i<kVertexBinning; ++i)
383   {
384     if (fdNdEta[i])
385       fdNdEta[i]->Write();
386     else
387       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
388
389     if (fdNdEtaPtCutOffCorrected[i])
390       fdNdEtaPtCutOffCorrected[i]->Write();
391     else
392       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
393   }
394
395   gDirectory->cd("../");
396 }
397
398 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
399 {
400   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
401
402   if (!dir)
403     dir = GetName();
404
405   gDirectory->cd(dir);
406
407   fData->LoadHistograms();
408
409   for (Int_t i=0; i<kVertexBinning; ++i)
410   {
411     fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
412     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
413   }
414
415   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
416
417   gDirectory->cd("../");
418 }
419
420 //____________________________________________________________________
421 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
422 {
423   // draws the histograms
424   
425   if (!simple)
426   {
427     if (fData)
428       fData->DrawHistograms(GetName());
429
430     TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
431     canvas->Divide(2, 1);
432
433     canvas->cd(1);
434     if (fdNdEtaPtCutOffCorrected[0])
435       fdNdEtaPtCutOffCorrected[0]->Draw();
436
437     if (fdNdEta[0])
438     {
439       fdNdEta[0]->SetLineColor(kRed);
440       fdNdEta[0]->Draw("SAME");
441     }
442
443     canvas->cd(2);
444     if (fPtDist)
445       fPtDist->Draw();
446   }
447
448     // histograms for different vertices?
449   if (kVertexBinning > 0)
450   {
451       // doesnt work, but i dont get it, giving up...
452     TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
453     TCanvas* canvas3 = 0;
454     if (!simple)
455       canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
456
457     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
458     //printf("%d\n", yPads);
459     //canvas2->Divide(2, yPads);
460
461     TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
462
463     for (Int_t i=0; i<kVertexBinning; ++i)
464     {
465       if (fdNdEtaPtCutOffCorrected[i])
466       {
467         canvas2->cd();
468
469         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
470         fdNdEtaPtCutOffCorrected[i]->Draw((i == 0) ? "" : "SAME");
471         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
472       }
473       if (canvas3 && fdNdEta[i])
474       {
475         canvas3->cd();
476
477         fdNdEta[i]->SetLineColor(i+1);
478         fdNdEta[i]->Draw((i == 0) ? "" : "SAME");
479       }
480     }
481
482     canvas2->cd();
483     legend->Draw();
484     canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
485
486     if (canvas3)
487     {
488       canvas3->cd();
489       legend->Draw();
490     }
491   }
492   
493   if (kVertexBinning == 3)
494   {
495      TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
496      TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
497      
498      if (clone && clone2)
499      {
500         TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
501     
502         clone->Divide(fdNdEtaPtCutOffCorrected[0]);
503         clone->GetYaxis()->SetRangeUser(0.95, 1.05);
504         clone->Draw();
505         
506         clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
507         clone2->Draw("SAME");
508
509         TLine* line = new TLine(-1, 1, 1, 1);
510         line->Draw();
511
512         canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
513      }
514    }
515 }
516
517 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
518 {
519   // Merges a list of dNdEtaAnalysis objects with this one.
520   // This is needed for PROOF.
521   // Returns the number of merged objects (including this)
522
523   if (!list)
524     return 0;
525
526   if (list->IsEmpty())
527     return 1;
528
529   TIterator* iter = list->MakeIterator();
530   TObject* obj;
531
532   // sub collections
533   const Int_t nCollections = 2 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
534   TList* collections[nCollections];
535   for (Int_t i=0; i<nCollections; ++i)
536     collections[i] = new TList;
537
538   Int_t count = 0;
539   while ((obj = iter->Next()))
540   {
541     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
542     if (entry == 0)
543       continue;
544
545     collections[0]->Add(entry->fData);
546     collections[1]->Add(entry->fPtDist);
547
548     for (Int_t i=0; i<kVertexBinning; ++i)
549     {
550       collections[2+i]->Add(entry->fdNdEta[i]);
551       collections[2+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
552     }
553
554     ++count;
555   }
556
557   fData->Merge(collections[0]);
558   fPtDist->Merge(collections[1]);
559   for (Int_t i=0; i<kVertexBinning; ++i)
560   {
561     fdNdEta[i]->Merge(collections[2+i]);
562     fdNdEta[i]->Merge(collections[2+kVertexBinning+i]);
563   }
564
565   for (Int_t i=0; i<nCollections; ++i)
566     delete collections[i];
567
568   return count+1;
569 }