PWG3/vertexingHF/AliCFVertexingHF3Prong.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2011, 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 //                                                               //
  20 // Class to compute variables for correction framework           //
  21 // for 3-body decays of D mesons (D+, Ds, Lc)                    //
  22 // in bins of cut variables                                      //
  23 // Origin:       Francesco Prino (prino@to.infn.it)              //
  24 //               Renu Bala       (bala@to.infn.it)               //
  25 //               Davide Caffarri (cafarri@pd.infn.it)            //
  26 ///////////////////////////////////////////////////////////////////
  27
  28 #include "AliAODMCParticle.h"
  29 #include "AliAODEvent.h"
  30 #include "TClonesArray.h"
  31 #include "AliCFVertexingHF.h"
  32 #include "AliESDtrack.h"
  33 #include "TDatabasePDG.h"
  34
  35 #include "AliCFVertexingHF3Prong.h"
  36 #include "AliCFContainer.h"
  37 #include "AliCFTaskVertexingHF.h"
  38
  39 ClassImp(AliCFVertexingHF3Prong)
  40
  41
  42 //_________________________________________
  43 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(Int_t decay):
  44 AliCFVertexingHF(),
  45   fDecay(decay)
  46  {
  47   //
  48   SetNProngs(3);
  49
  50   fPtAccCut=new Float_t[fProngs];
  51   fEtaAccCut=new Float_t[fProngs];
  52   for(Int_t iP=0; iP<fProngs; iP++){
  53           fPtAccCut[iP]=0.1;
  54           fEtaAccCut[iP]=0.9;
  55   }
  56
  57 }
  58 //_________________________________________
  59 AliCFVertexingHF3Prong::AliCFVertexingHF3Prong(TClonesArray *mcArray, UShort_t originDselection, Int_t decay):
  60   AliCFVertexingHF(mcArray, originDselection),
  61   fDecay(decay)
  62
  63 {
  64   //
  65   SetNProngs(3);
  66   fPtAccCut=new Float_t[fProngs];
  67   fEtaAccCut=new Float_t[fProngs];
  68   for(Int_t iP=0; iP<fProngs; iP++){
  69           fPtAccCut[iP]=0.1;
  70           fEtaAccCut[iP]=0.9;
  71   }
  72 }
  73
  74
  75 //_____________________________________
  76 AliCFVertexingHF3Prong& AliCFVertexingHF3Prong::operator=(const AliCFVertexingHF3Prong& c){
  77   //
  78   if  (this != &c) {
  79
  80     AliCFVertexingHF::operator=(c);
  81
  82   }
  83     return *this;
  84 }
  85
  86 //__________________________________________
  87 Bool_t AliCFVertexingHF3Prong::SetRecoCandidateParam(AliAODRecoDecayHF *recoCand){
  88   // Checks if candidate is signal and D meson is present in MC array
  89
  90   Bool_t bSignAssoc = kFALSE;
  91   fRecoCandidate = recoCand;
  92
  93   if (!fRecoCandidate) {
  94     AliError("fRecoCandidate not found, problem in assignement\n");
  95     return bSignAssoc;
  96   }
  97
  98   Int_t pdgCand = -1;
  99   Int_t pdgDaughter[3]={-1,-1,-1};
 100   if(fDecay==kDplustoKpipi){
 101     pdgCand=411;
 102     pdgDaughter[0]=321;
 103     pdgDaughter[1]=211;
 104     pdgDaughter[2]=211;
 105   }else if(fDecay==kDstoKKpi){
 106     pdgCand=431;
 107     pdgDaughter[0]=321;
 108     pdgDaughter[1]=321;
 109     pdgDaughter[2]=211;
 110   }else if(fDecay==kLctopKpi){
 111           pdgCand=4122;
 112           pdgDaughter[0]=2212;
 113           pdgDaughter[1]=321;
 114           pdgDaughter[2]=211;
 115   }else{
 116     AliError("WRONG DECAY SETTING");
 117     return bSignAssoc;
 118   }
 119
 120   Int_t mcLabel = fRecoCandidate->MatchToMC(pdgCand,fmcArray,3,pdgDaughter);
 121   if (mcLabel == -1) return bSignAssoc;
 122
 123   if (fRecoCandidate->NumberOfFakeDaughters()>0){
 124           fFake = 0;    // fake candidate
 125           if (fFakeSelection==1) return bSignAssoc;
 126   }
 127   if (fRecoCandidate->NumberOfFakeDaughters()==0){
 128           fFake = 2;    // non-fake candidate
 129           if (fFakeSelection==2) return bSignAssoc;
 130   }
 131
 132   SetMCLabel(mcLabel);
 133   fmcPartCandidate = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fmcLabel));
 134
 135   if (!fmcPartCandidate){
 136     AliDebug(3,"No part candidate");
 137     return bSignAssoc;
 138   }
 139
 140   bSignAssoc = kTRUE;
 141   return bSignAssoc;
 142 }
 143
 144 //______________________________________________
 145 Bool_t AliCFVertexingHF3Prong::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) {
 146         //
 147         // collecting all the necessary info from MC particle and fill vectorMC: 12 variables
 148         // pt_D
 149         // y_D
 150         // phi_D
 151         // ctau
 152         // cos point
 153         // pt_1
 154         // pt_2
 155         // pt_3
 156         // d0_1
 157         // d0_2
 158         // d0_3
 159         // zPrimVert
 160         // centrality
 161
 162         Bool_t bGenValues = kFALSE;
 163
 164         Int_t pdgCand = -1;
 165         if(fDecay==kDplustoKpipi){
 166                 pdgCand=411;
 167         }else if(fDecay==kDstoKKpi){
 168                 pdgCand=431;
 169         }else if(fDecay==kLctopKpi){
 170                 pdgCand=4122;
 171         }else{
 172                 AliError("WRONG DECAY SETTING");
 173                 return bGenValues;
 174         }
 175
 176         Double_t vertD[3] = {0,0,0};   // D origin
 177         fmcPartCandidate->XvYvZv(vertD);  // cm
 178
 179         Int_t nprongs = 3;
 180         Int_t daughter[3];
 181         Short_t charge = fmcPartCandidate->Charge();
 182
 183         // order the daughters as LS,OS,LS, e.g. D+ -> pi+ K- pi+
 184         // the 2 LS are ordered so that in pos. 0 there is the one with lower label value
 185         Int_t index=0;
 186         Int_t nDauLS=0;
 187         Int_t nDauOS=0;
 188
 189
 190         Int_t nDau=fmcPartCandidate->GetNDaughters();
 191         Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 192         if(nDau==3){
 193                 for(Int_t iDau=0; iDau<3; iDau++){
 194                         Int_t ind = labelFirstDau+iDau;
 195                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 196                         if(!part){
 197                                 AliError("Daughter particle not found in MC array");
 198                                 return bGenValues;
 199                         }
 200                         Short_t signDau=part->Charge();
 201                         if(signDau==charge){
 202                                 nDauLS++;
 203                                 daughter[index] = ind;
 204                                 index=2;
 205                         }else{
 206                                 daughter[1] = ind;
 207                                 nDauOS++;
 208                         }
 209                 }
 210         }else if(nDau==2){
 211                 for(Int_t iDau=0; iDau<2; iDau++){
 212                         Int_t ind = labelFirstDau+iDau;
 213                         AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 214                         if(!part){
 215                                 AliError("Daughter particle not found in MC array");
 216                                 return bGenValues;
 217                         }
 218                         Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 219                         if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 220                                 Short_t signDau=part->Charge();
 221                                 if(signDau==charge){
 222                                         nDauLS++;
 223                                         daughter[index] = ind;
 224                                         index=2;
 225                                 }else{
 226                                         daughter[1] = ind;
 227                                         nDauOS++;
 228                                 }
 229                         }else{
 230                                 Int_t nDauRes=part->GetNDaughters();
 231                                 if(nDauRes!=2){
 232                                         AliError("Wrong resonant decay");
 233                                         return bGenValues;
 234                                 }
 235                                 Int_t labelFirstDauRes = part->GetDaughter(0);
 236                                 for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 237                                         Int_t indDR = labelFirstDauRes+iDauRes;
 238                                         AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
 239                                         if(!partDR){
 240                                                 AliError("Daughter particle not found in MC array");
 241                                                 return bGenValues;
 242                                         }
 243                                         Short_t signDau=partDR->Charge();
 244                                         if(signDau==charge){
 245                                                 nDauLS++;
 246                                                 daughter[index] = ind;
 247                                                 index=2;
 248                                         }else{
 249                                                 daughter[1] = ind;
 250                                                 nDauOS++;
 251                                         }
 252                                 }
 253                         }
 254                 }
 255         }else{
 256                 AliError(Form("Wrong number of daughters %d",nDau));
 257                 return bGenValues;
 258         }
 259
 260         if(nDauLS!=2 || nDauOS!=1){
 261                 AliError(Form("Wrong decay channel: LS and OS daughters not OK: %d %d",nDauLS,nDauOS));
 262                 return bGenValues;
 263         }
 264         if(daughter[0]>daughter[2]){
 265                 Int_t tmp=daughter[0];
 266                 daughter[0]=daughter[2];
 267                 daughter[2]=tmp;
 268         }
 269
 270         // getting the momentum from the daughters and decay vertex
 271         Double_t px[3],py[3],pz[3],pt[3];
 272         Double_t vertDec[3] = {0,0,0};   // decay vertex
 273         for(Int_t iDau=0; iDau<3; iDau++){
 274                 AliAODMCParticle* part=dynamic_cast<AliAODMCParticle*>(fmcArray->At(daughter[iDau]));
 275                 if(!part){
 276                         AliError("Daughter particle not found in MC array");
 277                         return bGenValues;
 278                 }
 279                 px[iDau]=part->Px();
 280                 py[iDau]=part->Py();
 281                 pz[iDau]=part->Pz();
 282                 pt[iDau]=part->Pt();
 283                 if(iDau==0) part->XvYvZv(vertDec);
 284         }
 285
 286         Double_t d0[3] = {0.,0.,0.}; // dummy values!!!!
 287
 288         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vertD,vertDec,nprongs,charge,px,py,pz,d0);
 289         Double_t cT = decay->Ct(pdgCand);
 290
 291         switch (fConfiguration){
 292         case AliCFTaskVertexingHF::kSnail:
 293                 vectorMC[0] = fmcPartCandidate->Pt();
 294                 vectorMC[1] = fmcPartCandidate->Y() ;
 295                 vectorMC[2] = fmcPartCandidate->Phi();
 296                 vectorMC[3] = cT*1.E4 ;  // in micron
 297                 vectorMC[4] = 1.01;    // cos pointing angle, dummy value, meaningless in MC
 298                 vectorMC[5] = pt[0];
 299                 vectorMC[6] = pt[1];
 300                 vectorMC[7] = pt[2];
 301                 vectorMC[8] = fzMCVertex;    // z of reconstructed of primary vertex
 302                 vectorMC[9] = fCentValue; // reconstructed centrality value
 303                 vectorMC[10] = 1.;           // fake: always filling with 1 at MC level
 304                 vectorMC[11] = 1.01; // dummy value for cosPointingXY  multiplicity
 305                 vectorMC[12] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
 306                 vectorMC[13] = fMultiplicity; // reconstructed multiplicity
 307
 308                 if (fDecay==kLctopKpi){
 309                         vectorMC[11] = 0.; //dist12
 310                         vectorMC[12] = 0.; //dist23
 311                         vectorMC[13] = 0.; //sigmaVtx
 312                         vectorMC[14] = 0.; //sumd02
 313                         vectorMC[15] = 1.01; // dummy value for cosPointingXY  multiplicity
 314                         vectorMC[16] = 0.; // dummy value for NormalizedDecayLengthXY multiplicity
 315                         vectorMC[17] = fMultiplicity; // reconstructed multiplicity
 316                 }
 317                 break;
 318
 319         case AliCFTaskVertexingHF::kCheetah:
 320                 vectorMC[0] = fmcPartCandidate->Pt();
 321                 vectorMC[1] = fmcPartCandidate->Y() ;
 322                 vectorMC[2] = cT*1.E4; // in micron
 323                 vectorMC[3] = fmcPartCandidate->Phi();
 324                 vectorMC[4] = fzMCVertex;
 325                 vectorMC[5] = fCentValue;   // dummy value for dca, meaningless in MC
 326                 vectorMC[6] = 1. ;  // fake: always filling with 1 at MC level
 327                 vectorMC[7] = fMultiplicity;   // dummy value for d0pi, meaningless in MC, in micron
 328                 break;
 329         }
 330
 331         bGenValues = kTRUE;
 332         return bGenValues;
 333 }
 334
 335
 336 //____________________________________________
 337 Bool_t AliCFVertexingHF3Prong::GetRecoValuesFromCandidate(Double_t *vectorReco) const
 338 {
 339         // Fill vector (see above) with reconstructed quantities
 340         Bool_t bFillRecoValues=kFALSE;
 341
 342         Int_t pdgCand = -1;
 343         if(fDecay==kDplustoKpipi){
 344                 pdgCand=411;
 345         }else if(fDecay==kDstoKKpi){
 346                 pdgCand=431;
 347         }else if(fDecay==kLctopKpi){
 348                 pdgCand=4122;
 349                 // AliError("LambdaC not yet implemented");
 350                 // return bFillRecoValues;
 351         }else{
 352                 AliError("WRONG DECAY SETTING");
 353                 return bFillRecoValues;
 354         }
 355
 356         AliAODRecoDecayHF3Prong *decay3 = (AliAODRecoDecayHF3Prong*)fRecoCandidate;
 357         Short_t charge=decay3->Charge();
 358         Double_t rapidity=decay3->Y(pdgCand);
 359         Double_t cT=decay3->Ct(pdgCand);
 360         Double_t pt = decay3->Pt();
 361         Double_t cosPointingAngle = decay3->CosPointingAngle();
 362         Double_t phi = decay3->Phi();
 363         Double_t dist12= decay3->GetDist12toPrim();
 364         Double_t dist23 = decay3->GetDist23toPrim();
 365         Double_t sigmVert = decay3->GetSigmaVert();
 366         Double_t cosPointingAngleXY = decay3->CosPointingAngleXY();
 367         Double_t normDecayLengthXY = decay3->NormalizedDecayLengthXY();
 368
 369         Int_t daughtSorted[3];
 370         Int_t tmpIndex=0;
 371         Int_t nDauLS=0;
 372         Int_t nDauOS=0;
 373         for(Int_t iDau=0; iDau<3; iDau++){
 374                 AliAODTrack *trk = (AliAODTrack*)decay3->GetDaughter(iDau);
 375                 Int_t label = TMath::Abs(trk->GetLabel());
 376                 Short_t chargedau=trk->Charge();
 377                 if(chargedau==charge){
 378                         daughtSorted[tmpIndex]=label;
 379                         tmpIndex=2;
 380                         nDauLS++;
 381                 }else{
 382                         daughtSorted[1]=label;
 383                         nDauOS++;
 384                 }
 385         }
 386
 387         if(nDauLS!=2 || nDauOS!=1){
 388                 AliError("Wrong decay channel: number of OS and LS tracks not OK");
 389                 return bFillRecoValues;
 390         }
 391
 392         if(daughtSorted[0]>daughtSorted[2]){
 393                 Int_t tmp=daughtSorted[0];
 394                 daughtSorted[0]=daughtSorted[2];
 395                 daughtSorted[2]=tmp;
 396         }
 397
 398         Double_t d0prong0 = decay3->Getd0Prong(daughtSorted[0]);
 399         Double_t d0prong1 = decay3->Getd0Prong(daughtSorted[1]);
 400         Double_t d0prong2 = decay3->Getd0Prong(daughtSorted[2]);
 401
 402         switch (fConfiguration){
 403         case AliCFTaskVertexingHF::kSnail:
 404                 vectorReco[0] = pt;
 405                 vectorReco[1] = rapidity;
 406                 vectorReco[2] = phi;
 407                 vectorReco[3] = cT*1.E4;  // in micron
 408                 vectorReco[4] = cosPointingAngle;  // in micron
 409                 vectorReco[5] = decay3->PtProng(daughtSorted[0]);
 410                 vectorReco[6] = decay3->PtProng(daughtSorted[1]);
 411                 vectorReco[7] = decay3->PtProng(daughtSorted[2]);
 412                 vectorReco[8] = fzPrimVertex;    // z of reconstructed of primary vertex
 413                 vectorReco[9] = fCentValue; //reconstructed centrality value
 414                 vectorReco[10] = fFake;      // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2)
 415                 vectorReco[11] = cosPointingAngleXY;
 416                 vectorReco[12] = normDecayLengthXY; // in cm
 417                 vectorReco[13] = fMultiplicity; // reconstructed multiplicity
 418
 419                 if(fDecay==kLctopKpi){
 420                         Double_t sumd02 =(d0prong0*d0prong0 + d0prong1*d0prong1 + d0prong2*d0prong2);
 421                         vectorReco[11] = dist12*1.E4;
 422                         vectorReco[12] = dist23*1.E4;
 423                         vectorReco[13] = sigmVert*1.E4;
 424                         vectorReco[14] = sumd02*1.E8;
 425                         vectorReco[15] = cosPointingAngleXY;
 426                         vectorReco[16] = normDecayLengthXY; // in cm
 427                         vectorReco[17] = fMultiplicity; // reconstructed multiplicity
 428                 }
 429                 break;
 430         case AliCFTaskVertexingHF::kCheetah:
 431                 vectorReco[0] = pt;
 432                 vectorReco[1] = rapidity ;
 433                 vectorReco[2] = cT*1.E4; // in micron
 434                 vectorReco[3] = phi;
 435                 vectorReco[4] = fzPrimVertex;
 436                 vectorReco[5] = fCentValue;
 437                 vectorReco[6] = fFake ;
 438                 vectorReco[7] = fMultiplicity;
 439                 break;
 440         }
 441
 442         bFillRecoValues = kTRUE;
 443         return bFillRecoValues;
 444 }
 445
 446
 447 //_____________________________________________________________
 448 Bool_t AliCFVertexingHF3Prong::CheckMCChannelDecay() const
 449 {
 450   // Check the pdg codes of the daughters
 451   Bool_t checkCD = kFALSE;
 452
 453   Int_t pdgCand = -1;
 454   Int_t pdgDaughter[3]={-1,-1,-1};
 455   if(fDecay==kDplustoKpipi){
 456     pdgCand=411;
 457     pdgDaughter[0]=321;
 458     pdgDaughter[1]=211;
 459     pdgDaughter[2]=211;
 460   }else if(fDecay==kDstoKKpi){
 461     pdgCand=431;
 462     pdgDaughter[0]=321;
 463     pdgDaughter[1]=321;
 464     pdgDaughter[2]=211;
 465   }else if(fDecay==kLctopKpi){
 466           pdgCand=4122;
 467           pdgDaughter[0]=2212;
 468           pdgDaughter[1]=321;
 469           pdgDaughter[2]=211;
 470
 471   //  AliError("LambdaC not yet implemented");
 472   //  return checkCD;
 473   }else{
 474     AliError("WRONG DECAY SETTING");
 475     return checkCD;
 476   }
 477
 478
 479   Int_t daughter[3];
 480
 481   Int_t nDau=fmcPartCandidate->GetNDaughters();
 482   Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
 483   if(nDau==3){
 484     for(Int_t iDau=0; iDau<3; iDau++){
 485       Int_t ind = labelFirstDau+iDau;
 486       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 487       if(!part){
 488         AliError("Daughter particle not found in MC array");
 489         return checkCD;
 490       }
 491       daughter[iDau]=TMath::Abs(part->GetPdgCode());
 492     }
 493   }else if(nDau==2){
 494     Int_t nDauFound=0;
 495     for(Int_t iDau=0; iDau<2; iDau++){
 496       Int_t ind = labelFirstDau+iDau;
 497       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(ind));
 498       if(!part){
 499         AliError("Daughter particle not found in MC array");
 500         return checkCD;
 501       }
 502       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
 503       if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
 504         if(nDauFound>=3) return checkCD;
 505         daughter[nDauFound]=pdgCode;
 506         nDauFound++;
 507       }else{
 508         Int_t nDauRes=part->GetNDaughters();
 509         if(nDauRes!=2) return checkCD;
 510         Int_t labelFirstDauRes = part->GetDaughter(0);
 511         for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
 512           Int_t indDR = labelFirstDauRes+iDauRes;
 513           AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(fmcArray->At(indDR));
 514           if(!partDR){
 515             AliError("Daughter particle not found in MC array");
 516             return checkCD;
 517           }
 518           Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
 519           if(nDauFound>=3) return checkCD;
 520           daughter[nDauFound]=pdgCodeDR;
 521           nDauFound++;
 522         }
 523       }
 524     }
 525   }else{
 526     return checkCD;
 527   }
 528   for(Int_t iDau1=0; iDau1<3; iDau1++){
 529     for(Int_t iDau2=iDau1; iDau2<3; iDau2++){
 530       if(daughter[iDau1]<daughter[iDau2]){
 531         Int_t tmp=daughter[iDau1];
 532         daughter[iDau1]=daughter[iDau2];
 533         daughter[iDau2]=tmp;
 534       }
 535     }
 536   }
 537   for(Int_t iDau=0; iDau<3; iDau++){
 538     if(daughter[iDau]!=pdgDaughter[iDau]){
 539       AliDebug(2, "Wrong decay channel from MC, skipping!!");
 540       return checkCD;
 541     }
 542   }
 543
 544   checkCD = kTRUE;
 545   return checkCD;
 546
 547 }