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