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