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