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=0, ptCutMin=0, etaCutMax=0, ptCutMax=0;
 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         if (label1<=0 || label0 <= 0){
 793                 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
 794                 delete [] fLabelArray;
 795                 fLabelArray = 0x0;
 796                 return bLabelArray;
 797         }
 798         AliAODMCParticle* tmp0 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0));
 799         AliAODMCParticle* tmp1 = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label1));
 800         AliDebug(2, Form("label0 = %d (pdg = %d), label1 = %d (pdg = %d)", label0, tmp0->GetPdgCode(), label1, tmp1->GetPdgCode()));
 801
 802         if (label1 - label0 == fProngs-1){
 803                 for (Int_t iProng = 0; iProng<fProngs; iProng++){
 804                         mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
 805                         if (mcPartDaughter){
 806                                 fLabelArray[iProng] =  mcPartDaughter->GetLabel();
 807                         }
 808                         else{
 809                                 AliError("Failed casting the daughter particle, returning a NULL label array");
 810                                 delete [] fLabelArray;
 811                                 fLabelArray = 0x0;
 812                                 return bLabelArray;
 813                         }
 814                 }
 815
 816         }
 817         // resonant decay channel
 818         else if (label1 - label0 == fProngs-2 && fProngs > 2){
 819           AliDebug(3, "In the resonance decay channel");
 820                 Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 821                 Int_t foundDaughters = 0;
 822                 for(Int_t iDau=0; iDau<fProngs-1; iDau++){
 823                         Int_t iLabelDau = labelFirstDau+iDau;
 824                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDau));
 825                         if ( ! part ) {
 826                           AliError("Wrong particle type in fmcArray");
 827                           delete [] fLabelArray;
 828                           fLabelArray = 0x0;
 829                           return bLabelArray;
 830                         }
 831                         Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 832                         AliDebug(3, Form("Prong %d had pdg = %d", iDau, pdgCode));
 833                         if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 834                                 if (part) {
 835                                         fLabelArray[foundDaughters] = part->GetLabel();
 836                                         AliDebug(3, Form("part found at %d has label = %d", iLabelDau, part->GetLabel()));
 837                                         AliDebug(3, Form("fLabelArray[%d] = %d", foundDaughters, fLabelArray[foundDaughters]));
 838                                         foundDaughters++;
 839                                 }
 840                                 else{
 841                                         AliError("Error while casting particle! returning a NULL array");
 842                                         delete [] fLabelArray;
 843                                         fLabelArray = 0x0;
 844                                         return bLabelArray;
 845                                 }
 846                         }
 847                         // added K0S case - Start
 848                         else if (pdgCode==311) {
 849                           AliDebug(3, Form("K0S case, foundDaughters = %d", foundDaughters));
 850                           if (part->GetNDaughters()!=1) {
 851                             delete [] fLabelArray;
 852                             fLabelArray = 0x0;
 853                             return bLabelArray;
 854                           }
 855                           Int_t labelK0Dau = part->GetDaughter(0);
 856                           AliAODMCParticle* partK0S = dynamic_cast<AliAODMCParticle*>(fmcArray->At(labelK0Dau));
 857                           if(!partK0S){
 858                             AliError("Error while casting particle! returning a NULL array");
 859                             delete [] fLabelArray;
 860                             fLabelArray = 0x0;
 861                             return bLabelArray;
 862                           }
 863                           Int_t nDauRes = partK0S->GetNDaughters();
 864                           AliDebug(3, Form("nDauRes = %d", nDauRes));
 865                           if (nDauRes!=2 || partK0S->GetPdgCode() != 310) {
 866                             AliDebug(2, "No K0S on no 2-body decay");
 867                             delete [] fLabelArray;
 868                             fLabelArray = 0x0;
 869                             return bLabelArray;
 870                           }
 871                           Int_t labelFirstDauRes = partK0S->GetDaughter(0);
 872                           AliDebug(2, Form("Found K0S (%d)", labelK0Dau));
 873                           for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 874                             Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 875                             AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 876                             AliDebug(3, Form("daughter = %d, pointer = %p, with label = %d", iLabelDauRes, dauRes, dauRes->GetLabel()));
 877                             if (dauRes){
 878                               AliDebug(3, Form("PDG code = %d", dauRes->GetPdgCode()));
 879                               if (TMath::Abs(dauRes->GetPdgCode())!=211) {
 880                                 AliDebug(2,"K0S doesn't decay in 2 charged pions!");
 881                                 delete [] fLabelArray;
 882                                 fLabelArray = 0x0;
 883                                 return bLabelArray;
 884                               }
 885                               else {
 886                                 fLabelArray[foundDaughters] = iLabelDauRes;  // N.B.: do not use dauRes->GetLabel()!!!! it is wrong!!!
 887                                 AliDebug(3, Form("Setting fLabelArray[%d] = %d (before it was %d)", foundDaughters, iLabelDauRes, dauRes->GetLabel()));
 888                                 foundDaughters++;
 889                               }
 890                             }
 891                             else {
 892                               AliError("Error while casting resonant daughter! returning a NULL array");
 893                               delete [] fLabelArray;
 894                               fLabelArray = 0x0;
 895                               return bLabelArray;
 896                             }
 897                           }
 898                         }
 899                         // added K0S case - End
 900                         else{
 901                           Int_t nDauRes=part->GetNDaughters();
 902                           AliDebug(3, Form("nDauRes = %d", nDauRes));
 903                           if(nDauRes!=2) {
 904                             AliDebug(3, Form("nDauRes = %d, different from 2", nDauRes));
 905                             Int_t labelFirstDauResTest = part->GetDaughter(0);
 906                             for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 907                               Int_t iLabelDauResTest = labelFirstDauResTest+iDauRes;
 908                               AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauResTest));
 909                               if (dauRes){
 910                                 AliDebug(3, Form("pdg of daugh %d = %d", iDauRes, dauRes->GetPdgCode()));
 911                               }
 912                             }
 913                             delete [] fLabelArray;
 914                             fLabelArray = 0x0;
 915                             return bLabelArray;
 916                           }
 917                           Int_t labelFirstDauRes = part->GetDaughter(0);
 918                           for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
 919                             Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
 920                             AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
 921                             if (dauRes){
 922                               fLabelArray[foundDaughters] = dauRes->GetLabel();
 923                               foundDaughters++;
 924                             }
 925                             else{
 926                               AliError("Error while casting resonant daughter! returning a NULL array");
 927                               delete [] fLabelArray;
 928                               fLabelArray = 0x0;
 929                               return bLabelArray;
 930                             }
 931                           }
 932                         }
 933                 }
 934                 if (foundDaughters != fProngs){
 935                   AliDebug(3, Form("foundDaughters (%d) != fProngs (%d)", foundDaughters, fProngs));
 936                   delete [] fLabelArray;
 937                   fLabelArray = 0x0;
 938                   return bLabelArray;
 939                 }
 940         }
 941         // wrong correspondance label <--> prongs
 942         else{
 943           delete [] fLabelArray;
 944           fLabelArray = 0x0;
 945           return bLabelArray;
 946         }
 947         AliDebug(3, "Setting AccCuts");
 948         SetAccCut();   // setting the pt and eta acceptance cuts
 949         bLabelArray = kTRUE;
 950         return bLabelArray;
 951 }
 952
 953 //___________________________________________________________
 954
 955 void AliCFVertexingHF::SetPtAccCut(Float_t* ptAccCut)
 956 {
 957         //
 958         // setting the pt cut to be used in the Acceptance steps (MC+Reco)
 959         //
 960
 961         if (fProngs>0){
 962                 for (Int_t iP=0; iP<fProngs; iP++){
 963                         fPtAccCut[iP]=ptAccCut[iP];
 964                 }
 965         }
 966         return;
 967 }
 968
 969
 970
 971 //___________________________________________________________
 972
 973 void AliCFVertexingHF::SetEtaAccCut(Float_t* etaAccCut)
 974 {
 975         //
 976         // setting the eta cut to be used in the Acceptance steps (MC+Reco)
 977         //
 978
 979         if (fProngs>0){
 980                 for (Int_t iP=0; iP<fProngs; iP++){
 981                         fEtaAccCut[iP]=etaAccCut[iP];
 982                 }
 983         }
 984         return;
 985 }
 986 //___________________________________________________________
 987
 988 void AliCFVertexingHF::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
 989 {
 990         //
 991         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
 992         //
 993
 994         if (fProngs>0){
 995                 for (Int_t iP=0; iP<fProngs; iP++){
 996                         fPtAccCut[iP]=ptAccCut[iP];
 997                         fEtaAccCut[iP]=etaAccCut[iP];
 998                 }
 999         }
1000         return;
1001 }
1002
1003 //___________________________________________________________
1004
1005 void AliCFVertexingHF::SetAccCut()
1006 {
1007         //
1008         // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
1009         //
1010
1011         if (fProngs>0){
1012                 for (Int_t iP=0; iP<fProngs; iP++){
1013                         fPtAccCut[iP]=0.1;
1014                         fEtaAccCut[iP]=0.9;
1015                 }
1016         }
1017         return;
1018 }