]>
Commit | Line | Data |
---|---|---|
c4df3c2a | 1 | /************************************************************************* |
2 | * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | //////////////////////////////////////////////////////////////////////////// | |
16 | // Class to perform pt-spectra (and ptMin-spectra) extraction of mothers // | |
17 | // particles starting from measured pt-spectra of daughter particles // | |
18 | // that come from inclusive decays. // | |
19 | // E.g.: B->J/psi+X , B->e+X, B->D0+X, etc. // | |
20 | // // | |
21 | // In order to use this class, one first has to run a simulation // | |
22 | // (only kinematics) of the decay channel under study. The analysis // | |
23 | // can be runned using the class AliAnalysisTaskPtMothFromPtDaugh // | |
24 | // which loops on events to create a TNtupla that stores just // | |
25 | // kinematics informations for mothers and daughters of the decay // | |
26 | // under study (this is made in order to speed up). // | |
27 | // // | |
28 | // Therefore the standard inputs of this class are: // | |
29 | // (1) The TNtupla (created by the task using a TChain of galice.root) // | |
30 | // (2) pT-spectrum of the daughter particles // | |
31 | // // | |
32 | // Output would be the pT (and pTMin) spectrum of the mother, based // | |
33 | // on the correction factors computed from the Kinematics.root files // | |
34 | // // | |
35 | // // | |
36 | // Authors: Giuseppe E. Bruno & Fiorella Fionda // | |
37 | // (Giuseppe.Bruno@ba.infn.it) (Fiorella.Fionda@ba.infn.it) // | |
38 | //////////////////////////////////////////////////////////////////////////// | |
39 | ||
40 | #include "TH1F.h" | |
41 | #include "TNtuple.h" | |
42 | #include "TFile.h" | |
43 | #include "TParticle.h" | |
44 | #include "TArrayI.h" | |
45 | ||
46 | #include "AliPtMothFromPtDaugh.h" | |
47 | #include "AliStack.h" | |
48 | #include "AliLog.h" | |
49 | ||
50 | ClassImp(AliPtMothFromPtDaugh) | |
51 | //________________________________________________________________ | |
52 | AliPtMothFromPtDaugh::AliPtMothFromPtDaugh() : | |
53 | TNamed("AliPtMoth","AliPtMoth"), | |
54 | fDecayKine(0x0), | |
55 | fWij(0x0), | |
56 | fFi(0x0), | |
57 | fWijMin(0x0), | |
58 | fFiMin(0x0), | |
59 | fHistoPtDaughter(0x0), | |
60 | fHistoPtMothers(0x0), | |
61 | fHistoPtMinMothers(0x0), | |
62 | fMothers(0x0), | |
63 | fDaughter(0), | |
64 | fyMothMax(0), | |
65 | fyMothMin(0), | |
66 | fyDaughMax(0), | |
67 | fyDaughMin(0), | |
68 | fUseEta(kFALSE), | |
69 | fAnalysisMode(kUserAnalysis) | |
70 | { | |
71 | // | |
72 | // Default constructor | |
73 | // | |
74 | } | |
75 | ||
76 | //________________________________________________________________ | |
77 | AliPtMothFromPtDaugh::AliPtMothFromPtDaugh(const char* name, const char* title) : | |
78 | TNamed(name,title), | |
79 | fDecayKine(0x0), | |
80 | fWij(0x0), | |
81 | fFi(0x0), | |
82 | fWijMin(0x0), | |
83 | fFiMin(0x0), | |
84 | fHistoPtDaughter(0x0), | |
85 | fHistoPtMothers(0x0), | |
86 | fHistoPtMinMothers(0x0), | |
87 | fMothers(0x0), | |
88 | fDaughter(0), | |
89 | fyMothMax(0), | |
90 | fyMothMin(0), | |
91 | fyDaughMax(0), | |
92 | fyDaughMin(0), | |
93 | fUseEta(kFALSE), | |
94 | fAnalysisMode(kUserAnalysis) | |
95 | { | |
96 | // | |
97 | // Named constructor | |
98 | // | |
99 | } | |
100 | ||
101 | //________________________________________________________________ | |
102 | AliPtMothFromPtDaugh::~AliPtMothFromPtDaugh() | |
103 | { | |
104 | // | |
105 | // Default destructor | |
106 | // | |
107 | if(fDecayKine) {delete fDecayKine;} | |
108 | fDecayKine=0; | |
109 | if(fMothers) {delete fMothers;} | |
110 | fMothers=0; | |
111 | if(fHistoPtMothers) { delete fHistoPtMothers; } | |
112 | fHistoPtMothers=0; | |
113 | if(fHistoPtMinMothers) { delete fHistoPtMinMothers;} | |
114 | fHistoPtMinMothers=0; | |
115 | if(fHistoPtDaughter) { delete fHistoPtDaughter; fHistoPtDaughter=0;} | |
116 | } | |
117 | ||
118 | //______________________________________________________________________________________ | |
119 | AliPtMothFromPtDaugh::AliPtMothFromPtDaugh(const AliPtMothFromPtDaugh& extraction) : | |
120 | TNamed(extraction), | |
121 | fDecayKine(0), | |
122 | fWij(0), | |
123 | fFi(0), | |
124 | fWijMin(0), | |
125 | fFiMin(0), | |
126 | fHistoPtDaughter(0), | |
127 | fHistoPtMothers(0), | |
128 | fHistoPtMinMothers(0), | |
129 | fMothers(0), | |
130 | fDaughter(extraction.fDaughter), | |
131 | fyMothMax(extraction.fyMothMax), | |
132 | fyMothMin(extraction.fyMothMin), | |
133 | fyDaughMax(extraction.fyDaughMax), | |
134 | fyDaughMin(extraction.fyDaughMin), | |
135 | fUseEta(extraction.fUseEta), | |
136 | fAnalysisMode(extraction.fAnalysisMode) | |
137 | { | |
138 | // Copy constructor | |
139 | if(extraction.fHistoPtDaughter) fHistoPtDaughter =(TH1F*)extraction.fHistoPtDaughter->Clone("fHistoPtDaughter_copy"); | |
140 | if(extraction.fHistoPtMothers) fHistoPtMothers = (TH1F*)extraction.fHistoPtMothers->Clone("fHistoPtMothers_copy"); | |
141 | if(extraction.fHistoPtMinMothers) fHistoPtMinMothers = (TH1F*)extraction.fHistoPtMinMothers->Clone("fHistoPtMinMothers_copy"); | |
142 | ||
143 | if(extraction.fWij){ | |
144 | fWij = new Double_t*[2*(fHistoPtMothers->GetNbinsX())]; | |
145 | for(Int_t i=0;i<2*(fHistoPtMothers->GetNbinsX());i++) *(fWij+i)=new Double_t[fHistoPtDaughter->GetNbinsX()]; | |
146 | for(Int_t i=0;i<2*(fHistoPtMothers->GetNbinsX());i++){ | |
147 | for(Int_t j=0;j<fHistoPtDaughter->GetNbinsX();j++){ | |
148 | fWij[i][j]=extraction.fWij[i][j]; | |
149 | } | |
150 | } | |
151 | } | |
152 | ||
153 | if(extraction.fFi){ | |
154 | fFi=new Double_t[2*(fHistoPtMothers->GetNbinsX())]; | |
155 | for(Int_t i=0;i<2*(fHistoPtMothers->GetNbinsX());i++) fFi[i]=extraction.fFi[i]; | |
156 | } | |
157 | ||
158 | if(extraction.fWijMin){ | |
159 | fWijMin = new Double_t*[2*(fHistoPtMinMothers->GetNbinsX())]; | |
160 | for(Int_t i=0;i<2*(fHistoPtMinMothers->GetNbinsX());i++) *(fWijMin+i)=new Double_t[fHistoPtDaughter->GetNbinsX()]; | |
161 | for(Int_t i=0;i<2*(fHistoPtMinMothers->GetNbinsX());i++){ | |
162 | for(Int_t j=0;j<fHistoPtDaughter->GetNbinsX();j++){ | |
163 | fWijMin[i][j]=extraction.fWijMin[i][j]; | |
164 | } | |
165 | } | |
166 | } | |
167 | ||
168 | if(extraction.fFiMin){ | |
169 | fFiMin=new Double_t[2*(fHistoPtMinMothers->GetNbinsX())]; | |
170 | for(Int_t i=0;i<2*(fHistoPtMinMothers->GetNbinsX());i++) fFiMin[i]=extraction.fFiMin[i]; | |
171 | } | |
172 | ||
173 | if(extraction.fDecayKine) fDecayKine = (TNtuple*)extraction.fDecayKine->CloneTree(); | |
174 | ||
175 | if(extraction.fMothers) fMothers = new TArrayI(*(extraction.fMothers)); | |
176 | ||
177 | extraction.Copy(*this); | |
178 | } | |
179 | ||
180 | //______________________________________________________________________________________ | |
181 | AliPtMothFromPtDaugh& AliPtMothFromPtDaugh::operator=(const AliPtMothFromPtDaugh &extraction) | |
182 | { | |
183 | // operator assignment | |
184 | if (this!=&extraction) { | |
185 | this->~AliPtMothFromPtDaugh(); | |
186 | new(this) AliPtMothFromPtDaugh(extraction); | |
187 | } | |
188 | return *this; | |
189 | } | |
190 | ||
191 | //______________________________________________________________________________________ | |
192 | Bool_t AliPtMothFromPtDaugh::CreateWeights() | |
193 | { | |
194 | // Set mothers and daughters pdg codes if not | |
195 | // Put a control if mothers histograms binning-size, rapidity (or pseudoapidity) range | |
196 | // are setting. Read daughter histogram (histName) from the file (pathFileName) | |
197 | // Initialize dimensions for correction factors | |
198 | ||
199 | DeleteWeights(); | |
200 | if(!fMothers) {AliError("Set pdg codes of mothers by SetPdgMothers!\n"); return kFALSE;} | |
201 | if(!fDaughter) {AliError("Set pdg code of daughter by SetPdgDaughter!\n"); return kFALSE;} | |
202 | if(!fHistoPtDaughter) { AliError("Daughter histogram doesn't exist! \n"); return kFALSE;} | |
203 | ||
204 | //Set Rapidity or Pseudorapidity range for mothers if not | |
205 | if(!fyMothMax || !fyMothMin ){ AliError("Set rapidity range or pseudoRapidity range for mothers: use SetYmothers(ymin,ymax) or SetEtaMothers(etamin,etamax)"); return kFALSE;} | |
206 | if(!fyDaughMax || !fyDaughMin){ AliError("Set rapidity range or pseudoRapidity range for daughters:use SetYdaughters(ymin,ymax) or SetEtaDaughters(etamin,etamax)"); return kFALSE;} | |
207 | if(!fHistoPtMothers) {AliError("Call method SetBinsPtMoth to set pT-histogram "); return kFALSE;} | |
208 | if(!fHistoPtMinMothers){AliError("Call method SetBinsPtMinMoth to set pTmin-histogram "); return kFALSE;} | |
209 | ||
210 | Int_t nbinsM=(fHistoPtMothers->GetNbinsX()+2); | |
211 | Int_t nbinsD=(fHistoPtDaughter->GetNbinsX()+2); | |
212 | Int_t nbinsMmin=(fHistoPtMinMothers->GetNbinsX()+2); | |
213 | ||
214 | //Create pointers for weights to reconstruct daughter and mothers pT-spectra | |
215 | fWij=new Double_t*[2*nbinsM]; | |
216 | for(Int_t i=0;i<2*nbinsM;i++) | |
217 | {*(fWij+i)=new Double_t[nbinsD];} | |
218 | fFi=new Double_t[2*nbinsM]; | |
219 | ||
220 | fWijMin=new Double_t*[2*nbinsMmin]; | |
221 | for(Int_t i=0;i<2*nbinsMmin;i++) | |
222 | {*(fWijMin+i)=new Double_t[nbinsD];} | |
223 | fFiMin=new Double_t[2*nbinsMmin]; | |
224 | AliInfo(Form("Pt-mothers distribution: pt_min = %f pt_max=%f n_bins=%d \n", | |
225 | fHistoPtMothers->GetBinLowEdge(1),fHistoPtMothers->GetBinLowEdge(nbinsM-1),nbinsM-2)); | |
226 | AliInfo(Form("PtMinimum-mothers distribution: pt_min = %f pt_max=%f n_bins=%d \n", | |
227 | fHistoPtMinMothers->GetBinLowEdge(1),fHistoPtMinMothers->GetBinLowEdge(nbinsMmin-1),nbinsMmin-2)); | |
228 | AliInfo(Form("Pt-daughters distribution: pt_min = %f pt_max=%f n_bins=%d \n", | |
229 | fHistoPtDaughter->GetBinLowEdge(1),fHistoPtDaughter->GetBinLowEdge(nbinsD-1),nbinsD-2)); | |
230 | return kTRUE; | |
231 | } | |
232 | ||
233 | //______________________________________________________________________________________ | |
234 | void AliPtMothFromPtDaugh::DeleteWeights() | |
235 | { | |
236 | //delete correction factors | |
237 | //delete histogram of daughters | |
238 | if(!fHistoPtMothers || !fHistoPtMinMothers) {AliError("Mothers histograms don't exist! Cannot delete correction factors"); return;} | |
239 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
240 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
241 | if(fWij){ | |
242 | for(Int_t i=0; i<(2*nbinsM); i++) delete fWij[i]; | |
243 | delete [] fWij; fWij=0; | |
244 | } | |
245 | if(fFi) { delete fFi; fFi=0;} | |
246 | if(fWijMin){ | |
247 | for(Int_t i=0; i<(2*nbinsMmin); i++) delete fWijMin[i]; | |
248 | delete [] fWijMin; fWijMin=0; | |
249 | } | |
250 | if(fFiMin) { delete fFiMin; fFiMin=0;} | |
251 | ||
252 | return; | |
253 | } | |
254 | ||
255 | //______________________________________________________________________________________ | |
256 | Bool_t AliPtMothFromPtDaugh::ReadHistoPtDaught(const TH1F *hist) | |
257 | { | |
258 | //Initialize daughter histograms with hist | |
259 | if(!hist) {AliError("Set correct histograms of daughter! It doesn't exist!\n"); return kFALSE;} | |
260 | if(fHistoPtDaughter) delete fHistoPtDaughter; | |
261 | fHistoPtDaughter = (TH1F*)hist->Clone(); | |
262 | return kTRUE; | |
263 | } | |
264 | ||
265 | //______________________________________________________________________________________ | |
266 | Double_t* AliPtMothFromPtDaugh::SetBinsSize(Double_t ptmin, Double_t ptmax,Int_t nbins,Double_t alpha) | |
267 | { | |
268 | // return a pointer of double which contains the binning size for (mothers) histograms: | |
269 | // alpha have to be > 0 and: | |
270 | // alpha = 1 equal binning size | |
271 | // alpha < 1 increasing " | |
272 | // alpha > 1 decreasing " | |
273 | if(ptmin<0 || ptmax<0 || nbins<=0 || alpha<=0) {AliError("Set correct bin-size: ptmin>=0, ptmax>=0, nbins>0, alpha>0! \n"); return 0;} | |
274 | Double_t *edgebin = new Double_t[nbins+1]; | |
275 | Double_t ptmin1=TMath::Power(ptmin,alpha); | |
276 | Double_t ptmax1=TMath::Power(ptmax,alpha); | |
277 | Double_t size=(ptmax1-ptmin1)/nbins; | |
278 | for(Int_t i=0;i<nbins+1;i++) *(edgebin+i)=TMath::Power((ptmin1+i*size),(Double_t)1/alpha); | |
279 | return edgebin; | |
280 | } | |
281 | ||
282 | //______________________________________________________________________________________ | |
283 | void AliPtMothFromPtDaugh::SetBinsPtMoth(Double_t ptmin, Double_t ptmax,Int_t nbins,Double_t alpha) | |
284 | { | |
285 | // Set bin size for pt-spectrum of mothers using SetBinsSize: | |
286 | // alpha have to be > 0 and: | |
287 | // alpha = 1 equal binning size | |
288 | // alpha < 1 increasing " | |
289 | // alpha > 1 decreasing " | |
290 | Double_t* edges = SetBinsSize(ptmin,ptmax,nbins,alpha); | |
291 | SetBinsPtMoth(nbins,edges); | |
b5ed6afb | 292 | delete [] edges; |
c4df3c2a | 293 | return; |
294 | } | |
295 | ||
296 | //______________________________________________________________________________________ | |
297 | void AliPtMothFromPtDaugh::SetBinsPtMoth(Int_t nbins,Double_t *edgeBins) | |
298 | { | |
299 | //set bin size given by the pointer edgeBins for pt-spectrum of mothers: | |
300 | //the dimension of the pointer edgeBins is nbins+1 and the points | |
301 | //has to be written in increasing order | |
302 | if(nbins<0) {AliError("Numbers of bins should be > 0 !\n"); return;} | |
303 | if(fHistoPtMothers) {delete fHistoPtMothers; fHistoPtMothers=0;} | |
304 | fHistoPtMothers=new TH1F("fHistoPtMothers","Reconstructed p_{T}(Mothers)-spectrum",nbins,edgeBins); | |
305 | return; | |
306 | } | |
307 | ||
308 | //______________________________________________________________________________________ | |
309 | void AliPtMothFromPtDaugh::SetBinsPtMinMoth(Int_t nbins,Double_t *edgeBins) | |
310 | { | |
311 | //set bin size given by the pointer edgeBins for ptMin-spectrum of mothers: | |
312 | //the dimension of the pointer edgeBins is nbins+1 and the points | |
313 | //has to be written in increasing order | |
314 | if(nbins<0) {AliError("Numbers of bins should be > 0 !\n"); return;} | |
315 | if(fHistoPtMinMothers) {delete fHistoPtMinMothers; fHistoPtMinMothers=0;} | |
316 | fHistoPtMinMothers = new TH1F("fHistoPtMinMothers","Reconstructed p_{T}^{MIN}(Mothers)-spectrum",nbins,edgeBins); | |
317 | return; | |
318 | } | |
319 | ||
320 | //______________________________________________________________________________________ | |
321 | void AliPtMothFromPtDaugh::SetBinsPtMinMoth(Double_t ptmin, Double_t ptmax,Int_t nbins,Double_t alpha) | |
322 | { | |
323 | // Set bin size for ptMin-spectrum of mothers using SetBinsSize: | |
324 | // alpha have to be > 0 and: | |
325 | // alpha = 1 equal binning size | |
326 | // alpha < 1 increasing " | |
327 | // alpha > 1 decreasing " | |
328 | Double_t* edges = SetBinsSize(ptmin,ptmax,nbins,alpha); | |
329 | SetBinsPtMinMoth(nbins,edges); | |
b5ed6afb | 330 | delete [] edges; |
c4df3c2a | 331 | return; |
332 | } | |
333 | ||
334 | //______________________________________________________________________________________ | |
335 | void AliPtMothFromPtDaugh::SetPdgMothersPrivate(Int_t n_mothers,Int_t *pdgM) | |
336 | { | |
337 | // Set pdg codes of mothers given by the pointer pdgM for the analysis. | |
338 | // This is a private method. | |
339 | if(fMothers) { delete fMothers; fMothers = 0; } | |
340 | fMothers = new TArrayI(n_mothers,pdgM); | |
341 | return; | |
342 | } | |
343 | //______________________________________________________________________________________ | |
344 | void AliPtMothFromPtDaugh::SetPdgMothers(Int_t n_mothers,Int_t *pdgM) | |
345 | { | |
346 | // Set user pdg codes of mothers: first check | |
347 | // that the kUserAnalysis is the selected Analysis_mode. | |
348 | // If not print out a message of error. | |
349 | if(fAnalysisMode!=kUserAnalysis) { | |
350 | AliError("Nothing done: set the mode to kUserAnalysis first"); | |
351 | return; | |
352 | } | |
353 | //Set pdg codes of mothers given by the pointer pdgM for the analysis | |
354 | SetPdgMothersPrivate(n_mothers,pdgM); | |
355 | return; | |
356 | } | |
357 | ||
358 | //______________________________________________________________________________________ | |
359 | void AliPtMothFromPtDaugh::SetBeautyMothers() | |
360 | { | |
361 | // Set pdg codes of beauty particles: | |
362 | // B-mesons (1-24) B-barions (25-59) | |
363 | Int_t pdgBeauty[]={511,521,513,523,10511,10521,10513,10523,20513, | |
364 | 20523,515,525,531,10531,533,10533,205433,535,541,10541,543,10543, | |
365 | 20543,545,5122,5112,5212,5222,5114,5214,5224,5132,5232,5312,5322, | |
366 | 5314,5324,5332,5334,5142,5242,5412,5422,5414,5424,5342,5432,5434, | |
367 | 5442,5444,5512,5522,5514,5524,5532,5534,5542,5544,5554}; | |
368 | Int_t *pdgB=new Int_t[59]; | |
369 | for(Int_t i=0;i<59;i++) pdgB[i]=pdgBeauty[i]; | |
370 | SetPdgMothersPrivate(59,pdgB); | |
371 | } | |
372 | ||
373 | //______________________________________________________________________________________ | |
374 | void AliPtMothFromPtDaugh::InitDefaultAnalysis() | |
375 | { | |
376 | // Set mothers and daughter pdg codes depending from the selected analysis. | |
377 | // case kUserAnalysis: mothers and daughter are set by user (nothing to be done) | |
378 | // case kBtoJPSI: inclusive B-> J/Psi + X channels | |
379 | // case kBtoEle: inclusive B-> e + X channels | |
380 | // case kBtoMuon: inclusive B-> mu + X channels | |
381 | // case kBtoD0: inclusive B-> D0 + X channels | |
382 | ||
383 | switch(fAnalysisMode) | |
384 | { | |
385 | case kUserAnalysis: | |
386 | break; | |
387 | case kBtoJPSI: | |
388 | SetBeautyMothers(); | |
389 | fDaughter = 443; | |
390 | break; | |
391 | case kBtoEle: | |
392 | SetBeautyMothers(); | |
393 | fDaughter = 11; | |
394 | break; | |
395 | case kBtoMuon: | |
396 | SetBeautyMothers(); | |
397 | fDaughter = 13; | |
398 | break; | |
399 | case kBtoD0: | |
400 | SetBeautyMothers(); | |
401 | fDaughter = 421; | |
402 | break; | |
403 | } | |
404 | } | |
405 | ||
406 | //______________________________________________________________________________________ | |
407 | Double_t* AliPtMothFromPtDaugh::GetBinsSize(const TH1F *hist,Int_t &n) const | |
408 | { | |
409 | // return the binning size of the histogram hist | |
410 | // n return the number of bins | |
411 | Double_t* edges = (Double_t*)hist->GetXaxis()->GetXbins()->GetArray(); | |
412 | n=hist->GetNbinsX(); | |
413 | return edges; | |
414 | } | |
415 | ||
416 | //______________________________________________________________________________________ | |
417 | Bool_t AliPtMothFromPtDaugh::GetEtaMothers(Double_t &etaMin, Double_t &etaMax) const | |
418 | { | |
419 | // method to get the bounds of the pseudorapidity range | |
420 | // for mothers. Return kTRUE if pseudorapidity is used and put | |
421 | // pseudorapidity edges in etaMin and etaMax | |
422 | ||
423 | if(fUseEta == kFALSE){ | |
424 | AliError("You are using RAPIDITY range! \n"); | |
425 | etaMin = 0.; | |
426 | etaMax = 0.; | |
427 | return kFALSE; | |
428 | } | |
429 | etaMin = fyMothMin; | |
430 | etaMax = fyMothMax; | |
431 | return kTRUE; | |
432 | } | |
433 | ||
434 | //______________________________________________________________________________________ | |
435 | Bool_t AliPtMothFromPtDaugh::GetEtaDaughter(Double_t &etaMin, Double_t &etaMax) const | |
436 | { | |
437 | // method to get the bounds of the pseudorapidity range | |
438 | // for daughters. Return kTRUE if pseudorapidity is used and put | |
439 | // pseudorapidity edges in etaMin and etaMax | |
440 | ||
441 | if(fUseEta == kFALSE){ | |
442 | AliError("You are using RAPIDITY range! \n"); | |
443 | etaMin = 0.; | |
444 | etaMax = 0.; | |
445 | return kFALSE; | |
446 | } | |
447 | etaMin = fyDaughMin; | |
448 | etaMax = fyDaughMax; | |
449 | return kTRUE; | |
450 | } | |
451 | ||
452 | //______________________________________________________________________________________ | |
453 | Bool_t AliPtMothFromPtDaugh::GetYMothers(Double_t &yMin, Double_t &yMax) const | |
454 | { | |
455 | // method to get the bounds of the rapidity range | |
456 | // for mothers. Return kTRUE if rapidity is used and put | |
457 | // rapidity edges in yMin and yMax | |
458 | ||
459 | if(fUseEta == kTRUE){ | |
460 | AliError("You are using PSEUDORAPIDITY range! \n"); | |
461 | yMin = 0.; | |
462 | yMax = 0.; | |
463 | return kFALSE; | |
464 | } | |
465 | yMin = fyMothMin; | |
466 | yMax = fyMothMax; | |
467 | return kTRUE; | |
468 | } | |
469 | ||
470 | //______________________________________________________________________________________ | |
471 | Bool_t AliPtMothFromPtDaugh::GetYDaughter(Double_t &yMin, Double_t &yMax) const | |
472 | { | |
473 | // method to get the bounds of the rapidity range | |
474 | // for daughters. Return kTRUE if rapidity is used and put | |
475 | // rapidity edges in yMin and yMax | |
476 | ||
477 | if(fUseEta == kTRUE){ | |
478 | AliError("You are using PSEUDORAPIDITY range! \n"); | |
479 | yMin = 0.; | |
480 | yMax = 0.; | |
481 | return kFALSE; | |
482 | } | |
483 | yMin = fyDaughMin; | |
484 | yMax = fyDaughMax; | |
485 | return kTRUE; | |
486 | } | |
487 | ||
488 | ||
489 | //______________________________________________________________________________________ | |
490 | void AliPtMothFromPtDaugh::SetPdgDaugh(Int_t pdgD) | |
491 | { | |
492 | // Set pdg code for daughter particle. Check | |
493 | // that the kUserAnalysis is the selected Analysis_mode. | |
494 | // If not print out a message of error. | |
495 | switch(fAnalysisMode) | |
496 | { | |
497 | case kUserAnalysis: | |
498 | fDaughter = pdgD; | |
499 | break; | |
500 | case kBtoJPSI: | |
501 | if(pdgD!= 443) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} | |
502 | else {fDaughter = pdgD;} | |
503 | break; | |
504 | case kBtoEle: | |
505 | if(pdgD!= 11) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} | |
506 | else {fDaughter = pdgD;} | |
507 | break; | |
508 | case kBtoMuon: | |
509 | if(pdgD!= 13) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} | |
510 | else {fDaughter = pdgD;} | |
511 | break; | |
512 | case kBtoD0: | |
513 | if(pdgD!= 421) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} | |
514 | else {fDaughter = pdgD;} | |
515 | break; | |
516 | } | |
517 | return; | |
518 | } | |
519 | ||
520 | //______________________________________________________________________________________ | |
521 | Int_t* AliPtMothFromPtDaugh::GetPdgMothers(Int_t &n_mothers) const | |
522 | { | |
523 | // return the pointer to the array of pdg codes of mothers particles | |
524 | // if it exist. Put its dimension in n_mothers | |
525 | ||
526 | if(!fMothers) {AliError("Mothers pdg are not defined! \n"); return 0x0;} | |
527 | n_mothers = fMothers->GetSize(); | |
528 | return fMothers->GetArray(); | |
529 | } | |
530 | ||
531 | //______________________________________________________________________________________ | |
532 | Double_t AliPtMothFromPtDaugh::GetW(Int_t i,Int_t j) const | |
533 | { | |
534 | // Return value of correction factors Wij at the position i (pt-mothers bin index)- | |
535 | // j (pt-daughter bin index). Bin 0 is the underflow, bin nbins+1 the overflow. | |
536 | // If Wij don't exist or the indices i or j are out of the variability range return 0 | |
537 | ||
538 | if(!fWij) {AliError("Correction factors Wij are not been evaluated yet!\n"); return 0;} | |
539 | if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} | |
540 | if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} | |
541 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
542 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
543 | if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow ",0,nbinsM-1,nbinsM-1)); return 0;} | |
544 | ||
545 | if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} | |
546 | return fWij[i][j]; | |
547 | } | |
548 | ||
549 | //______________________________________________________________________________________ | |
550 | Double_t AliPtMothFromPtDaugh::GetStatErrW(Int_t i,Int_t j) const | |
551 | { | |
552 | // Return value of statistical error on correction factors Wij at the position | |
553 | // i (pt-mothers bin index)- j (pt-daughter bin index). Bin 0 is the underflow, | |
554 | // bin nbins+1 the overflow. If Wij don't exist or the indices i or j are out of the | |
555 | // variability range return 0 | |
556 | ||
557 | if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} | |
558 | if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} | |
559 | ||
560 | if(!fWij) {AliError("Correction factors Wij are not been evaluated yet!\n"); return 0;} | |
561 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
562 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
563 | if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsM-1,nbinsM-1)); return 0;} | |
564 | if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} | |
565 | return fWij[i+nbinsM][j]; | |
566 | } | |
567 | ||
568 | //______________________________________________________________________________________ | |
569 | Double_t AliPtMothFromPtDaugh::GetF(Int_t i) const | |
570 | { | |
571 | // Return value of correction factors Fi at the position i (pt-mothers bin index). | |
572 | // Bin 0 is the underflow, bin nbins+1 the overflow. | |
573 | // If Fi don't exist or the index i is out of the variability range return 0 | |
574 | ||
575 | if(!fFi) {AliError("Correction factors Fi are not been evaluated yet!\n"); return 0;} | |
576 | if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} | |
577 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
578 | if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsM-1,nbinsM-1)); return 0;} | |
579 | return fFi[i]; | |
580 | } | |
581 | ||
582 | //______________________________________________________________________________________ | |
583 | Double_t AliPtMothFromPtDaugh::GetStatErrF(Int_t i) const | |
584 | { | |
585 | // Return statistical error on correction factors Fi at the position i (pt-mothers bin index). | |
586 | // Bin 0 is the underflow, bin nbins+1 the overflow. | |
587 | // If Fi don't exist or the index i is out of the variability range return 0 | |
588 | ||
589 | if(!fFi) {AliError("Correction factors Fi are not been evaluated yet!\n"); return 0;} | |
590 | if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} | |
591 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
592 | if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsM-1,nbinsM-1)); return 0;} | |
593 | ||
594 | return fFi[i+nbinsM]; | |
595 | } | |
596 | ||
597 | //______________________________________________________________________________________ | |
598 | Double_t AliPtMothFromPtDaugh::GetWmin(Int_t i,Int_t j) const | |
599 | { | |
600 | // Return value of correction factors Wij_min at the position i (ptMin-mothers bin index)- | |
601 | // j (pt-daughter bin index). Bin 0 is the underflow, bin nbins+1 the overflow. | |
602 | // If Wij_min don't exist or the indices i or j are out of the variability range return 0 | |
603 | ||
604 | if(!fWijMin) {AliError("Correction factors Wij_min are not been evaluated yet!\n"); return 0;} | |
605 | if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} | |
606 | if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} | |
607 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
608 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
609 | if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} | |
610 | if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} | |
611 | return fWijMin[i][j]; | |
612 | } | |
613 | ||
614 | //______________________________________________________________________________________ | |
615 | Double_t AliPtMothFromPtDaugh::GetStatErrWmin(Int_t i,Int_t j) const | |
616 | { | |
617 | // Return value of statistical error on correction factors Wij_min at the position | |
618 | // i (ptMin-mothers bin index)- j (pt-daughter bin index). Bin 0 is the underflow, | |
619 | // bin nbins+1 the overflow. If Wij_min don't exist or the indices i or j are out of the | |
620 | // variability range return 0 | |
621 | ||
622 | if(!fWijMin) {AliError("Correction factors Wij_min are not been evaluated yet!\n"); return 0;} | |
623 | if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} | |
624 | if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} | |
625 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
626 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
627 | if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} | |
628 | if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} | |
629 | ||
630 | return fWijMin[i+nbinsMmin][j]; | |
631 | } | |
632 | ||
633 | //______________________________________________________________________________________ | |
634 | Double_t AliPtMothFromPtDaugh::GetFmin(Int_t i) const | |
635 | { | |
636 | // Return value of correction factors Fi_min at the position i (ptMin-mothers bin index). | |
637 | // Bin 0 is the underflow, bin nbins+1 the overflow. | |
638 | // If Fi_min don't exist or the index i is out of the variability range return 0 | |
639 | ||
640 | if(!fFiMin) {AliError("Correction factors Fi_min are not been evaluated yet!\n"); return 0;} | |
641 | if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} | |
642 | Int_t nbinsMmin=fHistoPtMothers->GetNbinsX()+2; | |
643 | if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} | |
644 | return fFiMin[i]; | |
645 | } | |
646 | ||
647 | //______________________________________________________________________________________ | |
648 | Double_t AliPtMothFromPtDaugh::GetStatErrFmin(Int_t i) const | |
649 | { | |
650 | // Return statistical error on correction factors Fi_min at the position i (ptMin-mothers bin index). | |
651 | // Bin 0 is the underflow, bin nbins+1 the overflow. | |
652 | // If Fi_min don't exist or the index i is out of the variability range return 0 | |
653 | ||
654 | if(!fFiMin) {AliError("Correction factors Fi_min are not been evaluated yet!\n"); return 0;} | |
655 | if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} | |
656 | Int_t nbinsMmin=fHistoPtMothers->GetNbinsX()+2; | |
657 | if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} | |
658 | ||
659 | return fFiMin[i+nbinsMmin]; | |
660 | } | |
661 | ||
662 | //______________________________________________________________________________________ | |
663 | Bool_t AliPtMothFromPtDaugh::IsMothers(Int_t pdgCode) | |
664 | { | |
665 | //return kTRUE if pdgCode is in the list of pdg codes of mothers | |
666 | Int_t dim = fMothers->GetSize(); | |
667 | for(Int_t i=0;i<dim;i++) { | |
668 | Int_t pdgMother = (Int_t)fMothers->GetAt(i); | |
669 | if(pdgCode == pdgMother) return kTRUE; } | |
670 | return kFALSE; | |
671 | } | |
672 | ||
673 | //______________________________________________________________________________________ | |
674 | Bool_t AliPtMothFromPtDaugh::IsSelectedDaugh(const TParticle *part, Int_t &labelDaugh, AliStack *const stack) | |
675 | { | |
676 | // return kTRUE if particle part has the selected daughter | |
677 | // if yes put the label of the track in labelDaugh | |
678 | TParticle *daugh; | |
679 | Int_t nDg=part->GetNDaughters(); | |
680 | if(nDg<=0) {AliError("I have no daugh!\n"); return kFALSE;} | |
681 | for(Int_t i=part->GetFirstDaughter();i<part->GetLastDaughter()+1;i++) | |
682 | { | |
683 | daugh=stack->Particle(i); | |
684 | if(TMath::Abs(daugh->GetPdgCode())==fDaughter) { | |
685 | labelDaugh=i; | |
686 | return kTRUE; | |
687 | } | |
688 | } | |
689 | return kFALSE; | |
690 | } | |
691 | ||
692 | //______________________________________________________________________________________ | |
693 | Bool_t AliPtMothFromPtDaugh::Rapidity(const TParticle *particle, Double_t &y) | |
694 | { | |
695 | // Evaluated rapidity of particle and put it in y. Return kFALSE if | |
696 | // cannot compute rapidity | |
697 | y=-999; | |
698 | if(particle->Energy()-particle->Pz()<=0) return kFALSE; | |
699 | y=0.5*TMath::Log((particle->Energy()+particle->Pz())/(particle->Energy()-particle->Pz())); | |
700 | return kTRUE; | |
701 | } | |
702 | ||
703 | //______________________________________________________________________________________ | |
704 | Int_t AliPtMothFromPtDaugh::GiveBinIndex(Double_t Ptpart,const TH1F *ptHist) const | |
705 | { | |
706 | // Return the bin index of pt respect to the binning size of ptHist | |
707 | // bin 0 is the underflow - nbins+1 is the overflow | |
708 | Int_t nbins=ptHist->GetNbinsX(); | |
709 | Int_t index=0; | |
710 | for(Int_t i=1;i<nbins+2;i++) | |
711 | { | |
712 | if(Ptpart<(ptHist->GetBinLowEdge(i))) | |
713 | { | |
714 | index=i-1; | |
715 | break; | |
716 | } | |
717 | } | |
718 | return index; | |
719 | } | |
720 | ||
721 | //______________________________________________________________________________________ | |
722 | Bool_t AliPtMothFromPtDaugh::CutDaugh(Double_t yD,Double_t ptD) | |
723 | { | |
724 | // put a control for rapidity yD and transverse momentum ptD of daughter | |
725 | // return kTRUE if fyDaughMin < yD < fyDaughMax and ptMinDaugh < ptD < ptMaxDaugh | |
726 | Double_t ptMinDaugh = fHistoPtDaughter->GetBinLowEdge(1); | |
727 | Double_t ptMaxDaugh = fHistoPtDaughter->GetBinLowEdge(fHistoPtDaughter->GetNbinsX()); | |
728 | if( yD < fyDaughMin || yD > fyDaughMax ) return kFALSE; | |
729 | if( ptD > ptMaxDaugh || ptD < ptMinDaugh ) return kFALSE; | |
730 | return kTRUE; | |
731 | } | |
732 | ||
733 | //______________________________________________________________________________________ | |
734 | Bool_t AliPtMothFromPtDaugh::EvaluateWij() | |
735 | { | |
736 | // Evaluate correction factors using to extract the ptRaw and | |
737 | // ptMinRaw distributions. Statistical errors on those are computed too | |
738 | ||
739 | if(!fHistoPtMothers || !fHistoPtDaughter || !fHistoPtMinMothers) | |
740 | {AliError("Control mother and daughter histograms!\n"); return kFALSE;} | |
741 | ||
742 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
743 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
744 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
745 | Float_t pdgM, pdgD, pxM,pyM,pzM,yM,pxD,pyD,pzD,yD,etaM,etaD,cutVarM,cutVarD; | |
746 | Double_t* entries=new Double_t[nbinsD]; | |
747 | for(Int_t ii=0;ii<nbinsD;ii++) {entries[ii]=0.;} | |
748 | Int_t i,j,iMin; | |
749 | if(!fDecayKine) {AliError("TNtupla is not defined!\n"); return kFALSE;} | |
750 | fDecayKine->SetBranchAddress("pdgM",&pdgM); | |
751 | fDecayKine->SetBranchAddress("pxM",&pxM); | |
752 | fDecayKine->SetBranchAddress("pyM",&pyM); | |
753 | fDecayKine->SetBranchAddress("pzM",&pzM); | |
754 | fDecayKine->SetBranchAddress("yM",&yM); | |
755 | fDecayKine->SetBranchAddress("etaM",&etaM); | |
756 | fDecayKine->SetBranchAddress("pdgD",&pdgD); | |
757 | fDecayKine->SetBranchAddress("pxD",&pxD); | |
758 | fDecayKine->SetBranchAddress("pyD",&pyD); | |
759 | fDecayKine->SetBranchAddress("pzD",&pzD); | |
760 | fDecayKine->SetBranchAddress("yD",&yD); | |
761 | fDecayKine->SetBranchAddress("etaD",&etaD); | |
762 | Double_t ptD,ptM; | |
763 | // Initialize correction factors for pT and pTmin if those exist | |
764 | if(!fWij) {AliError("Correction factors Wij have not been created!\n"); return kFALSE;} | |
765 | if(!fWijMin) {AliError("Correction factors Wij_min have not been created!\n"); return kFALSE;} | |
766 | for(Int_t ii=0;ii<(2*nbinsM);ii++){ | |
767 | for(Int_t jj=0;jj<nbinsD;jj++){ | |
768 | fWij[ii][jj]=0; | |
769 | } | |
770 | } | |
771 | for(Int_t ii=0;ii<(2*nbinsMmin);ii++){ | |
772 | for(Int_t jj=0;jj<nbinsD;jj++){ | |
773 | fWijMin[ii][jj]=0; | |
774 | } | |
775 | } | |
776 | Int_t nentries = (Int_t)fDecayKine->GetEntries(); | |
777 | Int_t fNcurrent=0; | |
778 | Int_t nb = (Int_t)fDecayKine->GetEvent(fNcurrent); | |
779 | for (Int_t iev=0; iev<nentries; iev++){ | |
780 | // check if rapidity or pseudorapidity range is set | |
781 | if(fUseEta == kFALSE) {cutVarD = yD; cutVarM = yM;} | |
782 | else {cutVarD = etaD; cutVarM = etaM;} | |
783 | ptD=TMath::Sqrt(pxD*pxD+pyD*pyD); | |
784 | ptM=TMath::Sqrt(pxM*pxM+pyM*pyM); | |
785 | pdgM = TMath::Abs(pdgM); | |
786 | if((TMath::Abs(pdgD))==fDaughter && IsMothers((Int_t)pdgM)) | |
787 | { | |
788 | j=GiveBinIndex(ptD,fHistoPtDaughter); | |
789 | i=GiveBinIndex(ptM,fHistoPtMothers); | |
790 | iMin=GiveBinIndex(ptM,fHistoPtMinMothers); | |
791 | if(!CutDaugh(cutVarD,ptD)) { fNcurrent++; nb = (Int_t)fDecayKine->GetEvent(fNcurrent); continue;} | |
792 | if(cutVarM>fyMothMin && cutVarM<fyMothMax){ | |
793 | fWij[i][j]+=1.; | |
794 | for(Int_t k=0;k<iMin+1;k++) {fWijMin[k][j]+=1.;} | |
795 | } | |
796 | entries[j]++; | |
797 | } | |
798 | fNcurrent++; | |
799 | nb = (Int_t)fDecayKine->GetEvent(fNcurrent); | |
800 | } | |
801 | for(Int_t jj=0;jj<nbinsD;jj++){ | |
802 | for(Int_t ii=0;ii<nbinsM;ii++){ | |
803 | if(entries[jj]>0){ | |
804 | fWij[ii][jj]=fWij[ii][jj]/entries[jj]; | |
805 | // evaluate statistical errors on fWij | |
806 | fWij[ii+nbinsM][jj]=(TMath::Sqrt(fWij[ii][jj]*(1-(fWij[ii][jj]/entries[jj]))))/entries[jj]; | |
807 | } | |
808 | else{ | |
809 | // if there are no entries in the bin-j of daughter distribution | |
810 | // set factor = -1 and error = 999 | |
811 | fWij[ii][jj]=-1; | |
812 | fWij[ii+nbinsM][jj]=999; | |
813 | } | |
814 | } | |
815 | for(Int_t ii=0;ii<nbinsMmin;ii++){ | |
816 | if(entries[jj]>0){ | |
817 | fWijMin[ii][jj]=fWijMin[ii][jj]/entries[jj]; | |
818 | //evaluate statistical errors on fWijMin | |
819 | fWijMin[ii+nbinsMmin][jj] = (TMath::Sqrt(fWijMin[ii][jj]*(1-(fWijMin[ii][jj]/entries[jj]))))/entries[jj]; | |
820 | } | |
821 | else{ | |
822 | //if there are no entries set correction factor = -1 and error = -999 | |
823 | fWijMin[ii][jj]=-1; | |
824 | fWijMin[ii+nbinsMmin][jj]=999; | |
825 | } | |
826 | } | |
827 | } | |
b5ed6afb | 828 | delete [] entries; |
c4df3c2a | 829 | return kTRUE; |
830 | } | |
831 | ||
832 | //______________________________________________________________________________________ | |
833 | Bool_t AliPtMothFromPtDaugh::EvaluateFi() | |
834 | { | |
835 | // Evaluate acceptance correction factors that are applied on the | |
836 | // raw distributions. Statistical errors on those are computed too | |
837 | ||
838 | if(!fHistoPtMothers || !fHistoPtMinMothers) | |
839 | {AliError("Control mother histograms!\n"); return kFALSE;} | |
840 | ||
841 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
842 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
843 | Double_t* entries=new Double_t[nbinsM]; | |
844 | Double_t* entries1=new Double_t[nbinsMmin]; | |
845 | if(!fFi) {AliError("Correction factors Fi have not been created!\n"); return kFALSE;} | |
846 | if(!fFiMin) {AliError("Correction factors Fi_min have not been created!\n"); return kFALSE;} | |
847 | //initialize the correction factor for pT and pTmin | |
848 | for(Int_t ii=0;ii<nbinsM;ii++){ | |
849 | fFi[ii]=0.; //for correction factors | |
850 | fFi[ii+nbinsM]=0.; //for statistical error on correction factors | |
851 | entries[ii]=0.; | |
852 | } | |
853 | for(Int_t ii=0;ii<nbinsMmin;ii++){ | |
854 | entries1[ii]=0.; | |
855 | fFiMin[ii]=0.; //for correction factors | |
856 | fFiMin[ii+nbinsMmin]=0.; //for statistical error on correction factors | |
857 | } | |
858 | Float_t pdgM, pdgD, pxM,pyM,pzM,yM,pxD,pyD,pzD,yD,etaM,etaD,cutVarD,cutVarM; | |
859 | Int_t i,iMin; | |
860 | if(!fDecayKine) {AliError("TNtupla is not defined!\n"); return kFALSE;} | |
861 | fDecayKine->SetBranchAddress("pdgM",&pdgM); | |
862 | fDecayKine->SetBranchAddress("pxM",&pxM); | |
863 | fDecayKine->SetBranchAddress("pyM",&pyM); | |
864 | fDecayKine->SetBranchAddress("pzM",&pzM); | |
865 | fDecayKine->SetBranchAddress("yM",&yM); | |
866 | fDecayKine->SetBranchAddress("etaM",&etaM); | |
867 | fDecayKine->SetBranchAddress("pdgD",&pdgD); | |
868 | fDecayKine->SetBranchAddress("pxD",&pxD); | |
869 | fDecayKine->SetBranchAddress("pyD",&pyD); | |
870 | fDecayKine->SetBranchAddress("pzD",&pzD); | |
871 | fDecayKine->SetBranchAddress("yD",&yD); | |
872 | fDecayKine->SetBranchAddress("etaD",&etaD); | |
873 | Double_t ptD,ptM; | |
874 | ||
875 | Int_t nentries = (Int_t)fDecayKine->GetEntries(); | |
876 | Int_t fNcurrent=0; | |
877 | Int_t nb = (Int_t)fDecayKine->GetEvent(fNcurrent); | |
878 | ||
879 | for (Int_t iev=0; iev<nentries; iev++){ | |
880 | pdgM = TMath::Abs(pdgM); | |
881 | if((TMath::Abs(pdgD))==fDaughter && IsMothers((Int_t)pdgM)) | |
882 | { | |
883 | //check if rapidity or pseudorapidity range is set | |
884 | if(fUseEta == kFALSE) {cutVarD = yD; cutVarM = yM;} | |
885 | else {cutVarD = etaD; cutVarM = etaM;} | |
886 | ptD=TMath::Sqrt(pxD*pxD+pyD*pyD); | |
887 | ptM=TMath::Sqrt(pxM*pxM+pyM*pyM); | |
888 | i=GiveBinIndex(ptM,fHistoPtMothers); | |
889 | iMin=GiveBinIndex(ptM,fHistoPtMinMothers); | |
890 | if(cutVarM<fyMothMin || cutVarM>fyMothMax){ | |
891 | fNcurrent++; nb = (Int_t)fDecayKine->GetEvent(fNcurrent); continue;} | |
892 | entries[i]++; | |
893 | for(Int_t k=0; k<iMin+1;k++) {entries1[k]+=1;} | |
894 | if(!CutDaugh(cutVarD,ptD)) {fNcurrent++; nb = (Int_t)fDecayKine->GetEvent(fNcurrent); continue;} | |
895 | fFi[i]+=1.; | |
896 | for(Int_t k=0; k<iMin+1;k++) {fFiMin[k]+=1.;} | |
897 | } | |
898 | fNcurrent++; | |
899 | nb = (Int_t) fDecayKine->GetEvent(fNcurrent); | |
900 | } | |
901 | ||
902 | for(Int_t ii=0;ii<nbinsM;ii++){ | |
903 | if(entries[ii]>0){ | |
904 | fFi[ii]/=entries[ii]; | |
905 | fFi[ii+nbinsM]=(TMath::Sqrt(fFi[ii]*(1-(fFi[ii]/entries[ii]))))/entries[ii]; | |
906 | } | |
907 | else{ | |
908 | fFi[ii]=-1; | |
909 | fFi[ii+nbinsM]=999; | |
910 | } | |
911 | } | |
912 | for(Int_t ii=0;ii<nbinsMmin;ii++){ | |
913 | if(entries1[ii]>0){ | |
914 | fFiMin[ii]/=entries1[ii]; | |
915 | fFiMin[ii+nbinsMmin]=(TMath::Sqrt(fFiMin[ii]*(1-(fFiMin[ii]/entries1[ii]))))/entries1[ii]; | |
916 | } | |
917 | else {fFiMin[ii]=-1; fFiMin[ii+nbinsMmin]=999;} | |
918 | } | |
b5ed6afb | 919 | delete [] entries; |
920 | delete [] entries1; | |
c4df3c2a | 921 | return kTRUE; |
922 | } | |
923 | ||
924 | //______________________________________________________________________________________ | |
925 | Bool_t AliPtMothFromPtDaugh::EvaluatePtMothRaw(TH1F *histoPt, TH1F *histoPtMin) | |
926 | { | |
927 | //Apply the fWij and fWijMin on the daughter distribution | |
928 | //in order to evaluate the pt and ptMin raw distributions for mothers | |
929 | ||
930 | if(!fHistoPtMothers || !fHistoPtDaughter || !fHistoPtMinMothers) | |
931 | {AliError("Control mother and daughter histograms!\n"); return kFALSE;} | |
932 | ||
933 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
934 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
935 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
936 | Double_t lowedge=0.; | |
937 | Double_t wfill=0.; | |
938 | Double_t wMinfill=0.; | |
939 | for(Int_t j=0;j<nbinsD;j++){ | |
940 | for(Int_t i=0;i<nbinsM;i++){ | |
941 | lowedge=fHistoPtMothers->GetBinCenter(i); | |
942 | if(fWij[i][j]>=0) wfill=fWij[i][j]*(fHistoPtDaughter->GetBinContent(j)); | |
943 | histoPt->Fill(lowedge,wfill); | |
944 | } | |
945 | for(Int_t i=0;i<nbinsMmin;i++){ | |
946 | lowedge=fHistoPtMinMothers->GetBinLowEdge(i); | |
947 | if(fWijMin[i][j]>=0) wMinfill=fWijMin[i][j]*(fHistoPtDaughter->GetBinContent(j)); | |
948 | histoPtMin->Fill(lowedge,wMinfill); | |
949 | } | |
950 | } | |
951 | return kTRUE; | |
952 | } | |
953 | ||
954 | //______________________________________________________________________________________ | |
955 | Bool_t AliPtMothFromPtDaugh::EvaluateErrPt(Double_t *erStat) | |
956 | { | |
957 | // Evaluate the statistical error on the pt-mothers distribution. | |
958 | // sigmaX: contribution that came from the measured distibution | |
959 | // sigmaWij: contribution that came from the fWij factors | |
960 | // sigmaFi: contribution that came from the fFi factors | |
961 | ||
962 | if(!fHistoPtMothers || !fHistoPtDaughter) | |
963 | {AliError("Control mother(pt) and daughter histograms!\n"); return kFALSE;} | |
964 | ||
965 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
966 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
967 | Double_t m=0; | |
968 | Double_t sigmaX, sigmaWij, sigmaFi; | |
969 | for(Int_t i=0;i<nbinsM;i++){ | |
970 | sigmaX=0.; | |
971 | sigmaWij=0.; | |
972 | sigmaFi=0; | |
973 | for(Int_t j=0;j<nbinsD;j++){ | |
974 | if(fWij[i][j]>=0){ | |
975 | sigmaX+=(((fWij[i][j])*(fWij[i][j]))*((fHistoPtDaughter->GetBinError(j))*(fHistoPtDaughter->GetBinError(j)))); | |
976 | sigmaWij+=((fHistoPtDaughter->GetBinContent(j))*(fHistoPtDaughter->GetBinContent(j)))*(fWij[i+nbinsM][j]*fWij[i+nbinsM][j]); | |
977 | sigmaFi+=(fWij[i][j])*(fHistoPtDaughter->GetBinContent(j)); | |
978 | } | |
979 | } | |
980 | if(fFi[i]>0) sigmaFi=((sigmaFi*sigmaFi)*(fFi[i+nbinsM]*fFi[i+nbinsM]))/(fFi[i]*fFi[i]); | |
981 | m=TMath::Sqrt(sigmaX+sigmaWij+sigmaFi); | |
982 | if(fFi[i]>0) erStat[i]=(1/fFi[i])*m; | |
983 | else erStat[i]=999; | |
984 | } | |
985 | return kTRUE; | |
986 | } | |
987 | ||
988 | //______________________________________________________________________________________ | |
989 | Bool_t AliPtMothFromPtDaugh::EvaluateErrPtMin(Double_t *erStat) | |
990 | { | |
991 | // Evaluate statistical error on ptMin mothers distribution | |
992 | // sigmaMinX: contribution that came from the measured distibution | |
993 | // sigmaMinWij: contribution that came from the fWijMin factors | |
994 | // sigmaMinFi: contribution that came from the fFiMin factors | |
995 | ||
996 | if(!fHistoPtDaughter || !fHistoPtMinMothers) | |
997 | {AliError("Control mother(ptMin) and daughter histograms!\n"); return kFALSE;} | |
998 | ||
999 | Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; | |
1000 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
1001 | Double_t m1=0; | |
1002 | Double_t sigmaMinX; | |
1003 | Double_t sigmaMinWij; | |
1004 | Double_t sigmaMinFi; | |
1005 | for(Int_t i=0;i<nbinsMmin;i++){ | |
1006 | sigmaMinX=0.; | |
1007 | sigmaMinWij=0.; | |
1008 | sigmaMinFi=0.; | |
1009 | for(Int_t j=0;j<nbinsD;j++){ | |
1010 | if(fWijMin[i][j]>=0){ | |
1011 | sigmaMinX+=(((fWijMin[i][j])*(fWijMin[i][j]))*((fHistoPtDaughter->GetBinError(j))*(fHistoPtDaughter->GetBinError(j)))); | |
1012 | sigmaMinWij+=((fHistoPtDaughter->GetBinContent(j))*(fHistoPtDaughter->GetBinContent(j)))*(fWijMin[i+nbinsMmin][j]*fWijMin[i+nbinsMmin][j]); | |
1013 | sigmaMinFi+=(fWijMin[i][j])*(fHistoPtDaughter->GetBinContent(j)); | |
1014 | } | |
1015 | } | |
1016 | if(fFiMin[i]>0) sigmaMinFi=((sigmaMinFi*sigmaMinFi)*(fFiMin[i+nbinsMmin]*fFiMin[i+nbinsMmin]))/(fFiMin[i]*fFiMin[i]); | |
1017 | m1=TMath::Sqrt(sigmaMinX+sigmaMinWij+sigmaMinFi); | |
1018 | if(fFiMin[i]>0) erStat[i]=(1/fFiMin[i])*m1; | |
1019 | else erStat[i]=999; | |
1020 | } | |
1021 | ||
1022 | return kTRUE; | |
1023 | } | |
1024 | ||
1025 | //______________________________________________________________________________________ | |
1026 | Bool_t AliPtMothFromPtDaugh::EvaluatePtMoth() | |
1027 | { | |
1028 | // Evaluate pt and ptMin distribution for mothers | |
1029 | // First evaluate the sigma raw distribution by calling EvaluatePtMothRaw | |
1030 | // then evaluate the pt and ptMin mothers distribution. | |
1031 | // Statistical errors on those distributions are evaluated too. | |
1032 | ||
1033 | if(!EvaluateWij()) return kFALSE; | |
1034 | if(!EvaluateFi()) return kFALSE; | |
1035 | ||
1036 | // reset pt and ptMin mothers histograms | |
1037 | fHistoPtMothers->Reset(); | |
1038 | fHistoPtMinMothers->Reset(); | |
1039 | ||
1040 | TH1F *histoPt=(TH1F*)fHistoPtMothers->Clone(); | |
1041 | TH1F *histoPtMin=(TH1F*)fHistoPtMinMothers->Clone(); | |
1042 | EvaluatePtMothRaw(histoPt,histoPtMin); | |
1043 | Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; | |
1044 | Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; | |
1045 | Double_t *erPtStat=new Double_t[nbinsM+2]; | |
1046 | EvaluateErrPt(erPtStat); | |
1047 | Double_t *erPtMinStat=new Double_t[nbinsMmin+2]; | |
1048 | EvaluateErrPtMin(erPtMinStat); | |
1049 | Double_t lowedge=0; | |
1050 | Double_t fwfill; | |
1051 | Double_t fMinfill; | |
1052 | for(Int_t i=0;i<nbinsM;i++){ | |
1053 | fwfill=0.; | |
1054 | lowedge=fHistoPtMothers->GetBinCenter(i); | |
1055 | if(fFi[i]>0){ | |
1056 | fwfill=(histoPt->GetBinContent(i))/fFi[i]; | |
1057 | fHistoPtMothers->Fill(lowedge,fwfill); | |
1058 | fHistoPtMothers->SetBinError(i,erPtStat[i]); | |
1059 | } | |
1060 | } | |
1061 | for(Int_t i=0;i<nbinsMmin;i++){ | |
1062 | fMinfill=0.; | |
1063 | lowedge=fHistoPtMinMothers->GetBinCenter(i); | |
1064 | if(fFiMin[i]>0){ | |
1065 | fMinfill=(histoPtMin->GetBinContent(i))/fFiMin[i]; | |
1066 | fHistoPtMinMothers->Fill(lowedge,fMinfill); | |
1067 | fHistoPtMinMothers->SetBinError(i,erPtMinStat[i]); | |
1068 | } | |
1069 | } | |
1070 | return kTRUE; | |
1071 | } | |
1072 | ||
1073 | //______________________________________________________________________________________ | |
1074 | void AliPtMothFromPtDaugh::WritePtMothHistoToFile(char *fileOutName) | |
1075 | { | |
1076 | // Write pt and ptMin histograms of mothers in a file | |
1077 | // with name fileOutName. Default name is "Mothers.root". | |
1078 | AliError(Form("Write mothers histograms in the file %s \n",fileOutName)); | |
1079 | if(!fHistoPtMothers) {AliError("Cannot write pt-mothers histogram! It doesn't exists!"); return;} | |
1080 | if(!fHistoPtMinMothers) { AliError("Cannot write ptMin-mothers histogram! It doesn't exists!"); return;} | |
1081 | TFile *outFile = TFile::Open(fileOutName,"RECREATE"); | |
1082 | outFile->cd(); | |
1083 | fHistoPtMothers->Write(); | |
1084 | fHistoPtMinMothers->Write(); | |
1085 | outFile->Close(); | |
1086 | return; | |
1087 | } |