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   }