PWG3/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 /////////////////////////////////////////////////////////////
  17 //
  18 // Class for AOD reconstructed heavy-flavour cascades
  19 //
  20 // Author: X-M. Zhang, zhangxm@iopp.ccnu.edu.cn
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TVector3.h>
  24 #include <TDatabasePDG.h>
  25 #include <TClonesArray.h>
  26 #include "AliAODMCParticle.h"
  27 #include "AliAODRecoDecay.h"
  28 #include "AliAODVertex.h"
  29 #include "AliAODRecoDecayHF2Prong.h"
  30 #include "AliAODRecoCascadeHF.h"
  31
  32 ClassImp(AliAODRecoCascadeHF)
  33 //-----------------------------------------------------------------------------
  34
  35 AliAODRecoCascadeHF::AliAODRecoCascadeHF() :
  36   AliAODRecoDecayHF2Prong()
  37 {
  38   //
  39   // Default Constructor
  40   //
  41 }
  42 //-----------------------------------------------------------------------------
  43 AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
  44                                          Double_t *px, Double_t *py, Double_t *pz,
  45                                          Double_t *d0, Double_t *d0err, Double_t dca) :
  46   AliAODRecoDecayHF2Prong(vtx2, px, py, pz, d0, d0err, dca)
  47 {
  48   //
  49   //  Constructor with AliAODVertex for decay vertex
  50   //
  51   SetCharge(charge);
  52 }
  53 //-----------------------------------------------------------------------------
  54 AliAODRecoCascadeHF::AliAODRecoCascadeHF(AliAODVertex *vtx2, Short_t charge,
  55                                          Double_t *d0, Double_t *d0err, Double_t dca) :
  56   AliAODRecoDecayHF2Prong(vtx2, d0, d0err, dca)
  57 {
  58   //
  59   //  Constructor with decay vertex and without prongs momenta
  60   //
  61   SetCharge(charge);
  62 }
  63 //-----------------------------------------------------------------------------
  64 AliAODRecoCascadeHF::AliAODRecoCascadeHF(const AliAODRecoCascadeHF &source) :
  65   AliAODRecoDecayHF2Prong(source)
  66 {
  67   //
  68   // Copy constructor
  69   //
  70 }
  71 //-----------------------------------------------------------------------------
  72 AliAODRecoCascadeHF &AliAODRecoCascadeHF::operator=(const AliAODRecoCascadeHF &source)
  73 {
  74   //
  75   // assignment operator
  76   //
  77   if(&source == this) return *this;
  78
  79   AliAODRecoDecayHF2Prong::operator=(source);
  80
  81   return *this;
  82 }
  83 //-----------------------------------------------------------------------------
  84 AliAODRecoCascadeHF::~AliAODRecoCascadeHF()
  85 {
  86   //
  87   // Default Destructor
  88   //
  89 }
  90 //-----------------------------------------------------------------------------
  91 Double_t AliAODRecoCascadeHF::InvMassDstarKpipi() const
  92 {
  93   //
  94   // 3 prong invariant mass of the D0 daughters and the soft pion
  95   //
  96   Double_t e[3];
  97   if (Charge()>0){
  98     e[0]=Get2Prong()->EProng(0,211);
  99     e[1]=Get2Prong()->EProng(1,321);
 100   }else{
 101     e[0]=Get2Prong()->EProng(0,321);
 102     e[1]=Get2Prong()->EProng(1,211);
 103   }
 104   e[2]=EProng(0,211);
 105
 106   Double_t esum = e[0]+e[1]+e[2];
 107   Double_t minv = TMath::Sqrt(esum*esum-P()*P());
 108
 109   return minv;
 110 }
 111 //----------------------------------------------------------------------------
 112 Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
 113                                      Int_t *pdgDg,Int_t *pdgDg2prong,
 114                                      TClonesArray *mcArray) const
 115 {
 116   //
 117   // Check if this candidate is matched to a MC signal
 118   // If no, return -1
 119   // If yes, return label (>=0) of the AliAODMCParticle
 120   //
 121
 122   Int_t ndg=GetNDaughters();
 123   if(!ndg) {
 124     AliError("No daughters available");
 125     return -1;
 126   }
 127
 128   AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
 129   Int_t lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
 130
 131   if(lab2Prong<0) return -1;
 132
 133   Int_t dgLabels[10];
 134
 135   // loop on daughters and write labels
 136   for(Int_t i=0; i<ndg; i++) {
 137     AliVTrack *trk = (AliVTrack*)GetDaughter(i);
 138     Int_t lab = trk->GetLabel();
 139     if(lab==-1) { // this daughter is the 2prong
 140       lab=lab2Prong;
 141     } else if(lab<-1) {
 142       printf("daughter with negative label\n");
 143       continue;
 144     }
 145     dgLabels[i] = lab;
 146   }
 147
 148   return AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
 149 }
 150 //-----------------------------------------------------------------------------
 151 Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
 152                                         const Double_t *cutsD0,
 153                                         Bool_t testD0) const
 154 {
 155   //
 156   // cutsDstar[0] = inv. mass half width of D* [GeV]
 157   // cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
 158   // cutsDstar[2] = PtMin of pi_s [GeV/c]
 159   // cutsDstar[3] = PtMax of pi_s [GeV/c]
 160   // cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
 161   //
 162   // cutsD0[0] = inv. mass half width [GeV]
 163   // cutsD0[1] = dca [cm]
 164   // cutsD0[2] = cosThetaStar
 165   // cutsD0[3] = pTK [GeV/c]
 166   // cutsD0[4] = pTPi [GeV/c]
 167   // cutsD0[5] = d0K [cm]   upper limit!
 168   // cutsD0[6] = d0Pi [cm]  upper limit!
 169   // cutsD0[7] = d0d0 [cm^2]
 170   // cutsD0[8] = cosThetaPoint
 171
 172
 173   // check that the D0 passes the cuts
 174   // (if we have a D*+, it has to pass as D0,
 175   //  if we have a D*-, it has to pass as D0bar)
 176
 177   if(testD0) {
 178     Int_t okD0=0,okD0bar=0;
 179     Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
 180     if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE;
 181   }
 182
 183   if( (PtProng(0)<cutsDstar[2]) || (PtProng(0)>cutsDstar[3]) ) return kFALSE;
 184
 185   Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 186   Double_t invmDstar = InvMassDstarKpipi();
 187   if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;
 188
 189   Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 190   if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;
 191
 192   Double_t theta = AngleD0dkpPisoft();
 193   if(theta>cutsDstar[4]) return kFALSE;
 194
 195   return kTRUE;
 196 }
 197 //-----------------------------------------------------------------------------
 198 Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0,
 199                                          Bool_t okLck0sp, Bool_t okLcLpi) const
 200 {
 201   // cuts on Lambdac candidates to V0+bachelor
 202   // (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
 203   // 0 = inv. mass half width in K0s hypothesis [GeV]
 204   // 1 = inv. mass half width in Lambda hypothesis [GeV]
 205   // 2 = inv. mass V0 in K0s hypothesis half width [GeV]
 206   // 3 = inv. mass V0 in Lambda hypothesis half width [GeV]
 207   // 4 = pT min Bachelor track [GeV/c]
 208   // 5 = pT min V0-Positive track [GeV/c]
 209   // 6 = pT min V0-Negative track [GeV/c]
 210   // 7 = dca cut on the V0 (cm)
 211   // 8 = dca cut on the cascade (cm)
 212
 213 //   if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() )
 214 //     { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
 215
 216   Double_t mLck0sp,mLcLpi;
 217   okLck0sp=1; okLcLpi=1;
 218
 219   Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
 220   Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
 221   Double_t mLPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
 222
 223   // k0s + p
 224   double mk0s = Getv0()->MassK0Short();
 225   mLck0sp = InvMassLctoK0sP();
 226
 227   // lambda + pi
 228   double mlambda = Getv0()->MassLambda();
 229   double malambda = Getv0()->MassAntiLambda();
 230   mLcLpi = InvMassLctoLambdaPi();
 231
 232   // cut on Lc mass
 233   //   with k0s p hypothesis
 234   if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
 235   //   with Lambda pi hypothesis
 236   if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
 237
 238   // cuts on the v0 mass
 239   if(TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
 240   if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] &&
 241       TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
 242
 243   if(!okLck0sp && !okLcLpi) return 0;
 244
 245   // cuts on the minimum pt of the tracks
 246   if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
 247   if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
 248   if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
 249
 250   // cut on the v0 dca
 251   if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[7]) return 0;
 252
 253   // cut on the cascade dca
 254   if( TMath::Abs(GetDCA(0))>cutsLctoV0[8] ||
 255       TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[8] ||
 256       TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[8] ) return 0;
 257
 258   return true;
 259
 260 }
 261 //-----------------------------------------------------------------------------
 262 Double_t AliAODRecoCascadeHF::AngleD0dkpPisoft() const {
 263   //
 264   // Angle of soft pion to D0 decay plane
 265   //
 266
 267   TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
 268   TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
 269   TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
 270
 271   TVector3 perp = p3Trk0.Cross(p3Trk1);
 272   Double_t theta = p3Trk2.Angle(perp);
 273   if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
 274   theta = TMath::Pi()/2. - theta;
 275
 276   return theta;
 277 }
 278 //-----------------------------------------------------------------------------
 279 Bool_t AliAODRecoCascadeHF::TrigonometricalCut() const {
 280   //
 281   // Trigonometrical constraint
 282   //
 283   TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
 284   TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
 285   TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
 286
 287   Double_t alpha = p3Trk0.Angle(p3Trk2);
 288   Double_t beta = p3Trk1.Angle(p3Trk2);
 289
 290   Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(AngleD0dkpPisoft());
 291   Double_t cosphi02 = TMath::Cos(beta) / TMath::Cos(AngleD0dkpPisoft());
 292
 293   Double_t phi01 = TMath::ACos(cosphi01);
 294   Double_t phi02 = TMath::ACos(cosphi02);
 295   Double_t phi00 = p3Trk0.Angle(p3Trk1);
 296
 297   if((phi01>phi00) || (phi02>phi00)) return kFALSE;
 298   return kTRUE;
 299 }