AliPhysics  323987d (323987d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAODRecoCascadeHF.cxx
Go to the documentation of this file.
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 
19 //
20 // Class for AOD reconstructed heavy-flavour cascades
21 //
22 // Author: X-M. Zhang, zhangxm@iopp.ccnu.edu.cn
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"
32 #include "AliAODRecoCascadeHF.h"
33 #include "AliRDHFCuts.h"
34 
38 
39 //-----------------------------------------------------------------------------
40 
43 {
47 }
48 //-----------------------------------------------------------------------------
50  Double_t *px, Double_t *py, Double_t *pz,
51  Double_t *d0, Double_t *d0err, Double_t dca) :
52  AliAODRecoDecayHF2Prong(vtx2, px, py, pz, d0, d0err, dca)
53 {
57  SetCharge(charge);
58 }
59 //-----------------------------------------------------------------------------
61  Double_t *d0, Double_t *d0err, Double_t dca) :
62  AliAODRecoDecayHF2Prong(vtx2, d0, d0err, dca)
63 {
67  SetCharge(charge);
68 }
69 //-----------------------------------------------------------------------------
72 {
76 }
77 //-----------------------------------------------------------------------------
79 {
80  //
82  //
83  if(&source == this) return *this;
84 
86 
87  return *this;
88 }
89 //-----------------------------------------------------------------------------
91 {
95 }
96 //-----------------------------------------------------------------------------
98 {
102  Double_t e[3];
103  if (Charge()>0){
104  e[0]=Get2Prong()->EProng(0,211);
105  e[1]=Get2Prong()->EProng(1,321);
106  }else{
107  e[0]=Get2Prong()->EProng(0,321);
108  e[1]=Get2Prong()->EProng(1,211);
109  }
110  e[2]=EProng(0,211);
111 
112  Double_t esum = e[0]+e[1]+e[2];
113  Double_t minv = TMath::Sqrt(esum*esum-P()*P());
114 
115  return minv;
116 }
117 //----------------------------------------------------------------------------
119  Int_t *pdgDg,Int_t *pdgDg2prong,
120  TClonesArray *mcArray, Bool_t isV0) const
121 {
127 
128  Int_t ndg=GetNDaughters();
129  if(ndg==0) {
130  AliError("No daughters available");
131  return -1;
132  }
133 
134  if ( isV0 &&
135  ( (pdgDg[1]==2212 && pdgDg[0]==310) ||
136  (pdgDg[1]==211 && pdgDg[0]==3122) ||
137  (pdgDg[1]==211 && pdgDg[0]==310) ||
138  (pdgDg[1]==321 && pdgDg[0]==310) ) ) {
139  AliWarning(Form("Please, pay attention: first element in AliAODRecoCascadeHF object must be the bachelor and second one V0. Skipping! (pdgDg[0] = %d, (pdgDg[1] = %d)", pdgDg[0], pdgDg[1]));
140  return -1;
141  }
142 
143  Int_t lab2Prong = -1;
144 
145  if (!isV0) {
146  AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
147  lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong);
148  } else {
149  AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
150  lab2Prong = theV0->MatchToMC(pdgabs2prong,mcArray,2,pdgDg2prong); // the V0
151  }
152 
153  if(lab2Prong<0) return -1;
154 
155  Int_t dgLabels[10]={0,0,0,0,0,0,0,0,0,0};
156 
157  if (!isV0) {
158  // loop on daughters and write labels
159  for(Int_t i=0; i<ndg; i++) {
160  AliVTrack *trk = dynamic_cast<AliVTrack*>(GetDaughter(i));
161  if(!trk) continue;
162  Int_t lab = trk->GetLabel();
163  if(lab==-1) { // this daughter is the 2prong
164  lab=lab2Prong;
165  } else if(lab<-1) continue;
166  dgLabels[i] = lab;
167  }
168  } else {
169  AliVTrack *trk = dynamic_cast<AliVTrack*>(GetBachelor()); // the bachelor
170  if (!trk) return -1;
171  dgLabels[0] = trk->GetLabel();//TMath::Abs(trk->GetLabel());
172  dgLabels[1] = lab2Prong;
173  }
174 
175  Int_t finalLabel = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels,2,2,pdgDg);
176 
177  if (finalLabel>=0) {
178  // Debug printouts for Lc->V0 bachelor, D+->K0s+pi, Ds->K0s+K cases
179 
180  if ( isV0 && (dgLabels[0]!=-1 && dgLabels[1]!=-1) ) {
181  AliAODv0 *theV0 = dynamic_cast<AliAODv0*>(Getv0());
182  Bool_t onTheFly = theV0->GetOnFlyStatus();
183  if (pdgDg[1]==310 && (pdgDg[0]==2212 || pdgDg[0]==211 || pdgDg[0]==321)) {
184  AliAODMCParticle*k0s = dynamic_cast<AliAODMCParticle*>(mcArray->At(lab2Prong));
185  if (k0s) {
186  Int_t labK0 = k0s->GetMother();
187  AliAODMCParticle*k0bar = dynamic_cast<AliAODMCParticle*>(mcArray->At(labK0));
188  if (k0bar) {
189  AliDebug(1, Form(" (onTheFly=%1d) LabelV0=%d (%d) -> LabelK0S=%d (%d -> %d %d)",
190  onTheFly, labK0, k0bar->GetPdgCode(), lab2Prong, pdgabs2prong, pdgDg2prong[0], pdgDg2prong[1]));
191  AliDebug(1, Form(" LabelCandidate=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
192  finalLabel, pdgabs, dgLabels[0], pdgDg[0], dgLabels[1], pdgDg[1]));
193  }
194  }
195  } else if (pdgDg[0]==211 && pdgDg[1]==3122) {
196  AliDebug(1,Form(" (onTheFly=%1d) LabelV0=%d (%d -> %d %d)",onTheFly,lab2Prong,pdgabs2prong,pdgDg2prong[0],pdgDg2prong[1]));
197  AliDebug(1,Form(" LabelLc=%d (%d) -> LabelBachelor=%d (%d) LabelV0=%d (%d)",
198  finalLabel,pdgabs,
199  dgLabels[0],pdgDg[0],dgLabels[1],pdgDg[1]));
200  }
201  }
202 
203  }
204 
205  return finalLabel;
206 
207 }
208 //-----------------------------------------------------------------------------
210  const Double_t *cutsD0,
211  Bool_t testD0) const
212 {
213  //
214  // cutsDstar[0] = inv. mass half width of D* [GeV]
215  // cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
216  // cutsDstar[2] = PtMin of pi_s [GeV/c]
217  // cutsDstar[3] = PtMax of pi_s [GeV/c]
218  // cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
219  //
220  // cutsD0[0] = inv. mass half width [GeV]
221  // cutsD0[1] = dca [cm]
222  // cutsD0[2] = cosThetaStar
223  // cutsD0[3] = pTK [GeV/c]
224  // cutsD0[4] = pTPi [GeV/c]
225  // cutsD0[5] = d0K [cm] upper limit!
226  // cutsD0[6] = d0Pi [cm] upper limit!
227  // cutsD0[7] = d0d0 [cm^2]
228  // cutsD0[8] = cosThetaPoint
229 
230 
231  // check that the D0 passes the cuts
232  // (if we have a D*+, it has to pass as D0,
233  // if we have a D*-, it has to pass as D0bar)
234 
235  if(testD0) {
236  Int_t okD0=0,okD0bar=0;
237  Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
238  if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE;
239  }
240 
241  if( (PtProng(0)<cutsDstar[2]) || (PtProng(0)>cutsDstar[3]) ) return kFALSE;
242 
243  Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
244  Double_t invmDstar = InvMassDstarKpipi();
245  if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;
246 
247  Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
248  if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;
249 
250  Double_t theta = AngleD0dkpPisoft();
251  if(theta>cutsDstar[4]) return kFALSE;
252 
253  return kTRUE;
254 }
255 //-----------------------------------------------------------------------------
257  Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const
258 {
270 
271  // if ( !Getv0() || !Getv0PositiveTrack() || !Getv0NegativeTrack() )
272  // { AliInfo(Form("Not adapted for ESDv0s, return true...")); return false; }
273 
274  Double_t mLck0sp,mLcLpi;
275  okLck0sp=1; okLcLpi=1; okLcLbarpi=1;
276 
277  Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
278  Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
279  Double_t mLPDG = TDatabasePDG::Instance()->GetParticle(3122)->Mass();
280 
281  // k0s + p
282  double mk0s = Getv0()->MassK0Short();
283  mLck0sp = InvMassLctoK0sP();
284 
285  // lambda + pi
286  double mlambda = Getv0()->MassLambda();
287  double malambda = Getv0()->MassAntiLambda();
288  mLcLpi = InvMassLctoLambdaPi();
289 
290  // cut on Lc mass
291  // with k0s p hypothesis
292  if(TMath::Abs(mLck0sp-mLcPDG)>cutsLctoV0[0]) okLck0sp = 0;
293  // with Lambda pi hypothesis
294  if(TMath::Abs(mLcLpi-mLcPDG)>cutsLctoV0[1]) okLcLpi = 0;
295  okLcLbarpi = okLcLpi;
296 
297  // cuts on the v0 mass
298  if( TMath::Abs(mk0s-mk0sPDG)>cutsLctoV0[2]) okLck0sp = 0;
299  //if( TMath::Abs(mlambda-mLPDG)>cutsLctoV0[3] &&
300  //TMath::Abs(malambda-mLPDG)>cutsLctoV0[3] ) okLcLpi = 0;
301  if( !(GetBachelor()->Charge()==+1 && TMath::Abs(mlambda-mLPDG)<=cutsLctoV0[3]) ) okLcLpi = 0;
302  if( !(GetBachelor()->Charge()==-1 && TMath::Abs(malambda-mLPDG)<=cutsLctoV0[3]) ) okLcLbarpi = 0;
303 
304  if(!okLck0sp && !okLcLpi && !okLcLbarpi) return 0;
305 
306  // cuts on the minimum pt of the tracks
307  if(TMath::Abs(GetBachelor()->Pt()) < cutsLctoV0[4]) return 0;
308  if(TMath::Abs(Getv0PositiveTrack()->Pt()) < cutsLctoV0[5]) return 0;
309  if(TMath::Abs(Getv0NegativeTrack()->Pt()) < cutsLctoV0[6]) return 0;
310 
311  // cut on the cascade dca
312  if( TMath::Abs(GetDCA(0))>cutsLctoV0[7] //||
313  //TMath::Abs(Getv0()->DcaPosToPrimVertex())>cutsLctoV0[7] ||
314  //TMath::Abs(Getv0()->DcaNegToPrimVertex())>cutsLctoV0[7]
315  ) return 0;
316 
317  // cut on the v0 dca
318  if(TMath::Abs(Getv0()->DcaV0Daughters()) > cutsLctoV0[8]) return 0;
319 
320  // cut on V0 cosine of pointing angle wrt PV
321  if (CosV0PointingAngle() < cutsLctoV0[9]) { // cosine of V0 pointing angle wrt primary vertex
322  AliDebug(4,Form(" V0 cosine of pointing angle doesn't pass the cut"));
323  return 0;
324  }
325 
326  // cut on bachelor transverse impact parameter wrt PV
327  if (TMath::Abs(Getd0Prong(0)) > cutsLctoV0[10]) { // bachelor transverse impact parameter wrt PV
328  AliDebug(4,Form(" bachelor transverse impact parameter doesn't pass the cut"));
329  return 0;
330  }
331 
332  // cut on V0 transverse impact parameter wrt PV
333  if (TMath::Abs(Getd0Prong(1)) > cutsLctoV0[11]) { // V0 transverse impact parameter wrt PV
334  AliDebug(4,Form(" V0 transverse impact parameter doesn't pass the cut"));
335  return 0;
336  }
337 
338  // cut on K0S invariant mass veto
339  if (TMath::Abs(Getv0()->MassK0Short()-mk0sPDG) < cutsLctoV0[12]) { // K0S invariant mass veto
340  AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
341  return 0;
342  }
343 
344  // cut on Lambda/LambdaBar invariant mass veto
345  if (TMath::Abs(Getv0()->MassLambda()-mLPDG) < cutsLctoV0[13] ||
346  TMath::Abs(Getv0()->MassAntiLambda()-mLPDG) < cutsLctoV0[13] ) { // Lambda/LambdaBar invariant mass veto
347  AliDebug(4,Form(" veto on K0S invariant mass doesn't pass the cut"));
348  return 0;
349  }
350 
351  // cut on gamma invariant mass veto
352  if (Getv0()->InvMass2Prongs(0,1,11,11) < cutsLctoV0[14]) { // K0S invariant mass veto
353  AliDebug(4,Form(" veto on gamma invariant mass doesn't pass the cut"));
354  return 0;
355  }
356 
357  // cut on V0 pT min
358  if (Getv0()->Pt() < cutsLctoV0[15]) { // V0 pT min
359  AliDebug(4,Form(" V0 track Pt=%2.2e > %2.2e",Getv0()->Pt(),cutsLctoV0[15]));
360  return 0;
361  }
362 
363  return true;
364 
365 }
366 //-----------------------------------------------------------------------------
371 
372  TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
373  TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
374  TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
375 
376  TVector3 perp = p3Trk0.Cross(p3Trk1);
377  Double_t theta = p3Trk2.Angle(perp);
378  if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
379  theta = TMath::Pi()/2. - theta;
380 
381  return theta;
382 }
383 //-----------------------------------------------------------------------------
388  TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
389  TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
390  TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
391 
392  Double_t alpha = p3Trk0.Angle(p3Trk2);
393  Double_t beta = p3Trk1.Angle(p3Trk2);
394 
395  Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(AngleD0dkpPisoft());
396  Double_t cosphi02 = TMath::Cos(beta) / TMath::Cos(AngleD0dkpPisoft());
397 
398  Double_t phi01 = TMath::ACos(cosphi01);
399  Double_t phi02 = TMath::ACos(cosphi02);
400  Double_t phi00 = p3Trk0.Angle(p3Trk1);
401 
402  if((phi01>phi00) || (phi02>phi00)) return kFALSE;
403  return kTRUE;
404 }
405 
406 //-----------------------------------------------------------------------------
408 {
412 
413  AliAODv0 *v0 = (AliAODv0*)Getv0();
414 
415  if (!v0)
416  return -1.;
417  AliAODVertex *vtxPrimary = GetPrimaryVtx();
418  Double_t posVtx[3] = {0.,0.,0.};
419  vtxPrimary->GetXYZ(posVtx);
420  return v0->DecayLengthV0(posVtx);
421 
422 }
423 //-----------------------------------------------------------------------------
425 {
429  AliAODv0 *v0 = (AliAODv0*)Getv0();
430 
431  if (!v0)
432  return -1.;
433  AliAODVertex *vtxPrimary = GetPrimaryVtx();
434  Double_t posVtx[3] = {0.,0.,0.};
435  vtxPrimary->GetXYZ(posVtx);
436  return v0->DecayLengthXY(posVtx);
437 
438 }
439 //-----------------------------------------------------------------------------
441 {
445 
446  AliAODv0 *v0 = (AliAODv0*)Getv0();
447 
448  if (!v0)
449  return -999.;
450 
451  AliAODVertex *vtxPrimary = GetPrimaryVtx();
452  Double_t posVtx[3] = {0.,0.,0.};
453  vtxPrimary->GetXYZ(posVtx);
454  return v0->CosPointingAngle(posVtx);
455 
456 }
457 //-----------------------------------------------------------------------------
459 {
463 
464  AliAODv0 *v0 = (AliAODv0*)Getv0();
465 
466  if (!v0)
467  return -999.;
468 
469  AliAODVertex *vtxPrimary = GetPrimaryVtx();
470  Double_t posVtx[3] = {0.,0.,0.};
471  vtxPrimary->GetXYZ(posVtx);
472  return v0->CosPointingAngleXY(posVtx);
473 
474 }
475 //-----------------------------------------------------------------------------
477 {
481 
482  AliAODv0 *v0 = (AliAODv0*)Getv0();
483 
484  if (!v0)
485  return -1.;
486  //AliAODVertex *vtxPrimary = GetPrimaryVtx();
487  //Double_t posVtx[3] = {0.,0.,0.};
488  //vtxPrimary->GetXYZ(posVtx);
489  //return v0->NormalizedDecayLength(posVtx);
490  return v0->NormalizedDecayLength(GetPrimaryVtx());
491 
492 }
493 //-----------------------------------------------------------------------------
495 {
499  AliAODv0 *v0 = (AliAODv0*)Getv0();
500 
501  if (!v0)
502  return -1.;
503  //AliAODVertex *vtxPrimary = GetPrimaryVtx();
504  //Double_t posVtx[3] = {0.,0.,0.};
505  //vtxPrimary->GetXYZ(posVtx);
506  //return v0->NormalizedDecayLengthXY(posVtx);
507  return v0->NormalizedDecayLengthXY(GetPrimaryVtx());
508 
509 }
510 //-----------------------------------------------------------------------------
516 
520 
521  if (!okCascLc && !okCascDplus && !okCascDs) {
522  // Cascade candidates don't have any flag: only Lc candidates are stored in the studied delta-AOD
523  //AliDebug(2, "Cascade candidate does not have any flag - candidate accepted if Lc hypothesis is required");
524  if (selFlag==AliRDHFCuts::kLctoV0Cuts) return kTRUE;
525  else return kFALSE;
526  }
527 
528  if (selFlag==AliRDHFCuts::kLctoV0Cuts && okCascLc) return kTRUE;
529  if (selFlag==AliRDHFCuts::kDplustoK0sCuts && okCascDplus) return kTRUE;
530  if (selFlag==AliRDHFCuts::kDstoK0sCuts && okCascDs) return kTRUE;
531 
532  return kFALSE;
533 }
534 
Int_t charge
#define P(T, U, S)
double Double_t
Definition: External.C:58
Double_t DeltaInvMass() const
Bool_t TrigonometricalCut() const
Bool_t HasSelectionBit(Int_t i) const
Double_t NormalizedV0DecayLength() const
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliAODTrack * Getv0NegativeTrack() const
AliAODv0 * Getv0() const
Double_t InvMassLctoLambdaPi() const
Double_t CosV0PointingAngleXY() const
Double_t NormalizedV0DecayLengthXY() const
Double_t DecayLengthXYV0() const
Bool_t SelectDstar(const Double_t *cutsDstar, const Double_t *cutsD0, Bool_t testD0=kTRUE) const
Double_t AngleD0dkpPisoft() const
Double_t InvMassLctoK0sP() const
AliAODTrack * Getv0PositiveTrack() const
int Int_t
Definition: External.C:63
AliAODTrack * GetBachelor() const
Bool_t SelectD0(const Double_t *cuts, Int_t &okD0, Int_t &okD0bar) const
Double_t CosV0PointingAngle() const
AliAODRecoDecayHF2Prong & operator=(const AliAODRecoDecayHF2Prong &source)
short Short_t
Definition: External.C:23
Bool_t SelectLctoV0(const Double_t *cutsLctoV0, Bool_t okLck0sp, Bool_t okLcLpi, Bool_t okLcLbarpi) const
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Bool_t CheckCascadeFlags(AliRDHFCuts::ESele selFlag=AliRDHFCuts::kLctoV0Cuts)
AliAODVertex * GetPrimaryVtx() const
Double_t InvMassDstarKpipi() const
bool Bool_t
Definition: External.C:53
Double_t DecayLengthV0() const
AliAODRecoDecayHF2Prong * Get2Prong() const
AliAODRecoCascadeHF & operator=(const AliAODRecoCascadeHF &source)