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   delete fmcArray;
 213   fmcArray = new TClonesArray(*(c.fmcArray));
 214   delete fRecoCandidate;
 215   fRecoCandidate = new AliAODRecoDecayHF(*(c.fRecoCandidate));
 216   delete fmcPartCandidate;
 217   fmcPartCandidate = new AliAODMCParticle(*(c.fmcPartCandidate));
 218   if (fProngs > 0){
 219     delete [] fLabelArray;
 220     delete [] fPtAccCut;
 221     delete [] fEtaAccCut;
 222     fLabelArray = new Int_t[fProngs];
 223     fPtAccCut = new Float_t[fProngs];
 224     fEtaAccCut = new Float_t[fProngs];
 225     if (c.fLabelArray) memcpy(fLabelArray,c.fLabelArray,fProngs*sizeof(Int_t));
 226     if (c.fPtAccCut) memcpy(fPtAccCut,c.fPtAccCut,fProngs*sizeof(Int_t));
 227     if (c.fEtaAccCut) memcpy(fEtaAccCut,c.fEtaAccCut,fProngs*sizeof(Int_t));
 228   }
 229 }
 230
 231 //___________________________________________________________
 232 void AliCFVertexingHF::SetDselection(UShort_t originDselection)
 233 {
 234         // setting the way the D0 will be selected
 235         // 0 --> only from c quarks
 236         // 1 --> only from b quarks
 237         // 2 --> from both c quarks and b quarks
 238
 239         fOriginDselection = originDselection;
 240
 241         if (fOriginDselection == 0) {
 242                 fKeepDfromB = kFALSE;
 243                 fKeepDfromBOnly = kFALSE;
 244         }
 245
 246         if (fOriginDselection == 1) {
 247                 fKeepDfromB = kTRUE;
 248                 fKeepDfromBOnly = kTRUE;
 249         }
 250
 251         if (fOriginDselection == 2) {
 252                 fKeepDfromB = kTRUE;
 253                 fKeepDfromBOnly = kFALSE;
 254         }
 255
 256         return;
 257 }
 258
 259 //______________________________________________________
 260 void AliCFVertexingHF::SetMCCandidateParam(Int_t label)
 261 {
 262         //
 263         // setting the parameters (candidate and n. daughters)
 264         //
 265
 266         fmcPartCandidate = dynamic_cast <AliAODMCParticle*> (fmcArray->At(label));
 267         if (fmcPartCandidate){
 268           fNDaughters = fmcPartCandidate->GetNDaughters();
 269         }
 270         else {
 271           AliError(Form("Dynamic cast failed, fNdaughters will remain set to %d",fNDaughters));
 272         }
 273         return;
 274 }
 275
 276 //____________________________________________________________
 277 Int_t AliCFVertexingHF::MCcquarkCounting(AliAODMCParticle* mcPart) const
 278 {
 279         //
 280         // counting the c-quarks
 281         //
 282
 283         Int_t cquarks = 0;
 284         if (mcPart) {
 285           if (mcPart->GetPdgCode() == 4) cquarks++;
 286           if (mcPart->GetPdgCode() == -4) cquarks++;
 287         }
 288         else {
 289                 AliWarning("Particle not found in tree, skipping\n");
 290                 return cquarks;
 291         }
 292
 293         return cquarks;
 294 }
 295
 296 //________________________________________________________
 297 Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClonesArray */*mcArray*/) const
 298 {
 299         //
 300         //checking the family
 301         //
 302
 303         Int_t pdgGranma = CheckOrigin();
 304
 305         if (pdgGranma == -99999){
 306                 AliDebug(2,"This particle does not have a quark in his genealogy\n");
 307                 return kFALSE;
 308         }
 309         if (pdgGranma == -9999){
 310                 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");
 311                 return kFALSE;
 312         }
 313
 314         if (pdgGranma == -999){
 315                 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");
 316                 return kFALSE;
 317         }
 318
 319         if (!CheckMCDaughters()) {
 320           AliDebug(3, "CheckMCDaughters false");
 321           return kFALSE;
 322         }
 323         if (!CheckMCChannelDecay()) {
 324                 AliDebug(3,"CheckMCChannelDecay false");
 325                 return kFALSE;
 326         }
 327         return kTRUE;
 328 }
 329
 330 //_________________________________________________________________________________________________
 331 Int_t AliCFVertexingHF::CheckOrigin() const
 332 {
 333         //
 334         // checking whether the mother of the particles come from a charm or a bottom quark
 335         //
 336
 337         Int_t pdgGranma = 0;
 338         Int_t mother = 0;
 339         mother = fmcPartCandidate->GetMother();
 340         Int_t istep = 0;
 341         Int_t abspdgGranma =0;
 342         Bool_t isFromB=kFALSE;
 343         Bool_t isQuarkFound=kFALSE;
 344         while (mother >=0 ){
 345                 istep++;
 346                 AliDebug(2,Form("mother at step %d = %d", istep, mother));
 347                 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
 348                 if (mcGranma){
 349                         pdgGranma = mcGranma->GetPdgCode();
 350                         AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
 351                         abspdgGranma = TMath::Abs(pdgGranma);
 352                         if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 353                           isFromB=kTRUE;
 354                         }
 355                         if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 356                         mother = mcGranma->GetMother();
 357                 }else{
 358                         AliError("Failed casting the mother particle!");
 359                         break;
 360                 }
 361         }
 362
 363         if(fRejectIfNoQuark && !isQuarkFound) return -99999;
 364         if(isFromB){
 365           if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
 366         }
 367         else{
 368           if (fKeepDfromBOnly) return -999;
 369         }
 370         return pdgGranma;
 371 }
 372
 373 //___________________________________________
 374 Bool_t AliCFVertexingHF::CheckMCDaughters()const
 375 {
 376         //
 377         // checking the daughters
 378         // at MC level
 379
 380         AliAODMCParticle *mcPartDaughter;
 381         Bool_t checkDaughters = kFALSE;
 382
 383         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 384         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 385         AliDebug(3,Form("label0 = %d, label1 = %d",label0,label1));
 386         if (label1<=0 || label0 <= 0){
 387                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 388                 return checkDaughters;
 389         }
 390
 391         if (fLabelArray == 0x0) {
 392                 return checkDaughters;
 393         }
 394
 395         for (Int_t iProng = 0; iProng<fProngs; iProng++){
 396                 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
 397                 if (!mcPartDaughter) {
 398                         AliWarning("At least one Daughter Particle not found in tree, skipping");
 399                         return checkDaughters;
 400                 }
 401         }
 402
 403         checkDaughters = kTRUE;
 404         return checkDaughters;
 405 }
 406
 407 //______________________________________________________
 408 Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
 409 {
 410         //
 411         // fill the container for Generator level selection
 412         //
 413
 414         Bool_t mcContainerFilled = kFALSE;
 415
 416         Double_t* vectorMC = new Double_t[fNVar];
 417         for (Int_t iVar = 0; iVar<fNVar; iVar++) vectorMC[iVar]= 9999.;
 418
 419         if (GetGeneratedValuesFromMCParticle(&vectorMC[0])){
 420                 for (Int_t iVar = 0; iVar<fNVar; iVar++){
 421                         containerInputMC[iVar] = vectorMC[iVar];
 422                 }
 423                 mcContainerFilled = kTRUE;
 424         }
 425         delete [] vectorMC;
 426         vectorMC = 0x0;
 427         return mcContainerFilled;
 428 }
 429
 430 //______________________________________________________
 431 Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput)
 432 {
 433         //
 434         // fill the container for Reconstrucred level selection
 435         //
 436
 437         Bool_t recoContainerFilled = kFALSE;
 438         Double_t* vectorValues = new Double_t[fNVar];
 439         Double_t* vectorReco = new Double_t[fNVar];
 440         for (Int_t iVar = 0; iVar<fNVar; iVar++) {
 441
 442                 vectorValues[iVar]= 9999.;
 443                 vectorReco[iVar]=9999.;
 444         }
 445
 446         if(fFillFromGenerated){
 447                 //filled with MC values
 448                 if (GetGeneratedValuesFromMCParticle(&vectorValues[0])){
 449                         for (Int_t iVar = 0; iVar<fNVar; iVar++){
 450                                 containerInput[iVar] = vectorValues[iVar];
 451                         }
 452                         recoContainerFilled = kTRUE;
 453                 }
 454         }
 455         else{
 456                 //filled with Reco values
 457
 458                 if (GetRecoValuesFromCandidate(&vectorReco[0])){
 459                         for (Int_t iVar = 0; iVar<fNVar; iVar++){
 460                           containerInput[iVar] = vectorReco[iVar];
 461                         }
 462                         recoContainerFilled = kTRUE;
 463                 }
 464         }
 465
 466         delete [] vectorValues;
 467         delete [] vectorReco;
 468         vectorValues = 0x0;
 469         vectorReco = 0x0;
 470
 471         return recoContainerFilled;
 472 }
 473
 474 //_____________________________________________________
 475 Bool_t AliCFVertexingHF::MCAcceptanceStep() const
 476 {
 477         //
 478         // checking the MC acceptance step
 479         //
 480
 481         Bool_t bMCAccStep = kFALSE;
 482
 483         AliAODMCParticle *mcPartDaughter;
 484         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 485         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 486         if (label1<=0 || label0 <= 0){
 487                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 488                 return bMCAccStep;
 489         }
 490
 491         if (fLabelArray == 0x0) {
 492                 return bMCAccStep;
 493         }
 494
 495         for (Int_t iProng = 0; iProng<fProngs; iProng++){
 496                 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
 497                 if (!mcPartDaughter) {
 498                         AliWarning("At least one Daughter Particle not found in tree, skipping");
 499                         return bMCAccStep;
 500                 }
 501                 Double_t eta = mcPartDaughter->Eta();
 502                 Double_t pt = mcPartDaughter->Pt();
 503
 504                 //set values of eta and pt in the constructor.
 505                 //              if (TMath::Abs(eta) > 0.9 || pt < 0.1){
 506                 if (TMath::Abs(eta) > fEtaAccCut[iProng] || pt < fPtAccCut[iProng]){
 507                         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]));
 508                         return bMCAccStep;
 509                 }
 510         }
 511         bMCAccStep = kTRUE;
 512         return bMCAccStep;
 513
 514 }
 515  //_____________________________________________________
 516 Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts **trackCuts) const
 517 {
 518         //
 519         // check on the kTPCrefit and kITSrefit conditions of the daughters
 520         //
 521         Bool_t bRefitStep = kFALSE;
 522
 523         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 524         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 525
 526         if (label1<=0 || label0 <= 0){
 527                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 528                 return bRefitStep;
 529         }
 530
 531         if (fLabelArray == 0x0) {
 532                 return bRefitStep;
 533         }
 534
 535         Int_t foundDaughters = 0;
 536         Int_t* temp = new Int_t[fProngs];
 537         for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
 538                 temp[ilabel] = fLabelArray[ilabel];
 539         }
 540
 541         //      if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
 542
 543         for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
 544                 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
 545                 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
 546                 Bool_t foundTrack = kFALSE;
 547                 Int_t prongindex;
 548                 for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
 549                         AliDebug(3,Form("fLabelArray[%d] = %d, track->GetLabel() = %d",ilabel,fLabelArray[ilabel],TMath::Abs(track->GetLabel())));
 550                         if ((track->GetLabel()<0)&&(fFakeSelection==1)) continue;
 551                         if ((track->GetLabel()>0)&&(fFakeSelection==2)) continue;
 552
 553                         if (TMath::Abs(track->GetLabel()) == temp[ilabel]) {
 554                                 foundTrack = kTRUE;
 555                                 temp[ilabel] = 0;
 556                                 prongindex=ilabel;
 557                                 break;
 558                         }
 559                 }
 560                 if (foundTrack){
 561                         foundDaughters++;
 562                         AliDebug(4,Form("daughter %d \n",foundDaughters));
 563                         if(trackCuts[prongindex]->GetRequireTPCRefit()){
 564                                 if(track->GetStatus()&AliESDtrack::kTPCrefit) {
 565                                         bRefitStep = kTRUE;
 566                                 }
 567                                 else {
 568                                         AliDebug(3, "Refit cut not passed , missing TPC refit\n");
 569                                         delete [] temp;
 570                                         temp = 0x0;
 571                                         return kFALSE;
 572                                 }
 573                         }
 574
 575                         if (trackCuts[prongindex]->GetRequireITSRefit()) {
 576                                 if(track->GetStatus()&AliESDtrack::kITSrefit){
 577                                         bRefitStep = kTRUE;
 578                                 }
 579                                 else {
 580                                         AliDebug(3, "Refit cut not passed , missing ITS refit\n");
 581                                         delete [] temp;
 582                                         temp = 0x0;
 583                                         return kFALSE;
 584                                 }
 585                         }
 586                 }
 587                 if (foundDaughters == fProngs){
 588                         break;
 589                 }
 590         }
 591         //}
 592         delete [] temp;
 593         temp = 0x0;
 594         if (foundDaughters== fProngs)  return bRefitStep;
 595         else return kFALSE;
 596 }
 597
 598 //____________________________________________________________________________
 599
 600 Bool_t AliCFVertexingHF::RecoStep()
 601 {
 602         //
 603         //check also vertex and ITS Refit and TPC Refit
 604         //
 605
 606         Bool_t bRecoStep = kFALSE;
 607         Int_t mcLabel = GetMCLabel();
 608
 609         if (mcLabel == -1) {
 610                 AliDebug(2,"No MC particle found");
 611                 return bRecoStep;
 612         }
 613         else{
 614                 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
 615                 if (!fmcPartCandidate){
 616                         AliWarning("Could not find associated MC in AOD MC tree");
 617                         return bRecoStep;
 618                 }
 619         }
 620
 621         Int_t pdgGranma = CheckOrigin();
 622
 623         if (pdgGranma == -99999){
 624                 AliDebug(2,"This particle does not have a quark in his genealogy\n");
 625                 return bRecoStep;
 626         }
 627         if (pdgGranma == -9999){
 628                 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");
 629                 return bRecoStep;
 630         }
 631
 632         if (pdgGranma == -999){
 633                 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");
 634                 return bRecoStep;
 635         }
 636
 637         bRecoStep=kTRUE;
 638         return bRecoStep;
 639 }
 640 //____________________________________________
 641 Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const
 642 {
 643         //
 644         // getting eta of the prong
 645         //
 646
 647         if (fRecoCandidate){
 648                 Double_t etaProng = fRecoCandidate->EtaProng(iProng);
 649                 return etaProng;
 650         }
 651         return 999999;
 652 }
 653 //______________________________________________________
 654 Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const
 655 {
 656         //
 657         // getting pt of the prong
 658         //
 659
 660         if (fRecoCandidate){
 661                 Double_t ptProng = fRecoCandidate->PtProng(iProng);
 662                 return ptProng;
 663         }
 664         return 999999;
 665
 666 }
 667
 668 //____________________________________________________________________
 669
 670 Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts **trackCuts) const
 671 {
 672         //
 673         // reco Acceptance step
 674         //
 675
 676         Bool_t bRecoAccStep = kFALSE;
 677
 678         Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
 679
 680         Float_t etaProng=0., ptProng=0.;
 681
 682         for (Int_t iProng =0; iProng<fProngs; iProng++){
 683
 684                 trackCuts[iProng]->GetEtaRange(etaCutMin, etaCutMax);
 685                 trackCuts[iProng]->GetPtRange(ptCutMin, ptCutMax);
 686                 etaProng = GetEtaProng(iProng);
 687                 ptProng = GetPtProng(iProng);
 688
 689                 Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
 690                 if (!acceptanceProng) {
 691                         AliDebug(2,"At least one reconstructed prong isn't in the acceptance\n");
 692                         return bRecoAccStep;
 693                 }
 694         }
 695
 696         bRecoAccStep=kTRUE;
 697         return bRecoAccStep;
 698 }
 699 //___________________________________________________________
 700
 701 Bool_t AliCFVertexingHF::FillUnfoldingMatrix(UInt_t pdg, Double_t fill[4]) const
 702 {
 703         //
 704         // filling the unfolding matrix
 705         //
 706
 707         if(fmcPartCandidate){
 708
 709                 fill[0] = GetPtCand();
 710                 fill[1] = GetYCand(pdg);
 711
 712                 fill[2] =  fmcPartCandidate->Pt();
 713                 fill[3] =  fmcPartCandidate->Y();
 714
 715                 return kTRUE;
 716         }
 717
 718         return kFALSE;
 719 }
 720 //___________________________________________________________
 721
 722 Int_t AliCFVertexingHF::CheckReflexion(Char_t isSign)
 723 {
 724         //
 725         // check for reflexion (particle/antiparticle)
 726         //
 727
 728         Int_t mcLabel = GetMCLabel();
 729
 730         if (mcLabel == -1) {
 731                 AliDebug(2,"No MC particle found");
 732                 return 0;
 733         }
 734         else{
 735                 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
 736                 if (!fmcPartCandidate){
 737                         AliWarning("Could not find associated MC in AOD MC tree");
 738                         return 0;
 739                 }
 740         }
 741
 742         if(fmcPartCandidate->GetPdgCode()>0) {
 743                 if (isSign == 1){ // I ask for antiparticle only
 744                         AliDebug(2,"candidate is particle, I ask for antiparticle only");
 745                         return 0;
 746                 }
 747                 return 1;  // particle
 748         }
 749         else if(fmcPartCandidate->GetPdgCode()<0) {
 750                 if (isSign == 0){ // I ask for particle only
 751                         AliDebug(2,"candidate is antiparticle, I ask for particle only");
 752                         return 0;
 753                 }
 754                 return 2;  // antiparticle
 755         }
 756         else return 0;  // ....shouldn't be...
 757
 758 }
 759 //___________________________________________________________
 760
 761 Bool_t AliCFVertexingHF::SetLabelArray()
 762 {
 763         //
 764         // setting the label arrays
 765         //
 766
 767         Bool_t bLabelArray = kFALSE;
 768
 769         fLabelArray = new Int_t[fProngs];
 770
 771         AliAODMCParticle *mcPartDaughter;
 772         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 773         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 774         AliDebug(2,Form("label0 = %d, label1 = %d",label0,label1));
 775         if (label1<=0 || label0 <= 0){
 776                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 777                 delete [] fLabelArray;
 778                 fLabelArray = 0x0;
 779                 return bLabelArray;
 780         }
 781
 782         if (label1 - label0 == fProngs-1){
 783                 for (Int_t iProng = 0; iProng<fProngs; iProng++){
 784                         mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
 785                         if (mcPartDaughter){
 786                                 fLabelArray[iProng] =  mcPartDaughter->GetLabel();
 787                         }
 788                         else{
 789                                 AliError("Failed casting the daughter particle, returning a NULL label array");
 790                                 delete [] fLabelArray;
 791                                 fLabelArray = 0x0;
 792                                 return bLabelArray;
 793                         }
 794                 }
 795
 796         }
 797         // resonant decay channel
 798         else if (label1 - label0 == fProngs-2 && fProngs > 2){
 799                 Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 800                 Int_t foundDaughters = 0;
 801                 for(Int_t iDau=0; iDau<fProngs-1; iDau++){
 802                         Int_t iLabelDau = labelFirstDau+iDau;
 803                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDau));
 804                         if ( ! part ) {
 805                           AliError("Wrong particle type in fmcArray");
 806                           delete [] fLabelArray;
 807                           fLabelArray = 0x0;
 808                           return bLabelArray;
 809                         }
 810                         Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 811                         if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 812                                 if (part) {
 813                                         fLabelArray[foundDaughters] = part->GetLabel();
 814                                         foundDaughters++;
 815                                 }
 816                                 else{
 817                                         AliError("Error while casting particle! returning a NULL array");
 818                                         delete [] fLabelArray;
 819                                         fLabelArray = 0x0;
 820                                         return bLabelArray;
 821                                 }
 822                         }
 823                         // added K0S case - Start
 824                         else if (pdgCode==311) {
 825                           if (part->GetNDaughters()!=1) {
 826                             delete [] fLabelArray;
 827                             fLabelArray = 0x0;
 828                             return bLabelArray;
 829                           }
 830                           Int_t labelK0Dau = part->GetDaughter(0);
 831                           AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelK0Dau));
 832                           if(!partK0S){
 833                             AliError("Error while casting particle! returning a NULL array");
 834                             delete [] fLabelArray;
 835                             fLabelArray = 0x0;
 836                             return bLabelArray;
 837                           }
 838                           Int_t nDauRes=partK0S->GetNDaughters();
 839                           if(nDauRes!=2 || partK0S->GetPdgCode()!=310) {
 840                             AliDebug(2,"No K0S on no 2-body decay");
 841                             delete [] fLabelArray;
 842                             fLabelArray = 0x0;
 843                             return bLabelArray;
 844                           }
 845                           Int_t labelFirstDauRes = partK0S->GetDaughter(0);
 846                           AliDebug(2,Form(" Found K0S (%d)",labelK0Dau));
 847                           for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 848                             Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 849                             AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 850                             if (dauRes){
 851                               if (TMath::Abs(dauRes->GetPdgCode())!=211) {
 852                                 AliDebug(2,"K0S doesn't decay in 2 charged pions!");
 853                                 delete [] fLabelArray;
 854                                 fLabelArray = 0x0;
 855                                 return bLabelArray;
 856                               }
 857                               else {
 858                                 fLabelArray[foundDaughters] = dauRes->GetLabel();
 859                                 foundDaughters++;
 860                               }
 861                             }
 862                             else {
 863                               AliError("Error while casting resonant daughter! returning a NULL array");
 864                               delete [] fLabelArray;
 865                               fLabelArray = 0x0;
 866                               return bLabelArray;
 867                             }
 868                           }
 869                         }
 870                         // added K0S case - End
 871                         else{
 872                                 Int_t nDauRes=part->GetNDaughters();
 873                                 if(nDauRes!=2) {
 874                                         delete [] fLabelArray;
 875                                         fLabelArray = 0x0;
 876                                         return bLabelArray;
 877                                 }
 878                                 Int_t labelFirstDauRes = part->GetDaughter(0);
 879                                 for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 880                                         Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 881                                         AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 882                                         if (dauRes){
 883                                                 fLabelArray[foundDaughters] = dauRes->GetLabel();
 884                                                 foundDaughters++;
 885                                         }
 886                                         else{
 887                                                 AliError("Error while casting resonant daughter! returning a NULL array");
 888                                                 delete [] fLabelArray;
 889                                                 fLabelArray = 0x0;
 890                                                 return bLabelArray;
 891                                         }
 892                                 }
 893                         }
 894                 }
 895                 if (foundDaughters != fProngs){
 896                         delete [] fLabelArray;
 897                         fLabelArray = 0x0;
 898                         return bLabelArray;
 899                 }
 900         }
 901         // wrong correspondance label <--> prongs
 902         else{
 903                 delete [] fLabelArray;
 904                 fLabelArray = 0x0;
 905                 return bLabelArray;
 906         }
 907         SetAccCut();   // setting the pt and eta acceptance cuts
 908         bLabelArray = kTRUE;
 909         return bLabelArray;
 910 }
 911
 912 //___________________________________________________________
 913
 914 void AliCFVertexingHF::SetPtAccCut(Float_t* ptAccCut)
 915 {
 916         //
 917         // setting the pt cut to be used in the Acceptance steps (MC+Reco)
 918         //
 919
 920         if (fProngs>0){
 921                 for (Int_t iP=0; iP<fProngs; iP++){
 922                         fPtAccCut[iP]=ptAccCut[iP];
 923                 }
 924         }
 925         return;
 926 }
 927
 928
 929
 930 //___________________________________________________________
 931
 932 void AliCFVertexingHF::SetEtaAccCut(Float_t* etaAccCut)
 933 {
 934         //
 935         // setting the eta cut to be used in the Acceptance steps (MC+Reco)
 936         //
 937
 938         if (fProngs>0){
 939                 for (Int_t iP=0; iP<fProngs; iP++){
 940                         fEtaAccCut[iP]=etaAccCut[iP];
 941                 }
 942         }
 943         return;
 944 }
 945 //___________________________________________________________
 946
 947 void AliCFVertexingHF::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
 948 {
 949         //
 950         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 951         //
 952
 953         if (fProngs>0){
 954                 for (Int_t iP=0; iP<fProngs; iP++){
 955                         fPtAccCut[iP]=ptAccCut[iP];
 956                         fEtaAccCut[iP]=etaAccCut[iP];
 957                 }
 958         }
 959         return;
 960 }
 961
 962 //___________________________________________________________
 963
 964 void AliCFVertexingHF::SetAccCut()
 965 {
 966         //
 967         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 968         //
 969
 970         if (fProngs>0){
 971                 for (Int_t iP=0; iP<fProngs; iP++){
 972                         fPtAccCut[iP]=0.1;
 973                         fEtaAccCut[iP]=0.9;
 974                 }
 975         }
 976         return;
 977 }