1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // Class for AOD reconstructed heavy-flavour cascades
22 // Author: X-M. Zhang, zhangxm@iopp.ccnu.edu.cn
23 /////////////////////////////////////////////////////////////
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"
34 ClassImp(AliAODRecoCascadeHF)
35 //-----------------------------------------------------------------------------
37 AliAODRecoCascadeHF::AliAODRecoCascadeHF() :
38 AliAODRecoDecayHF2Prong()
41 // Default Constructor
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)
51 // Constructor with AliAODVertex for decay vertex
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)
61 // Constructor with decay vertex and without prongs momenta
65 //-----------------------------------------------------------------------------
66 AliAODRecoCascadeHF::AliAODRecoCascadeHF(const AliAODRecoCascadeHF &source) :
67 AliAODRecoDecayHF2Prong(source)
73 //-----------------------------------------------------------------------------
74 AliAODRecoCascadeHF &AliAODRecoCascadeHF::operator=(const AliAODRecoCascadeHF &source)
77 // assignment operator
79 if(&source == this) return *this;
81 AliAODRecoDecayHF2Prong::operator=(source);
85 //-----------------------------------------------------------------------------
86 AliAODRecoCascadeHF::~AliAODRecoCascadeHF()
92 //-----------------------------------------------------------------------------
93 Double_t AliAODRecoCascadeHF::InvMassDstarKpipi() const
96 // 3 prong invariant mass of the D0 daughters and the soft pion
100 e[0]=Get2Prong()->EProng(0,211);
101 e[1]=Get2Prong()->EProng(1,321);
103 e[0]=Get2Prong()->EProng(0,321);
104 e[1]=Get2Prong()->EProng(1,211);
108 Double_t esum = e[0]+e[1]+e[2];
109 Double_t minv = TMath::Sqrt(esum*esum-P()*P());
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
119 // Check if this candidate is matched to a MC signal
121 // If yes, return label (>=0) of the AliAODMCParticle
124 Int_t ndg=GetNDaughters();
126 AliError("No daughters available");
130 Int_t lab2Prong = -1;
133 AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
134 lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
136 AliAODv0 *theV0 = Getv0();
137 lab2Prong = theV0->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
140 if(lab2Prong<0) return -1;
142 Int_t dgLabels[10]={0,0,0,0,0,0,0,0,0,0};
144 // loop on daughters and write labels
145 for(Int_t i=0; i<ndg; i++) {
146 AliVTrack *trk = (AliVTrack*)GetDaughter(i);
147 Int_t lab = trk->GetLabel();
148 if(lab==-1) { // this daughter is the 2prong
150 } else if(lab<-1) continue;
154 Int_t finalLabel = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
157 // debug printouts for Lc->V0 bachelor case
159 if ( isV0 && (dgLabels[0]!=-1 && dgLabels[1]!=-1) ) {
160 AliAODv0 *theV0 = Getv0();
161 Bool_t onTheFly = theV0->GetOnFlyStatus();
163 if ( (pdgDg[1]==2212 && pdgDg[0]==310) ||
164 (pdgDg[1]==211 && pdgDg[0]==3122) ) {
165 Int_t pdgDgtemp[2]={pdgDg[1],pdgDg[0]};
166 pdgDg[0]=pdgDgtemp[0];
167 pdgDg[1]=pdgDgtemp[1];
170 if (pdgDg[0]==2212 && pdgDg[1]==310) {
171 AliAODMCParticle*k0s = (AliAODMCParticle*)mcArray->At(lab2Prong);
172 Int_t labK0 = k0s->GetMother();
173 AliAODMCParticle*k0bar = (AliAODMCParticle*)mcArray->At(labK0);
174 AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d) -> LabelK0S=%d (%d -> %d %d)",onTheFly,labK0,k0bar->GetPdgCode(),lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
175 AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
177 dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
178 } else if (pdgDg[0]==211 && pdgDg[1]==3122) {
179 AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d -> %d %d)",onTheFly,lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
180 AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
182 dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
191 //-----------------------------------------------------------------------------
192 Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
193 const Double_t *cutsD0,
197 // cutsDstar[0] = inv. mass half width of D* [GeV]
198 // cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
199 // cutsDstar[2] = PtMin of pi_s [GeV/c]
200 // cutsDstar[3] = PtMax of pi_s [GeV/c]
201 // cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
203 // cutsD0[0] = inv. mass half width [GeV]
204 // cutsD0[1] = dca [cm]
205 // cutsD0[2] = cosThetaStar
206 // cutsD0[3] = pTK [GeV/c]
207 // cutsD0[4] = pTPi [GeV/c]
208 // cutsD0[5] = d0K [cm] upper limit!
209 // cutsD0[6] = d0Pi [cm] upper limit!
210 // cutsD0[7] = d0d0 [cm^2]
211 // cutsD0[8] = cosThetaPoint
214 // check that the D0 passes the cuts
215 // (if we have a D*+, it has to pass as D0,
216 // if we have a D*-, it has to pass as D0bar)
219 Int_t okD0=0,okD0bar=0;
220 Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
221 if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE;
224 if( (PtProng(0)<cutsDstar[2]) || (PtProng(0)>cutsDstar[3]) ) return kFALSE;
226 Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
227 Double_t invmDstar = InvMassDstarKpipi();
228 if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;
230 Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
231 if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;
233 Double_t theta = AngleD0dkpPisoft();
234 if(theta>cutsDstar[4]) return kFALSE;
238 //-----------------------------------------------------------------------------
239 Bool_t AliAODRecoCascadeHF::SelectLctoV0(const Double_t *cutsLctoV0,
240 Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const
242 // cuts on Lambdac candidates to V0+bachelor
243 // (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
244 // 0 = inv. mass half width in K0s hypothesis [GeV]
245 // 1 = inv. mass half width in Lambda hypothesis [GeV]
246 // 2 = inv. mass V0 in K0s hypothesis half width [GeV]
247 // 3 = inv. mass V0 in Lambda hypothesis half width [GeV]
248 // 4 = pT min Bachelor track [GeV/c]
249 // 5 = pT min V0-Positive track [GeV/c]
250 // 6 = pT min V0-Negative track [GeV/c]
251 // 7 = dca cut on the cascade (cm)
252 // 8 = dca cut on the V0 (cm)
254 // if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() )
255 // { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
257 Double_t mLck0sp,mLcLpi;
258 okLck0sp=1; okLcLpi=1; okLcLbarpi=1;
260 Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
261 Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
262 Double_t mLPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
265 double mk0s = Getv0()->MassK0Short();
266 mLck0sp = InvMassLctoK0sP();
269 double mlambda = Getv0()->MassLambda();
270 double malambda = Getv0()->MassAntiLambda();
271 mLcLpi = InvMassLctoLambdaPi();
274 // with k0s p hypothesis
275 if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
276 // with Lambda pi hypothesis
277 if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
278 okLcLbarpi = okLcLpi;
280 // cuts on the v0 mass
281 if( TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
282 //if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] &&
283 //TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
284 if( !(GetBachelor()->Charge()==+1 && TMath::Abs(mlambda-mLPDG)<=cutsLctoV0[3]) ) okLcLpi = 0;
285 if( !(GetBachelor()->Charge()==-1 && TMath::Abs(malambda-mLPDG)<=cutsLctoV0[3]) ) okLcLbarpi = 0;
287 if(!okLck0sp && !okLcLpi && !okLcLbarpi) return 0;
289 // cuts on the minimum pt of the tracks
290 if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
291 if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
292 if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
294 // cut on the cascade dca
295 if( TMath::Abs(GetDCA(0))>cutsLctoV0[7] //||
296 //TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[7] ||
297 //TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[7]
301 if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[8]) return 0;
303 // cut on V0 cosine of pointing angle wrt PV
304 if (CosV0PointingAngle() < cutsLctoV0[9]) { // cosine of V0 pointing angle wrt primary vertex
305 AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
309 // cut on bachelor transverse impact parameter wrt PV
310 if (TMath::Abs(Getd0Prong(0)) > cutsLctoV0[10]) { // bachelor transverse impact parameter wrt PV
311 AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
315 // cut on V0 transverse impact parameter wrt PV
316 if (TMath::Abs(Getd0Prong(1)) > cutsLctoV0[11]) { // V0 transverse impact parameter wrt PV
317 AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
321 // cut on K0S invariant mass veto
322 if (TMath::Abs(Getv0()->MassK0Short()-mk0sPDG) < cutsLctoV0[12]) { // K0S invariant mass veto
323 AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
327 // cut on Lambda/LambdaBar invariant mass veto
328 if (TMath::Abs(Getv0()->MassLambda()-mLPDG) < cutsLctoV0[13] ||
329 TMath::Abs(Getv0()->MassAntiLambda()-mLPDG) < cutsLctoV0[13] ) { // Lambda/LambdaBar invariant mass veto
330 AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
334 // cut on gamma invariant mass veto
335 if (Getv0()->InvMass2Prongs(0,1,11,11) < cutsLctoV0[14]) { // K0S invariant mass veto
336 AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
341 if (Getv0()->Pt() < cutsLctoV0[15]) { // V0 pT min
342 AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",Getv0()->Pt(),cutsLctoV0[15]));
349 //-----------------------------------------------------------------------------
350 Double_t AliAODRecoCascadeHF::AngleD0dkpPisoft() const {
352 // Angle of soft pion to D0 decay plane
355 TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
356 TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
357 TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
359 TVector3 perp = p3Trk0.Cross(p3Trk1);
360 Double_t theta = p3Trk2.Angle(perp);
361 if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
362 theta = TMath::Pi()/2. - theta;
366 //-----------------------------------------------------------------------------
367 Bool_t AliAODRecoCascadeHF::TrigonometricalCut() const {
369 // Trigonometrical constraint
371 TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
372 TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
373 TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
375 Double_t alpha = p3Trk0.Angle(p3Trk2);
376 Double_t beta = p3Trk1.Angle(p3Trk2);
378 Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(AngleD0dkpPisoft());
379 Double_t cosphi02 = TMath::Cos(beta) / TMath::Cos(AngleD0dkpPisoft());
381 Double_t phi01 = TMath::ACos(cosphi01);
382 Double_t phi02 = TMath::ACos(cosphi02);
383 Double_t phi00 = p3Trk0.Angle(p3Trk1);
385 if((phi01>phi00) || (phi02>phi00)) return kFALSE;
389 //-----------------------------------------------------------------------------
390 Double_t AliAODRecoCascadeHF::DecayLengthV0() const
393 // Returns V0 decay length wrt primary vertex
396 AliAODv0 *v0 = (AliAODv0*)Getv0();
400 AliAODVertex *vtxPrimary = GetPrimaryVtx();
401 Double_t posVtx[3] = {0.,0.,0.};
402 vtxPrimary->GetXYZ(posVtx);
403 return v0->DecayLengthV0(posVtx);
406 //-----------------------------------------------------------------------------
407 Double_t AliAODRecoCascadeHF::DecayLengthXYV0() const
410 // Returns transverse V0 decay length wrt primary vertex
412 AliAODv0 *v0 = (AliAODv0*)Getv0();
416 AliAODVertex *vtxPrimary = GetPrimaryVtx();
417 Double_t posVtx[3] = {0.,0.,0.};
418 vtxPrimary->GetXYZ(posVtx);
419 return v0->DecayLengthXY(posVtx);
422 //-----------------------------------------------------------------------------
423 Double_t AliAODRecoCascadeHF::CosV0PointingAngle() const
426 // Returns cosine of V0 pointing angle wrt primary vertex
429 AliAODv0 *v0 = (AliAODv0*)Getv0();
434 AliAODVertex *vtxPrimary = GetPrimaryVtx();
435 Double_t posVtx[3] = {0.,0.,0.};
436 vtxPrimary->GetXYZ(posVtx);
437 return v0->CosPointingAngle(posVtx);