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