Updates in LC->Kos+proton code (Annalisa)
[u/mrichter/AliRoot.git] / 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 }