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