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
 465         return recoContainerFilled;
 466 }
 467
 468 //_____________________________________________________
 469 Bool_t AliCFVertexingHF::MCAcceptanceStep() const
 470 {
 471         //
 472         // checking the MC acceptance step
 473         //
 474
 475         Bool_t bMCAccStep = kFALSE;
 476
 477         AliAODMCParticle *mcPartDaughter;
 478         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 479         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 480         if (label1==0 || label0 == 0){
 481                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 482                 return bMCAccStep;
 483         }
 484
 485         if (fLabelArray == 0x0) {
 486                 return bMCAccStep;
 487         }
 488
 489         for (Int_t iProng = 0; iProng<fProngs; iProng++){
 490                 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
 491                 if (!mcPartDaughter) {
 492                         AliWarning("At least one Daughter Particle not found in tree, skipping");
 493                         return bMCAccStep;
 494                 }
 495                 Double_t eta = mcPartDaughter->Eta();
 496                 Double_t pt = mcPartDaughter->Pt();
 497
 498                 //set values of eta and pt in the constructor.
 499                 //              if (TMath::Abs(eta) > 0.9 || pt < 0.1){
 500                 if (TMath::Abs(eta) > fEtaAccCut[iProng] || pt < fPtAccCut[iProng]){
 501                         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]));
 502                         return bMCAccStep;
 503                 }
 504         }
 505         bMCAccStep = kTRUE;
 506         return bMCAccStep;
 507
 508 }
 509  //_____________________________________________________
 510 Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts **trackCuts) const
 511 {
 512         //
 513         // check on the kTPCrefit and kITSrefit conditions of the daughters
 514         //
 515         Bool_t bRefitStep = kFALSE;
 516
 517         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 518         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 519
 520         if (label1==0 || label0 == 0){
 521                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 522                 return bRefitStep;
 523         }
 524
 525         if (fLabelArray == 0x0) {
 526                 return bRefitStep;
 527         }
 528
 529         Int_t foundDaughters = 0;
 530         Int_t* temp = new Int_t[fProngs];
 531         for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
 532                 temp[ilabel] = fLabelArray[ilabel];
 533         }
 534
 535         //      if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
 536
 537         for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
 538                 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
 539                 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
 540                 Bool_t foundTrack = kFALSE;
 541                 Int_t prongindex;
 542                 for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
 543                         AliDebug(3,Form("fLabelArray[%d] = %d, track->GetLabel() = %d",ilabel,fLabelArray[ilabel],TMath::Abs(track->GetLabel())));
 544                         if ((track->GetLabel()<0)&&(fFakeSelection==1)) continue;
 545                         if ((track->GetLabel()>0)&&(fFakeSelection==2)) continue;
 546
 547                         if (TMath::Abs(track->GetLabel()) == temp[ilabel]) {
 548                                 foundTrack = kTRUE;
 549                                 temp[ilabel] = 0;
 550                                 prongindex=ilabel;
 551                                 break;
 552                         }
 553                 }
 554                 if (foundTrack){
 555                         foundDaughters++;
 556                         AliDebug(4,Form("daughter %d \n",foundDaughters));
 557                         if(trackCuts[prongindex]->GetRequireTPCRefit()){
 558                                 if(track->GetStatus()&AliESDtrack::kTPCrefit) {
 559                                         bRefitStep = kTRUE;
 560                                 }
 561                                 else {
 562                                         AliDebug(3, "Refit cut not passed , missing TPC refit\n");
 563                                         delete [] temp;
 564                                         temp = 0x0;
 565                                         return kFALSE;
 566                                 }
 567                         }
 568
 569                         if (trackCuts[prongindex]->GetRequireITSRefit()) {
 570                                 if(track->GetStatus()&AliESDtrack::kITSrefit){
 571                                         bRefitStep = kTRUE;
 572                                 }
 573                                 else {
 574                                         AliDebug(3, "Refit cut not passed , missing ITS refit\n");
 575                                         delete [] temp;
 576                                         temp = 0x0;
 577                                         return kFALSE;
 578                                 }
 579                         }
 580                 }
 581                 if (foundDaughters == fProngs){
 582                         break;
 583                 }
 584         }
 585         //}
 586         delete [] temp;
 587         temp = 0x0;
 588         if (foundDaughters== fProngs)  return bRefitStep;
 589         else return kFALSE;
 590 }
 591
 592 //____________________________________________________________________________
 593
 594 Bool_t AliCFVertexingHF::RecoStep()
 595 {
 596         //
 597         //check also vertex and ITS Refit and TPC Refit
 598         //
 599
 600         Bool_t bRecoStep = kFALSE;
 601         Int_t mcLabel = GetMCLabel();
 602
 603         if (mcLabel == -1) {
 604                 AliDebug(2,"No MC particle found");
 605                 return bRecoStep;
 606         }
 607         else{
 608                 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
 609                 if (!fmcPartCandidate){
 610                         AliWarning("Could not find associated MC in AOD MC tree");
 611                         return bRecoStep;
 612                 }
 613         }
 614
 615         Int_t pdgGranma = CheckOrigin();
 616
 617         if (pdgGranma == -99999){
 618                 AliDebug(2,"This particle does not have a quark in his genealogy\n");
 619                 return bRecoStep;
 620         }
 621         if (pdgGranma == -9999){
 622                 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");
 623                 return bRecoStep;
 624         }
 625
 626         if (pdgGranma == -999){
 627                 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");
 628                 return bRecoStep;
 629         }
 630
 631         bRecoStep=kTRUE;
 632         return bRecoStep;
 633 }
 634 //____________________________________________
 635 Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const
 636 {
 637         //
 638         // getting eta of the prong
 639         //
 640
 641         if (fRecoCandidate){
 642                 Double_t etaProng = fRecoCandidate->EtaProng(iProng);
 643                 return etaProng;
 644         }
 645         return 999999;
 646 }
 647 //______________________________________________________
 648 Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const
 649 {
 650         //
 651         // getting pt of the prong
 652         //
 653
 654         if (fRecoCandidate){
 655                 Double_t ptProng = fRecoCandidate->PtProng(iProng);
 656                 return ptProng;
 657         }
 658         return 999999;
 659
 660 }
 661
 662 //____________________________________________________________________
 663
 664 Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts **trackCuts) const
 665 {
 666         //
 667         // reco Acceptance step
 668         //
 669
 670         Bool_t bRecoAccStep = kFALSE;
 671
 672         Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
 673
 674         Float_t etaProng=0., ptProng=0.;
 675
 676         for (Int_t iProng =0; iProng<fProngs; iProng++){
 677
 678                 trackCuts[iProng]->GetEtaRange(etaCutMin, etaCutMax);
 679                 trackCuts[iProng]->GetPtRange(ptCutMin, ptCutMax);
 680                 etaProng = GetEtaProng(iProng);
 681                 ptProng = GetPtProng(iProng);
 682
 683                 Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
 684                 if (!acceptanceProng) {
 685                         AliDebug(2,"At least one reconstructed prong isn't in the acceptance\n");
 686                         return bRecoAccStep;
 687                 }
 688         }
 689
 690         bRecoAccStep=kTRUE;
 691         return bRecoAccStep;
 692 }
 693 //___________________________________________________________
 694
 695 Bool_t AliCFVertexingHF::FillUnfoldingMatrix(Double_t fill[4]) const
 696 {
 697         //
 698         // filling the unfolding matrix
 699         //
 700
 701         if(fmcPartCandidate){
 702
 703                 fill[0] = GetPtCand();
 704                 fill[1] = GetYCand();
 705
 706                 fill[2] =  fmcPartCandidate->Pt();
 707                 fill[3] =  fmcPartCandidate->Y();
 708
 709                 return kTRUE;
 710         }
 711
 712         return kFALSE;
 713 }
 714 //___________________________________________________________
 715
 716 Int_t AliCFVertexingHF::CheckReflexion(Char_t isSign)
 717 {
 718         //
 719         // check for reflexion (particle/antiparticle)
 720         //
 721
 722         Int_t mcLabel = GetMCLabel();
 723
 724         if (mcLabel == -1) {
 725                 AliDebug(2,"No MC particle found");
 726                 return 0;
 727         }
 728         else{
 729                 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
 730                 if (!fmcPartCandidate){
 731                         AliWarning("Could not find associated MC in AOD MC tree");
 732                         return 0;
 733                 }
 734         }
 735
 736         if(fmcPartCandidate->GetPdgCode()>0) {
 737                 if (isSign == 1){ // I ask for antiparticle only
 738                         AliDebug(2,"candidate is particle, I ask for antiparticle only");
 739                         return 0;
 740                 }
 741                 return 1;  // particle
 742         }
 743         else if(fmcPartCandidate->GetPdgCode()<0) {
 744                 if (isSign == 0){ // I ask for particle only
 745                         AliDebug(2,"candidate is antiparticle, I ask for particle only");
 746                         return 0;
 747                 }
 748                 return 2;  // antiparticle
 749         }
 750         else return 0;  // ....shouldn't be...
 751
 752 }
 753 //___________________________________________________________
 754
 755 Bool_t AliCFVertexingHF::SetLabelArray()
 756 {
 757         //
 758         // setting the label arrays
 759         //
 760
 761         Bool_t bLabelArray = kFALSE;
 762
 763         fLabelArray = new Int_t[fProngs];
 764
 765         AliAODMCParticle *mcPartDaughter;
 766         Int_t label0 = fmcPartCandidate->GetDaughter(0);
 767         Int_t label1 = fmcPartCandidate->GetDaughter(1);
 768         AliDebug(2,Form("label0 = %d, label1 = %d",label0,label1));
 769         if (label1==0 || label0 == 0){
 770                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 771                 delete [] fLabelArray;
 772                 fLabelArray = 0x0;
 773                 return bLabelArray;
 774         }
 775
 776         if (label1 - label0 == fProngs-1){
 777                 for (Int_t iProng = 0; iProng<fProngs; iProng++){
 778                         mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
 779                         if (mcPartDaughter){
 780                                 fLabelArray[iProng] =  mcPartDaughter->GetLabel();
 781                         }
 782                         else{
 783                                 AliError("Failed casting the daughter particle, returning a NULL label array");
 784                                 delete [] fLabelArray;
 785                                 fLabelArray = 0x0;
 786                                 return bLabelArray;
 787                         }
 788                 }
 789
 790         }
 791         // resonant decay channel
 792         else if (label1 - label0 == fProngs-2 && fProngs > 2){
 793                 Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 794                 Int_t foundDaughters = 0;
 795                 for(Int_t iDau=0; iDau<fProngs-1; iDau++){
 796                         Int_t iLabelDau = labelFirstDau+iDau;
 797                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDau));
 798                         if ( ! part ) {
 799                           AliError("Wrong particle type in fmcArray");
 800                           delete [] fLabelArray;
 801                           fLabelArray = 0x0;
 802                           return bLabelArray;
 803                         }
 804                         Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 805                         if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 806                                 if (part) {
 807                                         fLabelArray[foundDaughters] = part->GetLabel();
 808                                         foundDaughters++;
 809                                 }
 810                                 else{
 811                                         AliError("Error while casting particle! returning a NULL array");
 812                                         delete [] fLabelArray;
 813                                         fLabelArray = 0x0;
 814                                         return bLabelArray;
 815                                 }
 816                         }
 817                         // added K0S case - Start
 818                         else if (pdgCode==311) {
 819                           if (part->GetNDaughters()!=1) {
 820                             delete [] fLabelArray;
 821                             fLabelArray = 0x0;
 822                             return bLabelArray;
 823                           }
 824                           Int_t labelK0Dau = part->GetDaughter(0);
 825                           AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelK0Dau));
 826                           if(!partK0S){
 827                             AliError("Error while casting particle! returning a NULL array");
 828                             delete [] fLabelArray;
 829                             fLabelArray = 0x0;
 830                             return bLabelArray;
 831                           }
 832                           Int_t nDauRes=partK0S->GetNDaughters();
 833                           if(nDauRes!=2 || partK0S->GetPdgCode()!=310) {
 834                             delete [] fLabelArray;
 835                             fLabelArray = 0x0;
 836                             return bLabelArray;
 837                           }
 838                           Int_t labelFirstDauRes = partK0S->GetDaughter(0);
 839                           AliDebug(2,Form(" Found K0S (%d)",labelK0Dau));
 840                           for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 841                             Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 842                             AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 843                             if (dauRes){
 844                               fLabelArray[foundDaughters] = dauRes->GetLabel();
 845                               foundDaughters++;
 846                             }
 847                             else {
 848                               AliError("Error while casting resonant daughter! returning a NULL array");
 849                               delete [] fLabelArray;
 850                               fLabelArray = 0x0;
 851                               return bLabelArray;
 852                             }
 853                           }
 854                         }
 855                         // added K0S case - End
 856                         else{
 857                                 Int_t nDauRes=part->GetNDaughters();
 858                                 if(nDauRes!=2) {
 859                                         delete [] fLabelArray;
 860                                         fLabelArray = 0x0;
 861                                         return bLabelArray;
 862                                 }
 863                                 Int_t labelFirstDauRes = part->GetDaughter(0);
 864                                 for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 865                                         Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 866                                         AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 867                                         if (dauRes){
 868                                                 fLabelArray[foundDaughters] = dauRes->GetLabel();
 869                                                 foundDaughters++;
 870                                         }
 871                                         else{
 872                                                 AliError("Error while casting resonant daughter! returning a NULL array");
 873                                                 delete [] fLabelArray;
 874                                                 fLabelArray = 0x0;
 875                                                 return bLabelArray;
 876                                         }
 877                                 }
 878                         }
 879                 }
 880                 if (foundDaughters != fProngs){
 881                         delete [] fLabelArray;
 882                         fLabelArray = 0x0;
 883                         return bLabelArray;
 884                 }
 885         }
 886         // wrong correspondance label <--> prongs
 887         else{
 888                 delete [] fLabelArray;
 889                 fLabelArray = 0x0;
 890                 return bLabelArray;
 891         }
 892         SetAccCut();   // setting the pt and eta acceptance cuts
 893         bLabelArray = kTRUE;
 894         return bLabelArray;
 895 }
 896
 897 //___________________________________________________________
 898
 899 void AliCFVertexingHF::SetPtAccCut(Float_t* ptAccCut)
 900 {
 901         //
 902         // setting the pt cut to be used in the Acceptance steps (MC+Reco)
 903         //
 904
 905         if (fProngs>0){
 906                 for (Int_t iP=0; iP<fProngs; iP++){
 907                         fPtAccCut[iP]=ptAccCut[iP];
 908                 }
 909         }
 910         return;
 911 }
 912
 913
 914
 915 //___________________________________________________________
 916
 917 void AliCFVertexingHF::SetEtaAccCut(Float_t* etaAccCut)
 918 {
 919         //
 920         // setting the eta cut to be used in the Acceptance steps (MC+Reco)
 921         //
 922
 923         if (fProngs>0){
 924                 for (Int_t iP=0; iP<fProngs; iP++){
 925                         fEtaAccCut[iP]=etaAccCut[iP];
 926                 }
 927         }
 928         return;
 929 }
 930 //___________________________________________________________
 931
 932 void AliCFVertexingHF::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
 933 {
 934         //
 935         // setting the pt and 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                         fPtAccCut[iP]=ptAccCut[iP];
 941                         fEtaAccCut[iP]=etaAccCut[iP];
 942                 }
 943         }
 944         return;
 945 }
 946
 947 //___________________________________________________________
 948
 949 void AliCFVertexingHF::SetAccCut()
 950 {
 951         //
 952         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 953         //
 954
 955         if (fProngs>0){
 956                 for (Int_t iP=0; iP<fProngs; iP++){
 957                         fPtAccCut[iP]=0.1;
 958                         fEtaAccCut[iP]=0.9;
 959                 }
 960         }
 961         return;
 962 }