]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/base/AliPtMothFromPtDaugh.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moving CounterCollection class to PWG3/base (Gines)
[u/mrichter/AliRoot.git] / PWG3 / base / AliPtMothFromPtDaugh.cxx
CommitLineData
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
50ClassImp(AliPtMothFromPtDaugh)
51//________________________________________________________________
52AliPtMothFromPtDaugh::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//________________________________________________________________
77AliPtMothFromPtDaugh::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//________________________________________________________________
102AliPtMothFromPtDaugh::~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//______________________________________________________________________________________
119AliPtMothFromPtDaugh::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//______________________________________________________________________________________
181AliPtMothFromPtDaugh& 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//______________________________________________________________________________________
192Bool_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//______________________________________________________________________________________
234void 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//______________________________________________________________________________________
256Bool_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//______________________________________________________________________________________
266Double_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//______________________________________________________________________________________
283void 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//______________________________________________________________________________________
297void 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//______________________________________________________________________________________
309void 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//______________________________________________________________________________________
321void 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//______________________________________________________________________________________
335void 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//______________________________________________________________________________________
344void 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//______________________________________________________________________________________
359void 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//______________________________________________________________________________________
374void 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//______________________________________________________________________________________
407Double_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//______________________________________________________________________________________
417Bool_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//______________________________________________________________________________________
435Bool_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//______________________________________________________________________________________
453Bool_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//______________________________________________________________________________________
471Bool_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//______________________________________________________________________________________
490void 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//______________________________________________________________________________________
521Int_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//______________________________________________________________________________________
532Double_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//______________________________________________________________________________________
550Double_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//______________________________________________________________________________________
569Double_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//______________________________________________________________________________________
583Double_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//______________________________________________________________________________________
598Double_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//______________________________________________________________________________________
615Double_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//______________________________________________________________________________________
634Double_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//______________________________________________________________________________________
648Double_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//______________________________________________________________________________________
663Bool_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//______________________________________________________________________________________
674Bool_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//______________________________________________________________________________________
693Bool_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//______________________________________________________________________________________
704Int_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//______________________________________________________________________________________
722Bool_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//______________________________________________________________________________________
734Bool_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//______________________________________________________________________________________
833Bool_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//______________________________________________________________________________________
925Bool_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//______________________________________________________________________________________
955Bool_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//______________________________________________________________________________________
989Bool_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//______________________________________________________________________________________
1026Bool_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//______________________________________________________________________________________
1074void 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 }