Use new method AliAODRecoDecay::MatchToMC that simplifies matching to MC
[u/mrichter/AliRoot.git] / 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
97   Double_t px[3],py[3],pz[3];
98   UInt_t pdg[3]={321,211,211};
99   pdg[0] = (Charge()>0 ? 211 : 321); // positive daughter of D0
100   px[0] = Get2Prong()->PxProng(0);
101   py[0] = Get2Prong()->PyProng(0);
102   pz[0] = Get2Prong()->PzProng(0);
103   pdg[1] = (Charge()>0 ? 321 : 211); // negative daughter of D0
104   px[1] = Get2Prong()->PxProng(1);
105   py[1] = Get2Prong()->PyProng(1);
106   pz[1] = Get2Prong()->PzProng(1);
107   pdg[2] = 211; // soft pion
108   px[2] = PxProng(0);
109   py[2] = PyProng(0);
110   pz[2] = PzProng(0);
111   Short_t dummycharge=0;
112   Double_t dummyd0[3]={0,0,0};
113   AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,3,dummycharge,px,py,pz,dummyd0);
114
115   Double_t minv = rd->InvMass(3,pdg);
116
117   delete rd; rd=NULL;
118
119   return minv;
120 }
121 //----------------------------------------------------------------------------
122 Int_t AliAODRecoCascadeHF::MatchToMC(Int_t pdgabs,Int_t pdgabs2prong,
123                                      TClonesArray *mcArray) const
124 {
125   //
126   // Check if this candidate is matched to a MC signal
127   // If no, return -1
128   // If yes, return label (>=0) of the AliAODMCParticle
129   // 
130
131   if(pdgabs!=413) {
132     printf("Only D*+ implemented\n"); 
133     return -1;
134   }
135
136   if(!GetNDaughters()) {
137     AliError("No daughters available");
138     return -1;
139   }
140   
141   AliAODRecoDecayHF2Prong *the2Prong = Get2Prong();
142   Int_t lab2Prong = the2Prong->MatchToMC(pdgabs2prong,mcArray);
143
144   if(lab2Prong<0) return -1;
145
146   Int_t *dgLabels = new Int_t[GetNDaughters()];
147
148   // loop on daughters and write labels
149   for(Int_t i=0; i<GetNDaughters(); i++) {
150     AliVTrack *trk = (AliVTrack*)GetDaughter(i);
151     Int_t lab = trk->GetLabel();
152     if(lab==-1) { // this daughter is the 2prong
153       lab=lab2Prong;
154     } else if(lab<-1) {
155       printf("daughter with negative label\n");
156       continue;
157     }
158     dgLabels[i] = lab;
159   }
160
161   Int_t labMother = AliAODRecoDecay::MatchToMC(pdgabs,mcArray,dgLabels);
162
163   delete [] dgLabels; dgLabels=NULL;
164
165   return labMother;
166 }
167 //-----------------------------------------------------------------------------
168 Bool_t AliAODRecoCascadeHF::SelectDstar(const Double_t *cutsDstar,
169                                         const Double_t *cutsD0,
170                                         Bool_t testD0) const
171 {
172   //
173   // cutsDstar[0] = inv. mass half width of D* [GeV]
174   // cutsDstar[1] = half width of (M_Kpipi-M_D0) [GeV]
175   // cutsDstar[2] = PtMin of pi_s [GeV/c]
176   // cutsDstar[3] = PtMax of pi_s [GeV/c]
177   // cutsDstar[4] = theta, angle between the pi_s and decay plane of the D0 [rad]
178   //
179   // cutsD0[0] = inv. mass half width [GeV]   
180   // cutsD0[1] = dca [cm]
181   // cutsD0[2] = cosThetaStar 
182   // cutsD0[3] = pTK [GeV/c]
183   // cutsD0[4] = pTPi [GeV/c]
184   // cutsD0[5] = d0K [cm]   upper limit!
185   // cutsD0[6] = d0Pi [cm]  upper limit!
186   // cutsD0[7] = d0d0 [cm^2]
187   // cutsD0[8] = cosThetaPoint
188
189
190   // check that the D0 passes the cuts
191   // (if we have a D*+, it has to pass as D0, 
192   //  if we have a D*-, it has to pass as D0bar)
193
194   if(testD0) {
195     Int_t okD0=0,okD0bar=0;
196     Get2Prong()->SelectD0(cutsD0,okD0,okD0bar);
197     if((Charge()==+1 && !okD0) || (Charge()==-1 && !okD0bar)) return kFALSE; 
198   }
199  
200   if( (PtProng(0)<cutsDstar[2]) || (PtProng(0)>cutsDstar[3]) ) return kFALSE;
201
202   Double_t mDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
203   Double_t invmDstar = InvMassDstarKpipi();
204   if(TMath::Abs(mDstar-invmDstar)>cutsDstar[0]) return kFALSE;
205
206   Double_t mD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
207   if(TMath::Abs((mDstar-mD0)-DeltaInvMass())>cutsDstar[1]) return kFALSE;
208
209   TVector3 p3Trk0(Get2Prong()->PxProng(0),Get2Prong()->PyProng(0),Get2Prong()->PzProng(0)); // from D0
210   TVector3 p3Trk1(Get2Prong()->PxProng(1),Get2Prong()->PyProng(1),Get2Prong()->PzProng(1)); // from D0
211   TVector3 p3Trk2(PxProng(0),PyProng(0),PzProng(0)); // pi_s
212
213   TVector3 perp = p3Trk0.Cross(p3Trk1);
214   Double_t theta = p3Trk2.Angle(perp);
215   if(theta>(TMath::Pi()-theta)) theta = TMath::Pi() - theta;
216   theta = TMath::Pi()/2. - theta;
217
218   if(theta>cutsDstar[4]) return kFALSE;
219
220   Double_t alpha = p3Trk0.Angle(p3Trk2);
221   Double_t belta = p3Trk1.Angle(p3Trk2);
222
223   Double_t cosphi01 = TMath::Cos(alpha) / TMath::Cos(theta);
224   Double_t cosphi02 = TMath::Cos(belta) / TMath::Cos(theta);
225
226   Double_t phi01 = TMath::ACos(cosphi01);
227   Double_t phi02 = TMath::ACos(cosphi02);
228   Double_t phi00 = p3Trk0.Angle(p3Trk1);
229
230   if((phi01>phi00) || (phi02>phi00)) return kFALSE;
231   
232   return kTRUE;
233 }
234 //-----------------------------------------------------------------------------