PWG3/vertexingHF/AliCFVertexingHF.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 for HF corrections as a function of many variables and step
  20 // Author : C. Zampolli, CERN
  21 // D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
  22 // Base class for HF Unfolding - agrelli@uu.nl
  23 //-----------------------------------------------------------------------
  24
  25 #include "TParticle.h"
  26 #include "TClonesArray.h"
  27 #include "AliAODMCParticle.h"
  28 #include "AliAODRecoDecayHF.h"
  29 #include "AliAODRecoDecayHF2Prong.h"
  30 #include "AliAODRecoDecayHF3Prong.h"
  31 #include "AliAODRecoDecayHF4Prong.h"
  32 #include "AliAODRecoCascadeHF.h"
  33 #include "AliAODMCHeader.h"
  34 #include "AliAODEvent.h"
  35 #include "AliLog.h"
  36 #include "AliESDtrackCuts.h"
  37 #include "AliESDtrack.h"
  38
  39 #include "AliCFVertexingHF.h"
  40
  41 //___________________________________________________________
  42 AliCFVertexingHF::AliCFVertexingHF() :
  43         fmcArray(0x0),
  44         fRecoCandidate(0),
  45         fmcPartCandidate(0x0),
  46         fNDaughters(0),
  47         fNVar(0),
  48         fzPrimVertex(0),
  49         fzMCVertex(0),
  50         fFillFromGenerated(0),
  51         fOriginDselection(0),
  52         fKeepDfromB(kFALSE),
  53         fKeepDfromBOnly(kFALSE),
  54         fmcLabel(0),
  55         fProngs(-1),
  56         fLabelArray(0x0),
  57         fCentValue(0.),
  58         fPtAccCut(0x0),
  59         fEtaAccCut(0x0),
  60         fFakeSelection(0),
  61         fFake(1.), // setting to MC value
  62         fRejectIfNoQuark(kFALSE)
  63 {
  64         //
  65         // constructor
  66         //
  67
  68
  69         return;
  70 }
  71
  72
  73
  74 //_____________________________________________________
  75 AliCFVertexingHF::AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection) :
  76         fmcArray(mcArray),
  77         fRecoCandidate(0),
  78         fmcPartCandidate(0x0),
  79         fNDaughters(0),
  80         fNVar(0),
  81         fzPrimVertex(0),
  82         fzMCVertex(0),
  83         fFillFromGenerated(0),
  84         fOriginDselection(0),
  85         fKeepDfromB(kFALSE),
  86         fKeepDfromBOnly(kFALSE),
  87         fmcLabel(0),
  88         fProngs(-1),
  89         fLabelArray(0x0),
  90         fCentValue(0.),
  91         fPtAccCut(0x0),
  92         fEtaAccCut(0x0),
  93         fFakeSelection(0),
  94         fFake(1.), // setting to MC value
  95         fRejectIfNoQuark(kFALSE)
  96 {
  97         //
  98         // constructor with mcArray
  99         //
 100
 101         SetDselection(originDselection);
 102         return;
 103 }
 104
 105 //_______________________________________________________
 106 AliCFVertexingHF::~AliCFVertexingHF()
 107 {
 108         //
 109         // destructor
 110         //
 111
 112         if (fmcArray) fmcArray = 0x0;
 113         if (fRecoCandidate) fRecoCandidate = 0x0;
 114         if (fmcPartCandidate) fmcPartCandidate = 0x0;
 115         if (fLabelArray){
 116                 delete [] fLabelArray;
 117                 fLabelArray = 0x0;
 118         }
 119         if (fPtAccCut){
 120                 delete [] fPtAccCut;
 121                 fPtAccCut = 0x0;
 122         }
 123         if (fEtaAccCut){
 124                 delete [] fEtaAccCut;
 125                 fEtaAccCut = 0x0;
 126         }
 127 }
 128
 129 //_____________________________________________________
 130 AliCFVertexingHF& AliCFVertexingHF::operator=(const AliCFVertexingHF& c)
 131 {
 132         //
 133         // assigment operator
 134         //
 135
 136         if (this!= &c){
 137                 TObject::operator=(c);
 138                 fmcArray = c.fmcArray;
 139                 fRecoCandidate = c.fRecoCandidate;
 140                 fmcPartCandidate = c.fmcPartCandidate;
 141                 fNDaughters = c.fNDaughters;
 142                 fNVar = c.fNVar;
 143                 fzPrimVertex = c.fzPrimVertex;
 144                 fzMCVertex = c.fzMCVertex;
 145                 fFillFromGenerated = c.fFillFromGenerated;
 146                 fOriginDselection = c.fOriginDselection;
 147                 fKeepDfromB = c.fKeepDfromB;
 148                 fKeepDfromBOnly = c.fKeepDfromBOnly;
 149                 fmcLabel = c.fmcLabel;
 150                 fProngs=c.fProngs;
 151                 fCentValue=c.fCentValue;
 152                 fFakeSelection=c.fFakeSelection;
 153                 fFake=c.fFake;
 154                 fRejectIfNoQuark=c.fRejectIfNoQuark;
 155                 if (fProngs > 0){
 156                         fLabelArray = new Int_t[fProngs];
 157                         fPtAccCut = new Float_t[fProngs];
 158                         fEtaAccCut = new Float_t[fProngs];
 159                         for(Int_t iP=0; iP<fProngs; iP++){
 160                                 fLabelArray[iP]=c.fLabelArray[iP];
 161                                 fPtAccCut[iP]=c.fPtAccCut[iP];
 162                                 fEtaAccCut[iP]=c.fEtaAccCut[iP];
 163                         }
 164                 }
 165         }
 166
 167         return *this;
 168 }
 169
 170 //____________________________________________________
 171 AliCFVertexingHF::AliCFVertexingHF(const AliCFVertexingHF &c) :
 172         TObject(c),
 173         fmcArray(c.fmcArray),
 174         fRecoCandidate(c.fRecoCandidate),
 175         fmcPartCandidate(c.fmcPartCandidate),
 176         fNDaughters(c.fNDaughters),
 177         fNVar(c.fNVar),
 178         fzPrimVertex(c.fzPrimVertex),
 179         fzMCVertex(c.fzMCVertex),
 180         fFillFromGenerated(c.fFillFromGenerated),
 181         fOriginDselection (c.fOriginDselection),
 182         fKeepDfromB (c.fKeepDfromB),
 183         fKeepDfromBOnly (c.fKeepDfromBOnly),
 184         fmcLabel(c.fmcLabel),
 185         fProngs(c.fProngs),
 186         fLabelArray(0x0),
 187         fCentValue(c.fCentValue),
 188         fPtAccCut(0x0),
 189         fEtaAccCut(0x0),
 190         fFakeSelection(c.fFakeSelection),
 191         fFake(c.fFake),
 192         fRejectIfNoQuark(c.fRejectIfNoQuark)
 193 {
 194         //
 195         //copy constructor
 196         //
 197         if (fProngs > 0){
 198                 fLabelArray = new Int_t[fProngs];
 199                 fPtAccCut = new Float_t[fProngs];
 200                 fEtaAccCut = new Float_t[fProngs];
 201                 if (c.fLabelArray) memcpy(fLabelArray,c.fLabelArray,fProngs*sizeof(Int_t));
 202                 if (c.fPtAccCut) memcpy(fPtAccCut,c.fPtAccCut,fProngs*sizeof(Int_t));
 203                 if (c.fEtaAccCut) memcpy(fEtaAccCut,c.fEtaAccCut,fProngs*sizeof(Int_t));
 204         }
 205 }
 206
 207 //___________________________________________________________
 208 void AliCFVertexingHF::SetDselection(UShort_t originDselection)
 209 {
 210         // setting the way the D0 will be selected
 211         // 0 --> only from c quarks
 212         // 1 --> only from b quarks
 213         // 2 --> from both c quarks and b quarks
 214
 215         fOriginDselection = originDselection;
 216
 217         if (fOriginDselection == 0) {
 218                 fKeepDfromB = kFALSE;
 219                 fKeepDfromBOnly = kFALSE;
 220         }
 221
 222         if (fOriginDselection == 1) {
 223                 fKeepDfromB = kTRUE;
 224                 fKeepDfromBOnly = kTRUE;
 225         }
 226
 227         if (fOriginDselection == 2) {
 228                 fKeepDfromB = kTRUE;
 229                 fKeepDfromBOnly = kFALSE;
 230         }
 231
 232         return;
 233 }
 234
 235 //______________________________________________________
 236 void AliCFVertexingHF::SetMCCandidateParam(Int_t label)
 237 {
 238         //
 239         // setting the parameters (candidate and n. daughters)
 240         //
 241
 242         fmcPartCandidate = dynamic_cast <AliAODMCParticle*> (fmcArray->At(label));
 243         if (fmcPartCandidate){
 244           fNDaughters = fmcPartCandidate->GetNDaughters();
 245         }
 246         else {
 247           AliError(Form("Dynamic cast failed, fNdaughters will remain set to %d",fNDaughters));
 248         }
 249         return;
 250 }
 251
 252 //____________________________________________________________
 253 Int_t AliCFVertexingHF::MCcquarkCounting(AliAODMCParticle* mcPart) const
 254 {
 255         //
 256         // counting the c-quarks
 257         //
 258
 259         Int_t cquarks = 0;
 260         if (mcPart) {
 261           if (mcPart->GetPdgCode() == 4) cquarks++;
 262           if (mcPart->GetPdgCode() == -4) cquarks++;
 263         }
 264         else {
 265                 AliWarning("Particle not found in tree, skipping\n");
 266                 return cquarks;
 267         }
 268
 269         return cquarks;
 270 }
 271
 272 //________________________________________________________
 273 Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClonesArray */*mcArray*/) const
 274 {
 275         //
 276         //checking the family
 277         //
 278
 279         Int_t pdgGranma = CheckOrigin();
 280
 281         if (pdgGranma == -99999){
 282                 AliDebug(2,"This particle does not have a quark in his genealogy\n");
 283                 return kFALSE;
 284         }
 285         if (pdgGranma == -9999){
 286                 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
 287                 return kFALSE;
 288         }
 289
 290         if (pdgGranma == -999){
 291                 AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
 292                 return kFALSE;
 293         }
 294
 295         if (!CheckMCDaughters()) {
 296           AliDebug(3, "CheckMCDaughters false");
 297           return kFALSE;
 298         }
 299         if (!CheckMCChannelDecay()) {
 300                 AliDebug(3,"CheckMCChannelDecay false");
 301                 return kFALSE;
 302         }
 303         return kTRUE;
 304 }
 305
 306 //_________________________________________________________________________________________________
 307 Int_t AliCFVertexingHF::CheckOrigin() const
 308 {
 309         //
 310         // checking whether the mother of the particles come from a charm or a bottom quark
 311         //
 312
 313         Int_t pdgGranma = 0;
 314         Int_t mother = 0;
 315         mother = fmcPartCandidate->GetMother();
 316         Int_t istep = 0;
 317         Int_t abspdgGranma =0;
 318         Bool_t isFromB=kFALSE;
 319         Bool_t isQuarkFound=kFALSE;
 320         while (mother >0 ){
 321                 istep++;
 322                 AliDebug(2,Form("mother at step %d = %d", istep, mother));
 323                 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
 324                 if (mcGranma){
 325                         pdgGranma = mcGranma->GetPdgCode();
 326                         AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
 327                         abspdgGranma = TMath::Abs(pdgGranma);
 328                         if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 329                           isFromB=kTRUE;
 330                         }
 331                         if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 332                         mother = mcGranma->GetMother();
 333                 }else{
 334                         AliError("Failed casting the mother particle!");
 335                         break;
 336                 }
 337         }
 338
 339         if(fRejectIfNoQuark && !isQuarkFound) return -99999;
 340         if(isFromB){
 341           if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
 342         }
 343         else{
 344           if (fKeepDfromBOnly) return -999;
 345         }
 346         return pdgGranma;
 347 }
 348
 349 //___________________________________________
 350 Bool_t AliCFVertexingHF::CheckMCDaughters()const
 351 {
 352         //
 353         // checking the daughters
 354         // at MC level
 355
 356         AliAODMCParticle *mcPartDaughter;
 357         Bool_t checkDaughters = kFALSE;
 358
 359         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 360         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 361         AliDebug(3,Form("label0 = %d, label1 = %d",label0,label1));
 362         if (label1==0 || label0 == 0){
 363                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 364                 return checkDaughters;
 365         }
 366
 367         if (fLabelArray == 0x0) {
 368                 return checkDaughters;
 369         }
 370
 371         for (Int_t iProng = 0; iProng<fProngs; iProng++){
 372                 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
 373                 if (!mcPartDaughter) {
 374                         AliWarning("At least one Daughter Particle not found in tree, skipping");
 375                         return checkDaughters;
 376                 }
 377         }
 378
 379         checkDaughters = kTRUE;
 380         return checkDaughters;
 381 }
 382
 383 //______________________________________________________
 384 Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
 385 {
 386         //
 387         // fill the container for Generator level selection
 388         //
 389
 390         Bool_t mcContainerFilled = kFALSE;
 391
 392         Double_t* vectorMC = new Double_t[fNVar];
 393         for (Int_t iVar = 0; iVar<fNVar; iVar++) vectorMC[iVar]= 9999.;
 394
 395         if (GetGeneratedValuesFromMCParticle(&vectorMC[0])){
 396                 for (Int_t iVar = 0; iVar<fNVar; iVar++){
 397                         containerInputMC[iVar] = vectorMC[iVar];
 398                 }
 399                 mcContainerFilled = kTRUE;
 400         }
 401         delete [] vectorMC;
 402         vectorMC = 0x0;
 403         return mcContainerFilled;
 404 }
 405
 406 //______________________________________________________
 407 Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput)
 408 {
 409         //
 410         // fill the container for Reconstrucred level selection
 411         //
 412
 413         Bool_t recoContainerFilled = kFALSE;
 414         Double_t* vectorValues = new Double_t[fNVar];
 415         Double_t* vectorReco = new Double_t[fNVar];
 416         for (Int_t iVar = 0; iVar<fNVar; iVar++) {
 417
 418                 vectorValues[iVar]= 9999.;
 419                 vectorReco[iVar]=9999.;
 420         }
 421
 422         if(fFillFromGenerated){
 423                 //filled with MC values
 424                 if (GetGeneratedValuesFromMCParticle(&vectorValues[0])){
 425                         for (Int_t iVar = 0; iVar<fNVar; iVar++){
 426                                 containerInput[iVar] = vectorValues[iVar];
 427                         }
 428                         recoContainerFilled = kTRUE;
 429                 }
 430         }
 431         else{
 432                 //filled with Reco values
 433
 434                 if (GetRecoValuesFromCandidate(&vectorReco[0])){
 435                         for (Int_t iVar = 0; iVar<fNVar; iVar++){
 436                           containerInput[iVar] = vectorReco[iVar];
 437                         }
 438                         recoContainerFilled = kTRUE;
 439                 }
 440         }
 441
 442         delete [] vectorValues;
 443         delete [] vectorReco;
 444         vectorValues = 0x0;
 445         vectorReco = 0x0;
 446         return recoContainerFilled;
 447 }
 448
 449 //_____________________________________________________
 450 Bool_t AliCFVertexingHF::MCAcceptanceStep() const
 451 {
 452         //
 453         // checking the MC acceptance step
 454         //
 455
 456         Bool_t bMCAccStep = kFALSE;
 457
 458         AliAODMCParticle *mcPartDaughter;
 459         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 460         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 461         if (label1==0 || label0 == 0){
 462                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 463                 return bMCAccStep;
 464         }
 465
 466         if (fLabelArray == 0x0) {
 467                 return bMCAccStep;
 468         }
 469
 470         for (Int_t iProng = 0; iProng<fProngs; iProng++){
 471                 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
 472                 if (!mcPartDaughter) {
 473                         AliWarning("At least one Daughter Particle not found in tree, skipping");
 474                         return bMCAccStep;
 475                 }
 476                 Double_t eta = mcPartDaughter->Eta();
 477                 Double_t pt = mcPartDaughter->Pt();
 478
 479                 //set values of eta and pt in the constructor.
 480                 //              if (TMath::Abs(eta) > 0.9 || pt < 0.1){
 481                 if (TMath::Abs(eta) > fEtaAccCut[iProng] || pt < fPtAccCut[iProng]){
 482                         AliDebug(3,Form("At least one daughter has eta or pt outside the required range (|eta| = %f, pt = %f, should be |eta| < %f, pt > %f \n", TMath::Abs(eta), pt, fEtaAccCut[iProng], fPtAccCut[iProng]));
 483                         return bMCAccStep;
 484                 }
 485         }
 486         bMCAccStep = kTRUE;
 487         return bMCAccStep;
 488
 489 }
 490  //_____________________________________________________
 491 Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts **trackCuts) const
 492 {
 493         //
 494         // check on the kTPCrefit and kITSrefit conditions of the daughters
 495         //
 496         Bool_t bRefitStep = kFALSE;
 497
 498         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 499         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 500
 501         if (label1==0 || label0 == 0){
 502                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 503                 return bRefitStep;
 504         }
 505
 506         if (fLabelArray == 0x0) {
 507                 return bRefitStep;
 508         }
 509
 510         Int_t foundDaughters = 0;
 511         Int_t* temp = new Int_t[fProngs];
 512         for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
 513                 temp[ilabel] = fLabelArray[ilabel];
 514         }
 515
 516         //      if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
 517
 518         for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
 519                 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
 520                 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
 521                 Bool_t foundTrack = kFALSE;
 522                 Int_t prongindex;
 523                 for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
 524                         AliDebug(3,Form("fLabelArray[%d] = %d, track->GetLabel() = %d",ilabel,fLabelArray[ilabel],TMath::Abs(track->GetLabel())));
 525                         if ((track->GetLabel()<0)&&(fFakeSelection==1)) continue;
 526                         if ((track->GetLabel()>0)&&(fFakeSelection==2)) continue;
 527
 528                         if (TMath::Abs(track->GetLabel()) == temp[ilabel]) {
 529                                 foundTrack = kTRUE;
 530                                 temp[ilabel] = 0;
 531                                 prongindex=ilabel;
 532                                 break;
 533                         }
 534                 }
 535                 if (foundTrack){
 536                         foundDaughters++;
 537                         AliDebug(4,Form("daughter %d \n",foundDaughters));
 538                         if(trackCuts[prongindex]->GetRequireTPCRefit()){
 539                                 if(track->GetStatus()&AliESDtrack::kTPCrefit) {
 540                                         bRefitStep = kTRUE;
 541                                 }
 542                                 else {
 543                                         AliDebug(3, "Refit cut not passed , missing TPC refit\n");
 544                                         delete [] temp;
 545                                         temp = 0x0;
 546                                         return kFALSE;
 547                                 }
 548                         }
 549
 550                         if (trackCuts[prongindex]->GetRequireITSRefit()) {
 551                                 if(track->GetStatus()&AliESDtrack::kITSrefit){
 552                                         bRefitStep = kTRUE;
 553                                 }
 554                                 else {
 555                                         AliDebug(3, "Refit cut not passed , missing ITS refit\n");
 556                                         delete [] temp;
 557                                         temp = 0x0;
 558                                         return kFALSE;
 559                                 }
 560                         }
 561                 }
 562                 if (foundDaughters == fProngs){
 563                         break;
 564                 }
 565         }
 566         //}
 567         delete [] temp;
 568         temp = 0x0;
 569         if (foundDaughters== fProngs)  return bRefitStep;
 570         else return kFALSE;
 571 }
 572
 573 //____________________________________________________________________________
 574
 575 Bool_t AliCFVertexingHF::RecoStep()
 576 {
 577         //
 578         //check also vertex and ITS Refit and TPC Refit
 579         //
 580
 581         Bool_t bRecoStep = kFALSE;
 582         Int_t mcLabel = GetMCLabel();
 583
 584         if (mcLabel == -1) {
 585                 AliDebug(2,"No MC particle found");
 586                 return bRecoStep;
 587         }
 588         else{
 589                 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
 590                 if (!fmcPartCandidate){
 591                         AliWarning("Could not find associated MC in AOD MC tree");
 592                         return bRecoStep;
 593                 }
 594         }
 595
 596         Int_t pdgGranma = CheckOrigin();
 597
 598         if (pdgGranma == -99999){
 599                 AliDebug(2,"This particle does not have a quark in his genealogy\n");
 600                 return bRecoStep;
 601         }
 602         if (pdgGranma == -9999){
 603                 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only prompt charm particles\n");
 604                 return bRecoStep;
 605         }
 606
 607         if (pdgGranma == -999){
 608                 AliDebug(2,"This particle come from a  prompt charm particle but according to the settings of the task, we want only the ones coming from B\n");
 609                 return bRecoStep;
 610         }
 611
 612         bRecoStep=kTRUE;
 613         return bRecoStep;
 614 }
 615 //____________________________________________
 616 Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const
 617 {
 618         //
 619         // getting eta of the prong
 620         //
 621
 622         if (fRecoCandidate){
 623                 Double_t etaProng = fRecoCandidate->EtaProng(iProng);
 624                 return etaProng;
 625         }
 626         return 999999;
 627 }
 628 //______________________________________________________
 629 Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const
 630 {
 631         //
 632         // getting pt of the prong
 633         //
 634
 635         if (fRecoCandidate){
 636                 Double_t ptProng = fRecoCandidate->PtProng(iProng);
 637                 return ptProng;
 638         }
 639         return 999999;
 640
 641 }
 642
 643 //____________________________________________________________________
 644
 645 Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts **trackCuts) const
 646 {
 647         //
 648         // reco Acceptance step
 649         //
 650
 651         Bool_t bRecoAccStep = kFALSE;
 652
 653         Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
 654
 655         Float_t etaProng=0., ptProng=0.;
 656
 657         for (Int_t iProng =0; iProng<fProngs; iProng++){
 658
 659                 trackCuts[iProng]->GetEtaRange(etaCutMin, etaCutMax);
 660                 trackCuts[iProng]->GetPtRange(ptCutMin, ptCutMax);
 661                 etaProng = GetEtaProng(iProng);
 662                 ptProng = GetPtProng(iProng);
 663
 664                 Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
 665                 if (!acceptanceProng) {
 666                         AliDebug(2,"At least one reconstructed prong isn't in the acceptance\n");
 667                         return bRecoAccStep;
 668                 }
 669         }
 670
 671         bRecoAccStep=kTRUE;
 672         return bRecoAccStep;
 673 }
 674 //___________________________________________________________
 675
 676 Bool_t AliCFVertexingHF::FillUnfoldingMatrix(Double_t fill[4]) const
 677 {
 678         //
 679         // filling the unfolding matrix
 680         //
 681
 682         if(fmcPartCandidate){
 683
 684                 fill[0] = GetPtCand();
 685                 fill[1] = GetYCand();
 686
 687                 fill[2] =  fmcPartCandidate->Pt();
 688                 fill[3] =  fmcPartCandidate->Y();
 689
 690                 return kTRUE;
 691         }
 692
 693         return kFALSE;
 694 }
 695 //___________________________________________________________
 696
 697 Int_t AliCFVertexingHF::CheckReflexion(Char_t isSign)
 698 {
 699         //
 700         // check for reflexion (particle/antiparticle)
 701         //
 702
 703         Int_t mcLabel = GetMCLabel();
 704
 705         if (mcLabel == -1) {
 706                 AliDebug(2,"No MC particle found");
 707                 return 0;
 708         }
 709         else{
 710                 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
 711                 if (!fmcPartCandidate){
 712                         AliWarning("Could not find associated MC in AOD MC tree");
 713                         return 0;
 714                 }
 715         }
 716
 717         if(fmcPartCandidate->GetPdgCode()>0) {
 718                 if (isSign == 1){ // I ask for antiparticle only
 719                         AliDebug(2,"candidate is particle, I ask for antiparticle only");
 720                         return 0;
 721                 }
 722                 return 1;  // particle
 723         }
 724         else if(fmcPartCandidate->GetPdgCode()<0) {
 725                 if (isSign == 0){ // I ask for particle only
 726                         AliDebug(2,"candidate is antiparticle, I ask for particle only");
 727                         return 0;
 728                 }
 729                 return 2;  // antiparticle
 730         }
 731         else return 0;  // ....shouldn't be...
 732
 733 }
 734 //___________________________________________________________
 735
 736 Bool_t AliCFVertexingHF::SetLabelArray()
 737 {
 738         //
 739         // setting the label arrays
 740         //
 741
 742         Bool_t bLabelArray = kFALSE;
 743
 744         fLabelArray = new Int_t[fProngs];
 745
 746         AliAODMCParticle *mcPartDaughter;
 747         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 748         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 749         AliDebug(2,Form("label0 = %d, label1 = %d",label0,label1));
 750         if (label1==0 || label0 == 0){
 751                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 752                 delete [] fLabelArray;
 753                 fLabelArray = 0x0;
 754                 return bLabelArray;
 755         }
 756
 757         if (label1 - label0 == fProngs-1){
 758                 for (Int_t iProng = 0; iProng<fProngs; iProng++){
 759                         mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
 760                         if (mcPartDaughter){
 761                                 fLabelArray[iProng] =  mcPartDaughter->GetLabel();
 762                         }
 763                         else{
 764                                 AliError("Failed casting the daughter particle, returning a NULL label array");
 765                                 delete [] fLabelArray;
 766                                 fLabelArray = 0x0;
 767                                 return bLabelArray;
 768                         }
 769                 }
 770
 771         }
 772         // resonant decay channel
 773         else if (label1 - label0 == fProngs-2 && fProngs > 2){
 774                 Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 775                 Int_t foundDaughters = 0;
 776                 for(Int_t iDau=0; iDau<fProngs-1; iDau++){
 777                         Int_t iLabelDau = labelFirstDau+iDau;
 778                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDau));
 779                         if ( ! part ) {
 780                           AliError("Wrong particle type in fmcArray");
 781                           delete [] fLabelArray;
 782                           fLabelArray = 0x0;
 783                           return bLabelArray;
 784                         }
 785                         Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 786                         if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 787                                 if (part) {
 788                                         fLabelArray[foundDaughters] = part->GetLabel();
 789                                         foundDaughters++;
 790                                 }
 791                                 else{
 792                                         AliError("Error while casting particle! returning a NULL array");
 793                                         delete [] fLabelArray;
 794                                         fLabelArray = 0x0;
 795                                         return bLabelArray;
 796                                 }
 797                         }
 798                         else{
 799                                 Int_t nDauRes=part->GetNDaughters();
 800                                 if(nDauRes!=2) {
 801                                         delete [] fLabelArray;
 802                                         fLabelArray = 0x0;
 803                                         return bLabelArray;
 804                                 }
 805                                 Int_t labelFirstDauRes = part->GetDaughter(0);
 806                                 for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 807                                         Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 808                                         AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 809                                         if (dauRes){
 810                                                 fLabelArray[foundDaughters] = dauRes->GetLabel();
 811                                                 foundDaughters++;
 812                                         }
 813                                         else{
 814                                                 AliError("Error while casting resonant daughter! returning a NULL array");
 815                                                 delete [] fLabelArray;
 816                                                 fLabelArray = 0x0;
 817                                                 return bLabelArray;
 818                                         }
 819                                 }
 820                         }
 821                 }
 822                 if (foundDaughters != fProngs){
 823                         delete [] fLabelArray;
 824                         fLabelArray = 0x0;
 825                         return bLabelArray;
 826                 }
 827         }
 828         // wrong correspondance label <--> prongs
 829         else{
 830                 delete [] fLabelArray;
 831                 fLabelArray = 0x0;
 832                 return bLabelArray;
 833         }
 834         SetAccCut();   // setting the pt and eta acceptance cuts
 835         bLabelArray = kTRUE;
 836         return bLabelArray;
 837 }
 838
 839 //___________________________________________________________
 840
 841 void AliCFVertexingHF::SetPtAccCut(Float_t* ptAccCut)
 842 {
 843         //
 844         // setting the pt cut to be used in the Acceptance steps (MC+Reco)
 845         //
 846
 847         if (fProngs>0){
 848                 for (Int_t iP=0; iP<fProngs; iP++){
 849                         fPtAccCut[iP]=ptAccCut[iP];
 850                 }
 851         }
 852         return;
 853 }
 854
 855
 856
 857 //___________________________________________________________
 858
 859 void AliCFVertexingHF::SetEtaAccCut(Float_t* etaAccCut)
 860 {
 861         //
 862         // setting the eta cut to be used in the Acceptance steps (MC+Reco)
 863         //
 864
 865         if (fProngs>0){
 866                 for (Int_t iP=0; iP<fProngs; iP++){
 867                         fEtaAccCut[iP]=etaAccCut[iP];
 868                 }
 869         }
 870         return;
 871 }
 872 //___________________________________________________________
 873
 874 void AliCFVertexingHF::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
 875 {
 876         //
 877         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 878         //
 879
 880         if (fProngs>0){
 881                 for (Int_t iP=0; iP<fProngs; iP++){
 882                         fPtAccCut[iP]=ptAccCut[iP];
 883                         fEtaAccCut[iP]=etaAccCut[iP];
 884                 }
 885         }
 886         return;
 887 }
 888
 889 //___________________________________________________________
 890
 891 void AliCFVertexingHF::SetAccCut()
 892 {
 893         //
 894         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 895         //
 896
 897         if (fProngs>0){
 898                 for (Int_t iP=0; iP<fProngs; iP++){
 899                         fPtAccCut[iP]=0.1;
 900                         fEtaAccCut[iP]=0.9;
 901                 }
 902         }
 903         return;
 904 }