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