]>
Commit | Line | Data |
---|---|---|
1 | /************************************************************************** | |
2 | * Copyright(c) 1998-2006, 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 | // Base class for AOD reconstructed decay | |
19 | // | |
20 | // Author: A.Dainese, andrea.dainese@lnl.infn.it | |
21 | ///////////////////////////////////////////////////////////// | |
22 | ||
23 | #include <TDatabasePDG.h> | |
24 | #include <TVector3.h> | |
25 | #include <TClonesArray.h> | |
26 | #include "AliLog.h" | |
27 | #include "AliVTrack.h" | |
28 | #include "AliAODMCParticle.h" | |
29 | #include "AliAODRecoDecay.h" | |
30 | ||
31 | ClassImp(AliAODRecoDecay) | |
32 | ||
33 | //-------------------------------------------------------------------------- | |
34 | AliAODRecoDecay::AliAODRecoDecay() : | |
35 | AliVTrack(), | |
36 | fSecondaryVtx(0x0), | |
37 | fOwnSecondaryVtx(0x0), | |
38 | fCharge(0), | |
39 | fNProngs(0), fNDCA(0), fNPID(0), | |
40 | fPx(0x0), fPy(0x0), fPz(0x0), | |
41 | fd0(0x0), | |
42 | fDCA(0x0), | |
43 | fPID(0x0), | |
44 | fEventNumber(-1),fRunNumber(-1) | |
45 | { | |
46 | // | |
47 | // Default Constructor | |
48 | // | |
49 | } | |
50 | //-------------------------------------------------------------------------- | |
51 | AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs, | |
52 | Short_t charge, | |
53 | Double_t *px,Double_t *py,Double_t *pz, | |
54 | Double_t *d0) : | |
55 | AliVTrack(), | |
56 | fSecondaryVtx(vtx2), | |
57 | fOwnSecondaryVtx(0x0), | |
58 | fCharge(charge), | |
59 | fNProngs(nprongs), fNDCA(0), fNPID(0), | |
60 | fPx(0x0), fPy(0x0), fPz(0x0), | |
61 | fd0(0x0), | |
62 | fDCA(0x0), | |
63 | fPID(0x0), | |
64 | fEventNumber(-1),fRunNumber(-1) | |
65 | { | |
66 | // | |
67 | // Constructor with AliAODVertex for decay vertex | |
68 | // | |
69 | ||
70 | fPx = new Double_t[GetNProngs()]; | |
71 | fPy = new Double_t[GetNProngs()]; | |
72 | fPz = new Double_t[GetNProngs()]; | |
73 | fd0 = new Double_t[GetNProngs()]; | |
74 | for(Int_t i=0; i<GetNProngs(); i++) { | |
75 | fPx[i] = px[i]; | |
76 | fPy[i] = py[i]; | |
77 | fPz[i] = pz[i]; | |
78 | fd0[i] = d0[i]; | |
79 | } | |
80 | } | |
81 | //-------------------------------------------------------------------------- | |
82 | AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs, | |
83 | Short_t charge, | |
84 | Double_t *d0) : | |
85 | AliVTrack(), | |
86 | fSecondaryVtx(vtx2), | |
87 | fOwnSecondaryVtx(0x0), | |
88 | fCharge(charge), | |
89 | fNProngs(nprongs), fNDCA(0), fNPID(0), | |
90 | fPx(0x0), fPy(0x0), fPz(0x0), | |
91 | fd0(0x0), | |
92 | fDCA(0x0), | |
93 | fPID(0x0), | |
94 | fEventNumber(-1),fRunNumber(-1) | |
95 | { | |
96 | // | |
97 | // Constructor with AliAODVertex for decay vertex and without prongs momenta | |
98 | // | |
99 | ||
100 | fd0 = new Double_t[GetNProngs()]; | |
101 | for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i]; | |
102 | } | |
103 | //-------------------------------------------------------------------------- | |
104 | AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) : | |
105 | AliVTrack(source), | |
106 | fSecondaryVtx(source.fSecondaryVtx), | |
107 | fOwnSecondaryVtx(source.fOwnSecondaryVtx), | |
108 | fCharge(source.fCharge), | |
109 | fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID), | |
110 | fPx(0x0), fPy(0x0), fPz(0x0), | |
111 | fd0(0x0), | |
112 | fDCA(0x0), | |
113 | fPID(0x0), | |
114 | fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber) | |
115 | { | |
116 | // | |
117 | // Copy constructor | |
118 | // | |
119 | if(source.GetNProngs()>0) { | |
120 | fd0 = new Double32_t[GetNProngs()]; | |
121 | memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t)); | |
122 | if(source.fPx) { | |
123 | fPx = new Double32_t[GetNProngs()]; | |
124 | fPy = new Double32_t[GetNProngs()]; | |
125 | fPz = new Double32_t[GetNProngs()]; | |
126 | memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t)); | |
127 | memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t)); | |
128 | memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t)); | |
129 | } | |
130 | if(source.fPID) { | |
131 | fPID = new Double32_t[fNPID]; | |
132 | memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t)); | |
133 | } | |
134 | if(source.fDCA) { | |
135 | fDCA = new Double32_t[fNDCA]; | |
136 | memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t)); | |
137 | } | |
138 | } | |
139 | } | |
140 | //-------------------------------------------------------------------------- | |
141 | AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source) | |
142 | { | |
143 | // | |
144 | // assignment operator | |
145 | // | |
146 | if(&source == this) return *this; | |
147 | fSecondaryVtx = source.fSecondaryVtx; | |
148 | fOwnSecondaryVtx = source.fOwnSecondaryVtx; | |
149 | fCharge = source.fCharge; | |
150 | fNProngs = source.fNProngs; | |
151 | fNDCA = source.fNDCA; | |
152 | fNPID = source.fNPID; | |
153 | fEventNumber = source.fEventNumber; | |
154 | fRunNumber = source.fRunNumber; | |
155 | if(source.GetNProngs()>0) { | |
156 | if(fd0)delete [] fd0; | |
157 | fd0 = new Double32_t[GetNProngs()]; | |
158 | memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t)); | |
159 | if(source.fPx) { | |
160 | if(fPx) delete [] fPx; | |
161 | fPx = new Double32_t[GetNProngs()]; | |
162 | if(fPy) delete [] fPy; | |
163 | fPy = new Double32_t[GetNProngs()]; | |
164 | if(fPz) delete [] fPz; | |
165 | fPz = new Double32_t[GetNProngs()]; | |
166 | memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t)); | |
167 | memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t)); | |
168 | memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t)); | |
169 | } | |
170 | if(source.fPID) { | |
171 | if(fPID) delete [] fPID; | |
172 | fPID = new Double32_t[fNPID]; | |
173 | memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t)); | |
174 | } | |
175 | if(source.fDCA) { | |
176 | if(fDCA) delete [] fDCA; | |
177 | fDCA = new Double32_t[fNDCA]; | |
178 | memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t)); | |
179 | } | |
180 | } | |
181 | return *this; | |
182 | } | |
183 | //-------------------------------------------------------------------------- | |
184 | AliAODRecoDecay::~AliAODRecoDecay() { | |
185 | // | |
186 | // Default Destructor | |
187 | // | |
188 | if(fPx) { delete [] fPx; fPx=NULL; } | |
189 | if(fPy) { delete [] fPy; fPy=NULL; } | |
190 | if(fPz) { delete [] fPz; fPz=NULL; } | |
191 | if(fd0) { delete [] fd0; fd0=NULL; } | |
192 | if(fPID) { delete [] fPID; fPID=NULL; } | |
193 | if(fDCA) { delete [] fDCA; fDCA=NULL; } | |
194 | if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; } | |
195 | } | |
196 | //---------------------------------------------------------------------------- | |
197 | Double_t AliAODRecoDecay::Alpha() const | |
198 | { | |
199 | // | |
200 | // Armenteros-Podolanski alpha for 2-prong decays | |
201 | // | |
202 | if(GetNProngs()!=2) { | |
203 | printf("Can be called only for 2-prong decays"); | |
204 | return (Double_t)-99999.; | |
205 | } | |
206 | return 1.-2./(1.+QlProng(0)/QlProng(1)); | |
207 | } | |
208 | //---------------------------------------------------------------------------- | |
209 | Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const | |
210 | { | |
211 | // | |
212 | // Decay length assuming it is produced at "point" [cm] | |
213 | // | |
214 | return TMath::Sqrt((point[0]-GetSecVtxX()) | |
215 | *(point[0]-GetSecVtxX()) | |
216 | +(point[1]-GetSecVtxY()) | |
217 | *(point[1]-GetSecVtxY()) | |
218 | +(point[2]-GetSecVtxZ()) | |
219 | *(point[2]-GetSecVtxZ())); | |
220 | } | |
221 | //---------------------------------------------------------------------------- | |
222 | Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const | |
223 | { | |
224 | // | |
225 | // Decay length in XY assuming it is produced at "point" [cm] | |
226 | // | |
227 | return TMath::Sqrt((point[0]-GetSecVtxX()) | |
228 | *(point[0]-GetSecVtxX()) | |
229 | +(point[1]-GetSecVtxY()) | |
230 | *(point[1]-GetSecVtxY())); | |
231 | } | |
232 | //---------------------------------------------------------------------------- | |
233 | Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const | |
234 | { | |
235 | // | |
236 | // Cosine of pointing angle in space assuming it is produced at "point" | |
237 | // | |
238 | TVector3 mom(Px(),Py(),Pz()); | |
239 | TVector3 fline(GetSecVtxX()-point[0], | |
240 | GetSecVtxY()-point[1], | |
241 | GetSecVtxZ()-point[2]); | |
242 | ||
243 | Double_t pta = mom.Angle(fline); | |
244 | ||
245 | return TMath::Cos(pta); | |
246 | } | |
247 | //---------------------------------------------------------------------------- | |
248 | Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const | |
249 | { | |
250 | // | |
251 | // Cosine of pointing angle in transverse plane assuming it is produced | |
252 | // at "point" | |
253 | // | |
254 | TVector3 momXY(Px(),Py(),0.); | |
255 | TVector3 flineXY(GetSecVtxX()-point[0], | |
256 | GetSecVtxY()-point[1], | |
257 | 0.); | |
258 | ||
259 | Double_t ptaXY = momXY.Angle(flineXY); | |
260 | ||
261 | return TMath::Cos(ptaXY); | |
262 | } | |
263 | //---------------------------------------------------------------------------- | |
264 | Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const | |
265 | { | |
266 | // | |
267 | // Only for 2-prong decays: | |
268 | // Cosine of decay angle (theta*) in the rest frame of the mother particle | |
269 | // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle, | |
270 | // pdgprong0 for prong 0 and pdgprong1 for prong1 | |
271 | // | |
272 | if(GetNProngs()!=2) { | |
273 | printf("Can be called only for 2-prong decays"); | |
274 | return (Double_t)-99999.; | |
275 | } | |
276 | Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass(); | |
277 | Double_t massp[2]; | |
278 | massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass(); | |
279 | massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass(); | |
280 | ||
281 | Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx); | |
282 | ||
283 | Double_t beta = P()/E(pdgvtx); | |
284 | Double_t gamma = E(pdgvtx)/massvtx; | |
285 | ||
286 | Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar; | |
287 | ||
288 | return cts; | |
289 | } | |
290 | //--------------------------------------------------------------------------- | |
291 | Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const | |
292 | { | |
293 | // | |
294 | // Decay time * c assuming it is produced at "point" [cm] | |
295 | // | |
296 | Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
297 | return DecayLength(point)*mass/P(); | |
298 | } | |
299 | //--------------------------------------------------------------------------- | |
300 | Double_t AliAODRecoDecay::E(UInt_t pdg) const | |
301 | { | |
302 | // | |
303 | // Energy | |
304 | // | |
305 | Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
306 | return TMath::Sqrt(mass*mass+P()*P()); | |
307 | } | |
308 | //--------------------------------------------------------------------------- | |
309 | Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const | |
310 | { | |
311 | // | |
312 | // Energy of ip-th prong | |
313 | // | |
314 | Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
315 | return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip)); | |
316 | } | |
317 | //-------------------------------------------------------------------------- | |
318 | Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const { | |
319 | // | |
320 | // This function returns the global covariance matrix of the track params | |
321 | // | |
322 | // Cov(x,x) ... : cv[0] | |
323 | // Cov(y,x) ... : cv[1] cv[2] | |
324 | // Cov(z,x) ... : cv[3] cv[4] cv[5] | |
325 | // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9] | |
326 | // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14] | |
327 | // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20] | |
328 | // | |
329 | // For XYZ we take the cov of the vertex, for PxPyPz we take the | |
330 | // sum of the covs of PxPyPz from the daughters, for the moment | |
331 | // we set the cov between position and momentum as the sum of | |
332 | // the same cov from the daughters. | |
333 | // | |
334 | ||
335 | Int_t j; | |
336 | for(j=0;j<21;j++) cv[j]=0.; | |
337 | ||
338 | if(!GetNDaughters()) { | |
339 | AliError("No daughters available"); | |
340 | return kFALSE; | |
341 | } | |
342 | ||
343 | Double_t v[6]; | |
344 | AliAODVertex *secv=GetSecondaryVtx(); | |
345 | if(!secv) { | |
346 | AliError("Vertex covariance matrix not available"); | |
347 | return kFALSE; | |
348 | } | |
349 | if(!secv->GetCovMatrix(v)) { | |
350 | AliError("Vertex covariance matrix not available"); | |
351 | return kFALSE; | |
352 | } | |
353 | ||
354 | Double_t p[21]; for(j=0;j<21;j++) p[j]=0.; | |
355 | Bool_t error=kFALSE; | |
356 | for(Int_t i=1; i<GetNDaughters(); i++) { | |
357 | AliVTrack *daugh = (AliVTrack*)GetDaughter(i); | |
358 | Double_t dcov[21]; | |
359 | if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE; | |
360 | for(j=0;j<21;j++) p[j] += dcov[j]; | |
361 | } | |
362 | if(error) { | |
363 | AliError("No covariance for at least one daughter") | |
364 | return kFALSE; | |
365 | } | |
366 | ||
367 | for(j=0; j<21; j++) { | |
368 | if(j<6) { | |
369 | cv[j] = v[j]; | |
370 | } else { | |
371 | cv[j] = p[j]; | |
372 | } | |
373 | } | |
374 | ||
375 | return kTRUE; | |
376 | } | |
377 | //---------------------------------------------------------------------------- | |
378 | UChar_t AliAODRecoDecay::GetITSClusterMap() const { | |
379 | // | |
380 | // We take the logical AND of the daughters cluster maps | |
381 | // (only if all daughters have the bit for given layer, we set the bit) | |
382 | // | |
383 | UChar_t map=0; | |
384 | ||
385 | if(!GetNDaughters()) { | |
386 | AliError("No daughters available"); | |
387 | return map; | |
388 | } | |
389 | ||
390 | for(Int_t l=0; l<12; l++) { // loop on ITS layers (from here we cannot know how many they are; let's put 12 to be conservative) | |
391 | Int_t bit = 1; | |
392 | for(Int_t i=0; i<GetNDaughters(); i++) { | |
393 | AliVTrack *daugh = (AliVTrack*)GetDaughter(i); | |
394 | if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0; | |
395 | } | |
396 | if(bit) SETBIT(map,l); | |
397 | } | |
398 | ||
399 | return map; | |
400 | } | |
401 | //-------------------------------------------------------------------------- | |
402 | ULong_t AliAODRecoDecay::GetStatus() const { | |
403 | // | |
404 | // Same as for ITSClusterMap | |
405 | // | |
406 | ULong_t status=0; | |
407 | ||
408 | if(!GetNDaughters()) { | |
409 | AliError("No daughters available"); | |
410 | return status; | |
411 | } | |
412 | ||
413 | AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0); | |
414 | status = status&(daugh0->GetStatus()); | |
415 | ||
416 | for(Int_t i=1; i<GetNDaughters(); i++) { | |
417 | AliVTrack *daugh = (AliVTrack*)GetDaughter(i); | |
418 | status = status&(daugh->GetStatus()); | |
419 | } | |
420 | ||
421 | return status; | |
422 | } | |
423 | //-------------------------------------------------------------------------- | |
424 | Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const | |
425 | { | |
426 | // | |
427 | // Impact parameter in the bending plane of the particle | |
428 | // w.r.t. to "point" | |
429 | // | |
430 | Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py(); | |
431 | k /= Pt()*Pt(); | |
432 | Double_t dx = GetSecVtxX()-point[0]+k*Px(); | |
433 | Double_t dy = GetSecVtxY()-point[1]+k*Py(); | |
434 | Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy); | |
435 | TVector3 mom(Px(),Py(),Pz()); | |
436 | TVector3 fline(GetSecVtxX()-point[0], | |
437 | GetSecVtxY()-point[1], | |
438 | GetSecVtxZ()-point[2]); | |
439 | TVector3 cross = mom.Cross(fline); | |
440 | return (cross.Z()>0. ? absImpPar : -absImpPar); | |
441 | } | |
442 | //---------------------------------------------------------------------------- | |
443 | Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const | |
444 | { | |
445 | // | |
446 | // Invariant mass for prongs mass hypotheses in pdg array | |
447 | // | |
448 | if(GetNProngs()!=npdg) { | |
449 | printf("npdg != GetNProngs()"); | |
450 | return (Double_t)-99999.; | |
451 | } | |
452 | Double_t energysum = 0.; | |
453 | ||
454 | for(Int_t i=0; i<GetNProngs(); i++) { | |
455 | energysum += EProng(i,pdg[i]); | |
456 | } | |
457 | ||
458 | Double_t mass = TMath::Sqrt(energysum*energysum-P()*P()); | |
459 | ||
460 | return mass; | |
461 | } | |
462 | //---------------------------------------------------------------------------- | |
463 | Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2, | |
464 | UInt_t pdg1,UInt_t pdg2) const | |
465 | { | |
466 | // | |
467 | // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2 | |
468 | // | |
469 | Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2); | |
470 | Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2)) | |
471 | +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2)) | |
472 | +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2)); | |
473 | Double_t mass = TMath::Sqrt(energysum*energysum-psum2); | |
474 | ||
475 | return mass; | |
476 | } | |
477 | //---------------------------------------------------------------------------- | |
478 | Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray) const | |
479 | { | |
480 | // | |
481 | // Check if this candidate is matched to a MC signal | |
482 | // If no, return -1 | |
483 | // If yes, return label (>=0) of the AliAODMCParticle | |
484 | // | |
485 | ||
486 | if(!GetNDaughters()) { | |
487 | AliError("No daughters available"); | |
488 | return -1; | |
489 | } | |
490 | ||
491 | Int_t *dgLabels = new Int_t[GetNDaughters()]; | |
492 | ||
493 | // loop on daughters and write the labels | |
494 | for(Int_t i=0; i<GetNDaughters(); i++) { | |
495 | AliVTrack *trk = (AliVTrack*)GetDaughter(i); | |
496 | dgLabels[i] = trk->GetLabel(); | |
497 | } | |
498 | ||
499 | Int_t labMother = MatchToMC(pdgabs,mcArray,dgLabels); | |
500 | ||
501 | delete [] dgLabels; dgLabels=NULL; | |
502 | ||
503 | return labMother; | |
504 | } | |
505 | //---------------------------------------------------------------------------- | |
506 | Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray, | |
507 | Int_t *dgLabels) const | |
508 | { | |
509 | // | |
510 | // Check if this candidate is matched to a MC signal | |
511 | // If no, return -1 | |
512 | // If yes, return label (>=0) of the AliAODMCParticle | |
513 | // | |
514 | ||
515 | if(!GetNDaughters()) { | |
516 | AliError("No daughters available"); | |
517 | return -1; | |
518 | } | |
519 | ||
520 | Int_t *labMom = new Int_t[GetNDaughters()]; | |
521 | Int_t i,lab,labMother,pdgMother; | |
522 | AliAODMCParticle *part=0; | |
523 | Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.; | |
524 | ||
525 | // loop on daughter labels | |
526 | for(i=0; i<GetNDaughters(); i++) { | |
527 | labMom[i]=-1; | |
528 | lab = dgLabels[i]; | |
529 | if(lab<0) { | |
530 | printf("daughter with negative label\n"); | |
531 | continue; | |
532 | } | |
533 | part = (AliAODMCParticle*)mcArray->At(lab); | |
534 | if(!part) { | |
535 | printf("no MC particle\n"); | |
536 | continue; | |
537 | } | |
538 | // keep sum of daughters' momenta, to check for mom conservation | |
539 | pxSumDgs += part->Px(); | |
540 | pySumDgs += part->Py(); | |
541 | pzSumDgs += part->Pz(); | |
542 | // for the J/psi, check that the daughters are electrons | |
543 | if(pdgabs==443 && TMath::Abs(part->GetPdgCode())!=11) continue; | |
544 | ||
545 | while(part->GetMother()>=0) { | |
546 | labMother=part->GetMother(); | |
547 | part = (AliAODMCParticle*)mcArray->At(labMother); | |
548 | if(!part) { | |
549 | printf("no MC mother particle\n"); | |
550 | break; | |
551 | } | |
552 | pdgMother = TMath::Abs(part->GetPdgCode()); | |
553 | if(pdgMother==pdgabs) { | |
554 | labMom[i]=labMother; | |
555 | break; | |
556 | } | |
557 | } | |
558 | } | |
559 | ||
560 | // check if the candidate is signal | |
561 | Bool_t isSignal=kTRUE; | |
562 | labMother=labMom[0]; | |
563 | // all labels have to be the same and !=-1 | |
564 | for(i=0; i<GetNDaughters(); i++) { | |
565 | if(labMom[i]==-1 || labMom[i]!=labMother) isSignal=kFALSE; | |
566 | } | |
567 | ||
568 | delete [] labMom; labMom=NULL; | |
569 | ||
570 | if(!isSignal) return -1; | |
571 | ||
572 | /* | |
573 | // check that the mother decayed in <GetNDaughters()> prongs | |
574 | Int_t ndg = TMath::Abs(part->GetDaughter(1)-part->GetDaughter(0))+1; | |
575 | if(ndg!=GetNDaughters()) return -1; | |
576 | AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(part->GetDaughter(1))); | |
577 | AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(part->GetDaughter(0))); | |
578 | printf("pdg %d %d %d %d %d %d\n",part->GetDaughter(1),part->GetDaughter(0),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode()); | |
579 | */ | |
580 | ||
581 | // the above works only for non-resonant decays, | |
582 | // it's better to check for mom conservation | |
583 | part = (AliAODMCParticle*)mcArray->At(labMother); | |
584 | Double_t pxMother = part->Px(); | |
585 | Double_t pyMother = part->Py(); | |
586 | Double_t pzMother = part->Pz(); | |
587 | // within 0.1% | |
588 | if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.001 && | |
589 | (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.001 && | |
590 | (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.001) | |
591 | return -1; | |
592 | ||
593 | return labMother; | |
594 | } | |
595 | //--------------------------------------------------------------------------- | |
596 | void AliAODRecoDecay::Print(Option_t* /*option*/) const | |
597 | { | |
598 | // | |
599 | // Print some information | |
600 | // | |
601 | printf("AliAODRecoDecay with %d prongs\n",GetNProngs()); | |
602 | printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ()); | |
603 | ||
604 | return; | |
605 | } | |
606 | //---------------------------------------------------------------------------- | |
607 | Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const | |
608 | { | |
609 | // | |
610 | // Relative angle between two prongs | |
611 | // | |
612 | TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1)); | |
613 | TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2)); | |
614 | ||
615 | Double_t angle = momA.Angle(momB); | |
616 | ||
617 | return angle; | |
618 | } | |
619 | //---------------------------------------------------------------------------- | |
620 | Double_t AliAODRecoDecay::QlProng(Int_t ip) const | |
621 | { | |
622 | // | |
623 | // Longitudinal momentum of prong w.r.t. to total momentum | |
624 | // | |
625 | TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip)); | |
626 | TVector3 momTot(Px(),Py(),Pz()); | |
627 | ||
628 | return mom.Dot(momTot)/momTot.Mag(); | |
629 | } | |
630 | //---------------------------------------------------------------------------- | |
631 | Double_t AliAODRecoDecay::QtProng(Int_t ip) const | |
632 | { | |
633 | // | |
634 | // Transverse momentum of prong w.r.t. to total momentum | |
635 | // | |
636 | TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip)); | |
637 | TVector3 momTot(Px(),Py(),Pz()); | |
638 | ||
639 | return mom.Perp(momTot); | |
640 | } | |
641 | //---------------------------------------------------------------------------- | |
642 | Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const | |
643 | { | |
644 | // | |
645 | // Longitudinal momentum of prong w.r.t. to flight line between "point" | |
646 | // and fSecondaryVtx | |
647 | // | |
648 | TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip)); | |
649 | TVector3 fline(GetSecVtxX()-point[0], | |
650 | GetSecVtxY()-point[1], | |
651 | GetSecVtxZ()-point[2]); | |
652 | ||
653 | return mom.Dot(fline)/fline.Mag(); | |
654 | } | |
655 | //---------------------------------------------------------------------------- | |
656 | Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const | |
657 | { | |
658 | // | |
659 | // Transverse momentum of prong w.r.t. to flight line between "point" and | |
660 | // fSecondaryVtx | |
661 | // | |
662 | TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip)); | |
663 | TVector3 fline(GetSecVtxX()-point[0], | |
664 | GetSecVtxY()-point[1], | |
665 | GetSecVtxZ()-point[2]); | |
666 | ||
667 | return mom.Perp(fline); | |
668 | } | |
669 | //-------------------------------------------------------------------------- |