10ede25e75df4fc6ed334696b9c3c7a3970c45d9
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrection.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaCorrection.h"
4
5 #include <TCanvas.h>
6 #include <TH3F.h>
7 #include <TH1D.h>
8
9 //____________________________________________________________________
10 ClassImp(AlidNdEtaCorrection)
11
12 //____________________________________________________________________
13 AlidNdEtaCorrection::AlidNdEtaCorrection(Char_t* name) 
14   : TNamed(name, name),
15   fNEvents(0),
16   fNTriggeredEvents(0)
17 {  
18   // constructor
19   //
20
21   Float_t binLimitsPt[] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 10.0, 100.0};
22
23   fTrack2ParticleCorrection = new AliCorrectionMatrix3D("nTrackToNPart", "nTrackToNPart", 40, -20, 20, 60, -6, 6, 14, binLimitsPt);
24
25   Float_t binLimitsN[]   = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
26                             10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
27   Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
28
29   fVertexRecoCorrection        = new AliCorrectionMatrix2D("vtxReco",       "vtxReco",10,binLimitsVtx ,22,binLimitsN);
30   fTriggerCorrection           = new AliCorrectionMatrix2D("trigger",       "trigger",10,binLimitsVtx ,22,binLimitsN);
31
32   fTriggerBiasCorrection       = new AliCorrectionMatrix2D("triggerBias",   "triggerBias",120,-6,6,100, 0, 10);
33
34   fTrack2ParticleCorrection ->SetAxisTitles("vtx z [cm]", "#eta", "p_{T} [GeV/c]");
35   fVertexRecoCorrection        ->SetAxisTitles("vtx z [cm]", "Ntracks");
36   fTriggerCorrection        ->SetAxisTitles("vtx z [cm]", "Ntracks");
37
38   fTriggerBiasCorrection       ->SetAxisTitles("#eta", "p_{T} [GeV/c]");
39 }
40
41 //____________________________________________________________________
42 void
43 AlidNdEtaCorrection::Finish() {
44   //
45   // finish method
46   //
47   // divide the histograms in the AliCorrectionMatrix2D objects to get the corrections
48
49   fTrack2ParticleCorrection->Divide();
50
51   TH3F* hist = fTrack2ParticleCorrection->GetCorrectionHistogram();
52   Int_t emptyBins = 0;
53   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
54     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
55       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
56         if (hist->GetBinContent(x, y, z) == 0)
57         {
58           printf("Empty bin in fTrack2ParticleCorrection at vtx = %f, eta = %f, pt = %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
59           ++emptyBins;
60         }
61
62   printf("INFO: In the central region fTrack2ParticleCorrection has %d empty bins\n", emptyBins);
63
64   fVertexRecoCorrection->Divide();
65   fTriggerCorrection->Divide();
66
67   if (fNEvents == 0)
68   {
69     printf("ERROR: fNEvents is empty. Cannot scale histogram. Skipping processing of fTriggerBiasCorrection\n");
70     return;
71   }
72   fTriggerBiasCorrection->GetMeasuredHistogram()->Scale(Double_t(fNTriggeredEvents)/Double_t(fNEvents));
73   fTriggerBiasCorrection->Divide();
74 }
75
76 //____________________________________________________________________
77 Long64_t 
78 AlidNdEtaCorrection::Merge(TCollection* list) {
79   // Merge a list of dNdEtaCorrection objects with this (needed for
80   // PROOF). 
81   // Returns the number of merged objects (including this).
82
83   if (!list)
84     return 0;
85   
86   if (list->IsEmpty())
87     return 1;
88
89   TIterator* iter = list->MakeIterator();
90   TObject* obj;
91
92   // collections of measured and generated histograms
93   TList* collectionNtrackToNparticle = new TList;
94   TList* collectionVertexReco        = new TList;
95   TList* collectionTriggerBias       = new TList;
96
97   Int_t count = 0;
98   while ((obj = iter->Next())) {
99     
100     AlidNdEtaCorrection* entry = dynamic_cast<AlidNdEtaCorrection*> (obj);
101     if (entry == 0) 
102       continue;
103
104     collectionNtrackToNparticle ->Add(entry->GetTrack2ParticleCorrection());
105     collectionVertexReco        ->Add(entry->GetVertexRecoCorrection());
106     collectionTriggerBias        ->Add(entry->GetTriggerBiasCorrection());
107
108     count++;
109   }
110   fTrack2ParticleCorrection ->Merge(collectionNtrackToNparticle);
111   fVertexRecoCorrection        ->Merge(collectionVertexReco);
112   fTriggerBiasCorrection        ->Merge(collectionTriggerBias);
113   
114   delete collectionNtrackToNparticle;
115   delete collectionVertexReco;
116   delete collectionTriggerBias;
117
118   return count+1;
119 }
120
121
122 //____________________________________________________________________
123 Bool_t
124 AlidNdEtaCorrection::LoadHistograms(Char_t* fileName, Char_t* dir) {
125   //
126   // loads the histograms
127   //
128
129   fTrack2ParticleCorrection ->LoadHistograms(fileName, dir);
130   fVertexRecoCorrection        ->LoadHistograms(fileName, dir);
131   fTriggerCorrection        ->LoadHistograms(fileName, dir);
132   fTriggerBiasCorrection       ->LoadHistograms(fileName, dir);
133
134   return kTRUE;
135 }
136
137
138 //____________________________________________________________________
139 void
140 AlidNdEtaCorrection::SaveHistograms() {
141   //
142   // save the histograms
143   //
144
145   gDirectory->mkdir(fName.Data());
146   gDirectory->cd(fName.Data());
147
148   fTrack2ParticleCorrection->SaveHistograms();
149   fVertexRecoCorrection->SaveHistograms();
150   fTriggerCorrection->SaveHistograms();
151   fTriggerBiasCorrection->SaveHistograms();
152
153   gDirectory->cd("../");
154 }
155
156 //____________________________________________________________________
157 void AlidNdEtaCorrection::DrawHistograms()
158 {
159   //
160   // call the draw histogram method of the two AliCorrectionMatrix2D objects
161
162   fTrack2ParticleCorrection ->DrawHistograms();
163   fVertexRecoCorrection        ->DrawHistograms();
164   fTriggerCorrection        ->DrawHistograms();
165   fTriggerBiasCorrection       ->DrawHistograms();
166
167 }
168
169 //____________________________________________________________________
170 Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
171 {
172   // calculates the fraction of particles measured (some are missed due to the pt cut off)
173   // uses the generated particle histogram from fTrack2ParticleCorrection
174
175   const TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
176
177   // find eta borders, if eta is negative assume -0.8 ... 0.8
178   Int_t etaBegin = 0;
179   Int_t etaEnd = 0;
180   if (eta < 0)
181   {
182     etaBegin = generated->GetYaxis()->FindBin(-0.8);
183     etaEnd = generated->GetYaxis()->FindBin(0.8);
184   }
185   else
186   {
187     etaBegin = generated->GetYaxis()->FindBin(eta);
188     etaEnd = etaBegin;
189   }
190
191   Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
192   Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
193
194   TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", generated->GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
195   ptProj->GetXaxis()->SetTitle(generated->GetZaxis()->GetTitle());
196
197   Int_t ptBin = ptProj->FindBin(ptCutOff);
198   Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
199   Float_t all = ptProj->Integral();
200
201   if (all == 0)
202     return -1;
203
204   Float_t fraction = abovePtCut / all;
205
206   if (debug)
207   {
208     new TCanvas;
209     ptProj->Draw();
210   }
211   else
212     delete ptProj;
213
214   return fraction;
215 }
216