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