b693bfbda82e5613421388f155f3f34bd4c6a9d7
[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, AliPWG0Helper::AnalysisMode analysisMode) :
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), analysisMode);
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   //new TCanvas; vtxVsEta->DrawCopy();
287   //vtxVsEta->Rebin2D(1, 4);
288
289   const Float_t vertexRange = 9.99;
290
291   for (Int_t iEta=1; iEta<=vtxVsEta->GetNbinsY(); iEta++)
292   {
293     // do we have several histograms for different vertex positions?
294     Int_t vertexBinGlobalBegin = vertexHist->GetXaxis()->FindBin(-vertexRange);
295     Int_t vertexBinWidth = (vertexHist->GetXaxis()->FindBin(vertexRange) - vertexBinGlobalBegin + 1) / (kVertexBinning-1);
296     //printf("vertexBinGlobalBegin = %d, vertexBinWidth = %d\n", vertexBinGlobalBegin, vertexBinWidth);
297     for (Int_t vertexPos=0; vertexPos<kVertexBinning; ++vertexPos)
298     {
299
300       Int_t vertexBinBegin = vertexBinGlobalBegin;
301       Int_t vertexBinEnd   = vertexBinGlobalBegin + vertexBinWidth * (kVertexBinning-1);
302
303       // the first histogram is always for the whole vertex range
304       if (vertexPos > 0)
305       {
306         vertexBinBegin = vertexBinGlobalBegin + vertexBinWidth * (vertexPos-1);
307         vertexBinEnd = vertexBinBegin + vertexBinWidth;
308       }
309
310       //printf("vertexBinBegin = %d, vertexBinEnd = %d\n", vertexBinBegin, vertexBinEnd);
311
312       Float_t totalEvents = vertexHist->Integral(vertexBinBegin, vertexBinEnd - 1);
313       if (totalEvents == 0)
314       {
315         printf("WARNING: No events for hist %d %d %d\n", vertexPos, vertexBinBegin, vertexBinEnd);
316         continue;
317       }
318
319       Float_t sum = 0;
320       Float_t sumError2 = 0;
321       for (Int_t iVtx = vertexBinBegin; iVtx < vertexBinEnd; iVtx++)
322       {      
323         if (vtxVsEta->GetBinContent(iVtx, iEta) != 0)
324         {
325           sum = sum + vtxVsEta->GetBinContent(iVtx, iEta);
326
327         if (sumError2 > 10e30)
328           printf("WARNING: sum of error2 is dangerously large - be prepared for crash... ");
329
330           sumError2 = sumError2 + TMath::Power(vtxVsEta->GetBinError(iVtx, iEta),2);
331         }
332       }
333
334       Float_t ptCutOffCorrection = 1;
335       if (correction && ptCut > 0)
336         ptCutOffCorrection = correction->GetMeasuredFraction(correctionType, ptCut, vtxVsEta->GetYaxis()->GetBinCenter(iEta));
337
338       if (ptCutOffCorrection <= 0)
339       {
340         printf("UNEXPECTED: ptCutOffCorrection is %f for hist %d %d %d\n", ptCutOffCorrection, vertexPos, vertexBinBegin, vertexBinEnd);
341         continue;
342       }
343
344       //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);
345
346       Int_t bin = fdNdEta[vertexPos]->FindBin(vtxVsEta->GetYaxis()->GetBinCenter(iEta));
347       if (bin > 0 && bin < fdNdEta[vertexPos]->GetNbinsX())
348       {
349         Float_t dndeta = sum / totalEvents;
350         Float_t error  = TMath::Sqrt(sumError2) / totalEvents;
351
352         dndeta = dndeta / fdNdEta[vertexPos]->GetBinWidth(bin);
353         error  = error / fdNdEta[vertexPos]->GetBinWidth(bin);
354
355         fdNdEta[vertexPos]->SetBinContent(bin, dndeta);
356         fdNdEta[vertexPos]->SetBinError(bin, error);
357
358         dndeta /= ptCutOffCorrection;
359         error  /= ptCutOffCorrection;
360
361         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinContent(bin, dndeta);
362         fdNdEtaPtCutOffCorrected[vertexPos]->SetBinError(bin, error);
363
364         //Printf("Bin %d has dN/deta = %f", bin, dndeta);
365       }
366     }
367   }
368
369   //new TCanvas; fdNdEta[0]->DrawCopy();
370 }
371
372 //____________________________________________________________________
373 void dNdEtaAnalysis::SaveHistograms()
374 {
375   // save the histograms to a directory with the name of this class (retrieved from TNamed)
376
377   gDirectory->mkdir(GetName());
378   gDirectory->cd(GetName());
379
380   if (fData)
381   {
382     fData->SaveHistograms();
383   }
384   else
385     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fData is 0\n");
386
387   if (fPtDist)
388     fPtDist       ->Write();
389   else
390     printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fPtDist is 0\n");
391
392   for (Int_t i=0; i<kVertexBinning; ++i)
393   {
394     if (fdNdEta[i])
395       fdNdEta[i]->Write();
396     else
397       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEta[%d] is 0\n", i);
398
399     if (fdNdEtaPtCutOffCorrected[i])
400       fdNdEtaPtCutOffCorrected[i]->Write();
401     else
402       printf("dNdEtaAnalysis::SaveHistograms: UNEXPECTED: fdNdEtaPtCutOffCorrected[%d] is 0\n", i);
403   }
404
405   gDirectory->cd("../");
406 }
407
408 void dNdEtaAnalysis::LoadHistograms(const Char_t* dir)
409 {
410   // loads the histograms from a directory with the name of this class (retrieved from TNamed)
411
412   if (!dir)
413     dir = GetName();
414
415   gDirectory->cd(dir);
416
417   fData->LoadHistograms();
418
419   for (Int_t i=0; i<kVertexBinning; ++i)
420   {
421     fdNdEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEta[i]->GetName()));
422     fdNdEtaPtCutOffCorrected[i] = dynamic_cast<TH1F*> (gDirectory->Get(fdNdEtaPtCutOffCorrected[i]->GetName()));
423   }
424
425   fPtDist = dynamic_cast<TH1F*> (gDirectory->Get(fPtDist->GetName()));
426
427   gDirectory->cd("../");
428 }
429
430 //____________________________________________________________________
431 void dNdEtaAnalysis::DrawHistograms(Bool_t simple)
432 {
433   // draws the histograms
434   
435   if (!simple)
436   {
437     if (fData)
438       fData->DrawHistograms(GetName());
439
440     TCanvas* canvas = new TCanvas(Form("%s_dNdEtaAnalysis", GetName()), Form("%s_dNdEtaAnalysis", GetName()), 800, 400);
441     canvas->Divide(2, 1);
442
443     canvas->cd(1);
444     if (fdNdEtaPtCutOffCorrected[0])
445       fdNdEtaPtCutOffCorrected[0]->DrawCopy();
446
447     if (fdNdEta[0])
448     {
449       fdNdEta[0]->SetLineColor(kRed);
450       fdNdEta[0]->DrawCopy("SAME");
451     }
452
453     canvas->cd(2);
454     if (fPtDist)
455       fPtDist->DrawCopy();
456   }
457
458     // histograms for different vertices?
459   if (kVertexBinning > 0)
460   {
461       // doesnt work, but i dont get it, giving up...
462     TCanvas* canvas2 = new TCanvas(Form("%s_dNdEtaAnalysisVtx", GetName()), Form("%s_dNdEtaAnalysisVtx", GetName()), 450, 450);
463     TCanvas* canvas3 = 0;
464     if (!simple)
465       canvas3 = new TCanvas(Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), Form("%s_dNdEtaAnalysisVtx_noptcutoff", GetName()), 450, 450);
466
467     //Int_t yPads = (Int_t) TMath::Ceil(((Double_t) kVertexBinning - 1) / 2);
468     //printf("%d\n", yPads);
469     //canvas2->Divide(2, yPads);
470
471     TLegend* legend = new TLegend(0.4, 0.2, 0.6, 0.4);
472
473     for (Int_t i=0; i<kVertexBinning; ++i)
474     {
475       if (fdNdEtaPtCutOffCorrected[i])
476       {
477         canvas2->cd();
478
479         fdNdEtaPtCutOffCorrected[i]->SetLineColor(i+1);
480         fdNdEtaPtCutOffCorrected[i]->DrawCopy((i == 0) ? "" : "SAME");
481         legend->AddEntry(fdNdEtaPtCutOffCorrected[i], (i == 0) ? "Vtx All" : Form("Vtx Bin %d", i-1));
482       }
483       if (canvas3 && fdNdEta[i])
484       {
485         canvas3->cd();
486
487         fdNdEta[i]->SetLineColor(i+1);
488         fdNdEta[i]->DrawCopy((i == 0) ? "" : "SAME");
489       }
490     }
491
492     canvas2->cd();
493     legend->Draw();
494     canvas2->SaveAs(Form("%s_%s.gif", canvas2->GetName(), GetName()));
495
496     if (canvas3)
497     {
498       canvas3->cd();
499       legend->Draw();
500     }
501   }
502   
503   if (kVertexBinning == 3)
504   {
505      TH1* clone = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[1]->Clone("clone"));
506      TH1* clone2 = dynamic_cast<TH1*> (fdNdEtaPtCutOffCorrected[2]->Clone("clone2"));
507      
508      if (clone && clone2)
509      {
510         TCanvas* canvas4 = new TCanvas(Form("%s_dNdEtaAnalysisVtxRatios", GetName()), Form("%s_dNdEtaAnalysisVtxRatios", GetName()), 450, 450);
511     
512         clone->Divide(fdNdEtaPtCutOffCorrected[0]);
513         clone->GetYaxis()->SetRangeUser(0.95, 1.05);
514         clone->DrawCopy();
515         
516         clone2->Divide(fdNdEtaPtCutOffCorrected[0]);
517         clone2->DrawCopy("SAME");
518
519         TLine* line = new TLine(-1, 1, 1, 1);
520         line->Draw();
521
522         canvas4->SaveAs(Form("%s_%s.gif", canvas4->GetName(), GetName()));
523      }
524    }
525 }
526
527 Long64_t dNdEtaAnalysis::Merge(TCollection* list)
528 {
529   // Merges a list of dNdEtaAnalysis objects with this one.
530   // This is needed for PROOF.
531   // Returns the number of merged objects (including this)
532
533   if (!list)
534     return 0;
535
536   if (list->IsEmpty())
537     return 1;
538
539   TIterator* iter = list->MakeIterator();
540   TObject* obj;
541
542   // sub collections
543   const Int_t nCollections = 2 * kVertexBinning + 2; // 2 standalone hists, two arrays of size kVertexBinning
544   TList* collections[nCollections];
545   for (Int_t i=0; i<nCollections; ++i)
546     collections[i] = new TList;
547
548   Int_t count = 0;
549   while ((obj = iter->Next()))
550   {
551     dNdEtaAnalysis* entry = dynamic_cast<dNdEtaAnalysis*> (obj);
552     if (entry == 0)
553       continue;
554
555     collections[0]->Add(entry->fData);
556     collections[1]->Add(entry->fPtDist);
557
558     for (Int_t i=0; i<kVertexBinning; ++i)
559     {
560       collections[2+i]->Add(entry->fdNdEta[i]);
561       collections[2+kVertexBinning+i]->Add(entry->fdNdEtaPtCutOffCorrected[i]);
562     }
563
564     ++count;
565   }
566
567   fData->Merge(collections[0]);
568   fPtDist->Merge(collections[1]);
569   for (Int_t i=0; i<kVertexBinning; ++i)
570   {
571     fdNdEta[i]->Merge(collections[2+i]);
572     fdNdEta[i]->Merge(collections[2+kVertexBinning+i]);
573   }
574
575   for (Int_t i=0; i<nCollections; ++i)
576     delete collections[i];
577
578   return count+1;
579 }