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