PWG3/base/AliPtMothFromPtDaugh.cxx

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