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