PWGHF/vertexingHF/AliAODRecoCascadeHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 for AOD reconstructed heavy-flavour cascades
  21 //
  22 // Author: X-M. Zhang, zhangxm@iopp.ccnu.edu.cn
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TVector3.h>
  26 #include <TDatabasePDG.h>
  27 #include <TClonesArray.h>
  28 #include "AliAODMCParticle.h"
  29 #include "AliAODRecoDecay.h"
  30 #include "AliAODVertex.h"
  31 #include "AliAODRecoDecayHF2Prong.h"
  32 #include "AliAODRecoCascadeHF.h"
  33
  34 ClassImp(AliAODRecoCascadeHF)
  35 //-----------------------------------------------------------------------------
  36
  37 AliAODRecoCascadeHF::AliAODRecoCascadeHF() :
  38   AliAODRecoDecayHF2Prong()
  39 {
  40   //
  41   // Default Constructor
  42   //
  43 }
  44 //-----------------------------------------------------------------------------
  45 AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
  46                                          Double_t *px, Double_t *py, Double_t *pz,
  47                                          Double_t *d0, Double_t *d0err, Double_t dca) :
  48   AliAODRecoDecayHF2Prong(vtx2, px, py, pz, d0, d0err, dca)
  49 {
  50   //
  51   //  Constructor with AliAODVertex for decay vertex
  52   //
  53   SetCharge(charge);
  54 }
  55 //-----------------------------------------------------------------------------
  56 AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
  57                                          Double_t *d0, Double_t *d0err, Double_t dca) :
  58   AliAODRecoDecayHF2Prong(vtx2, d0, d0err, dca)
  59 {
  60   //
  61   //  Constructor with decay vertex and without prongs momenta
  62   //
  63   SetCharge(charge);
  64 }
  65 //-----------------------------------------------------------------------------
  66 AliAODRecoCascadeHF::AliAODRecoCascadeHF(const AliAODRecoCascadeHF &source) :
  67   AliAODRecoDecayHF2Prong(source)
  68 {
  69   //
  70   // Copy constructor
  71   //
  72 }
  73 //-----------------------------------------------------------------------------
  74 AliAODRecoCascadeHF &AliAODRecoCascadeHF::operator=(const AliAODRecoCascadeHF &source)
  75 {
  76   //
  77   // assignment operator
  78   //
  79   if(&source == this) return *this;
  80
  81   AliAODRecoDecayHF2Prong::operator=(source);
  82
  83   return *this;
  84 }
  85 //-----------------------------------------------------------------------------
  86 AliAODRecoCascadeHF::~AliAODRecoCascadeHF()
  87 {
  88   //
  89   // Default Destructor
  90   //
  91 }
  92 //-----------------------------------------------------------------------------
  93 Double_t AliAODRecoCascadeHF::InvMassDstarKpipi() const
  94 {
  95   //
  96   // 3 prong invariant mass of the D0 daughters and the soft pion
  97   //
  98   Double_t e[3];
  99   if (Charge()>0){
 100     e[0]=Get2Prong()->EProng(0,211);
 101     e[1]=Get2Prong()->EProng(1,321);
 102   }else{
 103     e[0]=Get2Prong()->EProng(0,321);
 104     e[1]=Get2Prong()->EProng(1,211);
 105   }
 106   e[2]=EProng(0,211);
 107
 108   Double_t esum = e[0]+e[1]+e[2];
 109   Double_t minv = TMath::Sqrt(esum*esum-P()*P());
 110
 111   return minv;
 112 }
 113 //----------------------------------------------------------------------------
 114 Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
 115                                      Int_t *pdgDg,Int_t *pdgDg2prong,
 116                                      TClonesArray *mcArray, Bool_t isV0) const
 117 {
 118   //
 119   // Check if this candidate is matched to a MC signal
 120   // If no, return -1
 121   // If yes, return label (>=0) of the AliAODMCParticle
 122   //
 123
 124   Int_t ndg=GetNDaughters();
 125   if(ndg==0) {
 126     AliError("No daughters available");
 127     return -1;
 128   }
 129
 130   if ( isV0 &&
 131        ( (pdgDg[1]==2212 && pdgDg[0]==310) ||
 132          (pdgDg[1]==211 && pdgDg[0]==3122) ) ) {
 133     AliWarning("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping!");
 134     return -1;
 135   }
 136
 137   Int_t lab2Prong = -1;
 138
 139   if (!isV0) {
 140     AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
 141     lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
 142   } else {
 143     AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
 144     lab2Prong = theV0->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong); // the V0
 145   }
 146
 147   if(lab2Prong<0) return -1;
 148
 149   Int_t dgLabels[10]={0,0,0,0,0,0,0,0,0,0};
 150
 151   if (!isV0) {
 152     // loop on daughters and write labels
 153     for(Int_t i=0; i<ndg; i++) {
 154       AliVTrack *trk = dynamic_cast<AliVTrack*>(GetDaughter(i));
 155       if(!trk) continue;
 156       Int_t lab = trk->GetLabel();
 157       if(lab==-1) { // this daughter is the 2prong
 158         lab=lab2Prong;
 159       } else if(lab<-1) continue;
 160       dgLabels[i] = lab;
 161     }
 162   } else {
 163     AliVTrack *trk = dynamic_cast<AliVTrack*>(GetBachelor()); // the bachelor
 164     if (!trk) return -1;
 165     dgLabels[0] = trk->GetLabel();//TMath::Abs(trk->GetLabel());
 166     dgLabels[1] = lab2Prong;
 167   }
 168
 169   Int_t finalLabel = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
 170
 171   if (finalLabel>=0){
 172     // debug printouts for Lc->V0 bachelor case
 173
 174     if ( isV0 && (dgLabels[0]!=-1 && dgLabels[1]!=-1) ) {
 175       AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
 176       Bool_t onTheFly = theV0->GetOnFlyStatus();
 177       if (pdgDg[0]==2212 && pdgDg[1]==310) {
 178         AliAODMCParticle*k0s = dynamic_cast<AliAODMCParticle*>(mcArray->At(lab2Prong));
 179         if(k0s){
 180           Int_t labK0 = k0s->GetMother();
 181           AliAODMCParticle*k0bar = dynamic_cast<AliAODMCParticle*>(mcArray->At(labK0));
 182           if(k0bar){
 183             AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d) -> LabelK0S=%d (%d -> %d %d)",onTheFly,labK0,k0bar->GetPdgCode(),lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
 184             AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
 185                             finalLabel,pdgabs,
 186                             dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
 187           }
 188         }
 189       } else if (pdgDg[0]==211 && pdgDg[1]==3122) {
 190         AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d -> %d %d)",onTheFly,lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
 191         AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
 192                         finalLabel,pdgabs,
 193                       dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
 194       }
 195
 196     }
 197   }
 198
 199   return finalLabel;
 200
 201 }
 202 //-----------------------------------------------------------------------------
 203 Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
 204                                         const Double_t *cutsD0,
 205                                         Bool_t testD0) const
 206 {
 207   //
 208   // cutsDstar[0] = inv. mass half width of D* [GeV]
 209   // cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
 210   // cutsDstar[2] = PtMin of pi_s [GeV/c]
 211   // cutsDstar[3] = PtMax of pi_s [GeV/c]
 212   // cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
 213   //
 214   // cutsD0[0] = inv. mass half width [GeV]
 215   // cutsD0[1] = dca [cm]
 216   // cutsD0[2] = cosThetaStar
 217   // cutsD0[3] = pTK [GeV/c]
 218   // cutsD0[4] = pTPi [GeV/c]
 219   // cutsD0[5] = d0K [cm]   upper limit!
 220   // cutsD0[6] = d0Pi [cm]  upper limit!
 221   // cutsD0[7] = d0d0 [cm^2]
 222   // cutsD0[8] = cosThetaPoint
 223
 224
 225   // check that the D0 passes the cuts
 226   // (if we have a D*+, it has to pass as D0,
 227   //  if we have a D*-, it has to pass as D0bar)
 228
 229   if(testD0) {
 230     Int_t okD0=0,okD0bar=0;
 231     Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
 232     if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE;
 233   }
 234
 235   if( (PtProng(0)<cutsDstar[2]) || (PtProng(0)>cutsDstar[3]) ) return kFALSE;
 236
 237   Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 238   Double_t invmDstar = InvMassDstarKpipi();
 239   if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;
 240
 241   Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 242   if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;
 243
 244   Double_t theta = AngleD0dkpPisoft();
 245   if(theta>cutsDstar[4]) return kFALSE;
 246
 247   return kTRUE;
 248 }
 249 //-----------------------------------------------------------------------------
 250 Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0,
 251                                          Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const
 252 {
 253   // cuts on Lambdac candidates to V0+bachelor
 254   // (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
 255   // 0 = inv. mass half width in K0s hypothesis [GeV]
 256   // 1 = inv. mass half width in Lambda hypothesis [GeV]
 257   // 2 = inv. mass V0 in K0s hypothesis half width [GeV]
 258   // 3 = inv. mass V0 in Lambda hypothesis half width [GeV]
 259   // 4 = pT min Bachelor track [GeV/c]
 260   // 5 = pT min V0-Positive track [GeV/c]
 261   // 6 = pT min V0-Negative track [GeV/c]
 262   // 7 = dca cut on the cascade (cm)
 263   // 8 = dca cut on the V0 (cm)
 264
 265   //   if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() )
 266   //     { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
 267
 268   Double_t mLck0sp,mLcLpi;
 269   okLck0sp=1; okLcLpi=1; okLcLbarpi=1;
 270
 271   Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
 272   Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
 273   Double_t mLPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
 274
 275   // k0s + p
 276   double mk0s = Getv0()->MassK0Short();
 277   mLck0sp = InvMassLctoK0sP();
 278
 279   // lambda + pi
 280   double mlambda = Getv0()->MassLambda();
 281   double malambda = Getv0()->MassAntiLambda();
 282   mLcLpi = InvMassLctoLambdaPi();
 283
 284   // cut on Lc mass
 285   //   with k0s p hypothesis
 286   if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
 287   //   with Lambda pi hypothesis
 288   if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
 289   okLcLbarpi = okLcLpi;
 290
 291   // cuts on the v0 mass
 292   if( TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
 293   //if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] &&
 294   //TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
 295   if( !(GetBachelor()->Charge()==+1 && TMath::Abs(mlambda-mLPDG)<=cutsLctoV0[3]) ) okLcLpi = 0;
 296   if( !(GetBachelor()->Charge()==-1 && TMath::Abs(malambda-mLPDG)<=cutsLctoV0[3]) ) okLcLbarpi = 0;
 297
 298   if(!okLck0sp && !okLcLpi && !okLcLbarpi) return 0;
 299
 300   // cuts on the minimum pt of the tracks
 301   if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
 302   if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
 303   if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
 304
 305   // cut on the cascade dca
 306   if( TMath::Abs(GetDCA(0))>cutsLctoV0[7] //||
 307       //TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[7] ||
 308       //TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[7]
 309       ) return 0;
 310
 311   // cut on the v0 dca
 312   if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[8]) return 0;
 313
 314   // cut on V0 cosine of pointing angle wrt PV
 315   if (CosV0PointingAngle() < cutsLctoV0[9]) { // cosine of V0 pointing angle wrt primary vertex
 316     AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
 317     return 0;
 318   }
 319
 320   // cut on bachelor transverse impact parameter wrt PV
 321   if (TMath::Abs(Getd0Prong(0)) > cutsLctoV0[10]) { // bachelor transverse impact parameter wrt PV
 322     AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
 323     return 0;
 324   }
 325
 326   // cut on V0 transverse impact parameter wrt PV
 327   if (TMath::Abs(Getd0Prong(1)) > cutsLctoV0[11]) { // V0 transverse impact parameter wrt PV
 328     AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
 329     return 0;
 330   }
 331
 332   // cut on K0S invariant mass veto
 333   if (TMath::Abs(Getv0()->MassK0Short()-mk0sPDG) < cutsLctoV0[12]) { // K0S invariant mass veto
 334     AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
 335     return 0;
 336   }
 337
 338   // cut on Lambda/LambdaBar invariant mass veto
 339   if (TMath::Abs(Getv0()->MassLambda()-mLPDG) < cutsLctoV0[13] ||
 340       TMath::Abs(Getv0()->MassAntiLambda()-mLPDG) < cutsLctoV0[13] ) { // Lambda/LambdaBar invariant mass veto
 341     AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
 342     return 0;
 343   }
 344
 345   // cut on gamma invariant mass veto
 346   if (Getv0()->InvMass2Prongs(0,1,11,11) < cutsLctoV0[14]) { // K0S invariant mass veto
 347     AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
 348     return 0;
 349   }
 350
 351   // cut on V0 pT min
 352   if (Getv0()->Pt() < cutsLctoV0[15]) { // V0 pT min
 353     AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",Getv0()->Pt(),cutsLctoV0[15]));
 354     return 0;
 355   }
 356
 357   return true;
 358
 359 }
 360 //-----------------------------------------------------------------------------
 361 Double_t AliAODRecoCascadeHF::AngleD0dkpPisoft() const {
 362   //
 363   // Angle of soft pion to D0 decay plane
 364   //
 365
 366   TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
 367   TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
 368   TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
 369
 370   TVector3 perp = p3Trk0.Cross(p3Trk1);
 371   Double_t theta = p3Trk2.Angle(perp);
 372   if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
 373   theta = TMath::Pi()/2. - theta;
 374
 375   return theta;
 376 }
 377 //-----------------------------------------------------------------------------
 378 Bool_t AliAODRecoCascadeHF::TrigonometricalCut() const {
 379   //
 380   // Trigonometrical constraint
 381   //
 382   TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
 383   TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
 384   TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
 385
 386   Double_t alpha = p3Trk0.Angle(p3Trk2);
 387   Double_t beta = p3Trk1.Angle(p3Trk2);
 388
 389   Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(AngleD0dkpPisoft());
 390   Double_t cosphi02 = TMath::Cos(beta) / TMath::Cos(AngleD0dkpPisoft());
 391
 392   Double_t phi01 = TMath::ACos(cosphi01);
 393   Double_t phi02 = TMath::ACos(cosphi02);
 394   Double_t phi00 = p3Trk0.Angle(p3Trk1);
 395
 396   if((phi01>phi00) || (phi02>phi00)) return kFALSE;
 397   return kTRUE;
 398 }
 399
 400 //-----------------------------------------------------------------------------
 401 Double_t AliAODRecoCascadeHF::DecayLengthV0() const
 402 {
 403   //
 404   // Returns V0 decay length wrt primary vertex
 405   //
 406
 407   AliAODv0 *v0 = (AliAODv0*)Getv0();
 408
 409   if (!v0)
 410     return -1.;
 411   AliAODVertex *vtxPrimary = GetPrimaryVtx();
 412   Double_t posVtx[3] = {0.,0.,0.};
 413   vtxPrimary->GetXYZ(posVtx);
 414   return v0->DecayLengthV0(posVtx);
 415
 416 }
 417 //-----------------------------------------------------------------------------
 418 Double_t AliAODRecoCascadeHF::DecayLengthXYV0() const
 419 {
 420   //
 421   // Returns transverse V0 decay length wrt primary vertex
 422   //
 423   AliAODv0 *v0 = (AliAODv0*)Getv0();
 424
 425   if (!v0)
 426     return -1.;
 427   AliAODVertex *vtxPrimary = GetPrimaryVtx();
 428   Double_t posVtx[3] = {0.,0.,0.};
 429   vtxPrimary->GetXYZ(posVtx);
 430   return v0->DecayLengthXY(posVtx);
 431
 432 }
 433 //-----------------------------------------------------------------------------
 434 Double_t AliAODRecoCascadeHF::CosV0PointingAngle() const
 435 {
 436   //
 437   // Returns cosine of V0 pointing angle wrt primary vertex
 438   //
 439
 440   AliAODv0 *v0 = (AliAODv0*)Getv0();
 441
 442   if (!v0)
 443     return -999.;
 444
 445   AliAODVertex *vtxPrimary = GetPrimaryVtx();
 446   Double_t posVtx[3] = {0.,0.,0.};
 447   vtxPrimary->GetXYZ(posVtx);
 448   return v0->CosPointingAngle(posVtx);
 449
 450 }
 451 //-----------------------------------------------------------------------------
 452 Double_t AliAODRecoCascadeHF::CosV0PointingAngleXY() const
 453 {
 454   //
 455   // Returns XY cosine of V0 pointing angle wrt primary vertex
 456   //
 457
 458   AliAODv0 *v0 = (AliAODv0*)Getv0();
 459
 460   if (!v0)
 461     return -999.;
 462
 463   AliAODVertex *vtxPrimary = GetPrimaryVtx();
 464   Double_t posVtx[3] = {0.,0.,0.};
 465   vtxPrimary->GetXYZ(posVtx);
 466   return v0->CosPointingAngleXY(posVtx);
 467
 468 }
 469 //-----------------------------------------------------------------------------
 470 Double_t AliAODRecoCascadeHF::NormalizedV0DecayLength() const
 471 {
 472   //
 473   // Returns V0 normalized decay length wrt primary vertex
 474   //
 475
 476   AliAODv0 *v0 = (AliAODv0*)Getv0();
 477
 478   if (!v0)
 479     return -1.;
 480   //AliAODVertex *vtxPrimary = GetPrimaryVtx();
 481   //Double_t posVtx[3] = {0.,0.,0.};
 482   //vtxPrimary->GetXYZ(posVtx);
 483   //return v0->NormalizedDecayLength(posVtx);
 484   return v0->NormalizedDecayLength(GetPrimaryVtx());
 485
 486 }
 487 //-----------------------------------------------------------------------------
 488 Double_t AliAODRecoCascadeHF::NormalizedV0DecayLengthXY() const
 489 {
 490   //
 491   // Returns transverse V0 normalized decay length wrt primary vertex
 492   //
 493   AliAODv0 *v0 = (AliAODv0*)Getv0();
 494
 495   if (!v0)
 496     return -1.;
 497   //AliAODVertex *vtxPrimary = GetPrimaryVtx();
 498   //Double_t posVtx[3] = {0.,0.,0.};
 499   //vtxPrimary->GetXYZ(posVtx);
 500   //return v0->NormalizedDecayLengthXY(posVtx);
 501   return v0->NormalizedDecayLengthXY(GetPrimaryVtx());
 502
 503 }