1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDtrackInfo.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // Reconstruction QA //
23 // Alex Bercuci <A.Bercuci@gsi.de> //
24 // Markus Fasel <M.Fasel@gsi.de> //
26 ////////////////////////////////////////////////////////////////////////////
28 #include "TDatabasePDG.h"
32 #include "AliTrackReference.h"
33 #include "AliExternalTrackParam.h"
35 #include "AliTRDseedV1.h"
36 #include "AliTRDtrackV1.h"
37 #include "AliTRDgeometry.h"
38 #include "AliTRDtrackerV1.h"
40 #include "AliTRDtrackInfo.h"
42 ClassImp(AliTRDtrackInfo)
43 ClassImp(AliTRDtrackInfo::AliMCinfo)
44 ClassImp(AliTRDtrackInfo::AliESDinfo)
45 Double_t AliTRDtrackInfo::AliMCinfo::fgKalmanStep = 2.;
46 Bool_t AliTRDtrackInfo::AliMCinfo::fgKalmanUpdate = kTRUE;
48 //___________________________________________________
49 AliTRDtrackInfo::AliTRDtrackInfo():
57 // Default constructor
62 //___________________________________________________
63 AliTRDtrackInfo::AliTRDtrackInfo(const AliTRDtrackInfo &trdInfo):
64 TObject((const TObject&)trdInfo)
65 ,fNClusters(trdInfo.fNClusters)
74 if(trdInfo.fMC) fMC = new AliMCinfo(*trdInfo.fMC);
76 if(trdInfo.fTRDtrack){
77 fTRDtrack = new AliTRDtrackV1(*trdInfo.fTRDtrack);
78 if(trdInfo.fTRDtrack->IsOwner()) fTRDtrack->SetOwner();
82 //___________________________________________________
83 AliTRDtrackInfo::AliMCinfo::AliMCinfo()
89 memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
92 //___________________________________________________
93 AliTRDtrackInfo::AliMCinfo::AliMCinfo(const AliMCinfo &mc)
96 ,fNTrackRefs(mc.fNTrackRefs)
102 memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
103 for(Int_t ien = 0; ien < 12; ien++){
104 if(mc.fTrackRefs[ien])
105 fTrackRefs[ien] = new AliTrackReference(*(mc.fTrackRefs[ien]));
109 //________________________________________________________
110 Int_t AliTRDtrackInfo::AliMCinfo::GetPID() const
112 // Translate pdg code to PID index
116 case kPositron: return AliPID::kElectron;
118 case kMuonMinus: return AliPID::kMuon;
120 case kPiMinus: return AliPID::kPion;
122 case kKMinus: return AliPID::kKaon;
124 case kProtonBar: return AliPID::kProton;
129 //___________________________________________________
130 AliTRDtrackInfo::AliESDinfo::AliESDinfo()
147 memset(fTRDr, 0, AliPID::kSPECIES*sizeof(Double32_t));
148 memset(fTRDv0pid, 0, AliPID::kSPECIES*sizeof(Int_t));
151 //___________________________________________________
152 AliTRDtrackInfo::AliESDinfo::AliESDinfo(const AliESDinfo &esd)
155 ,fStatus(esd.fStatus)
156 ,fKinkIndex(esd.fKinkIndex)
157 ,fTPCncls(esd.fTPCncls)
159 ,fTRDpidQuality(esd.fTRDpidQuality)
160 ,fTRDnSlices(esd.fTRDnSlices)
169 memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
170 memcpy(fTRDv0pid, esd.fTRDv0pid, AliPID::kSPECIES*sizeof(Int_t));
173 fTRDslices = new Double32_t[fTRDnSlices];
174 memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
176 if(esd.fOP) fOP = new AliExternalTrackParam(*esd.fOP);
177 if(esd.fTPCout) fTPCout = new AliExternalTrackParam(*esd.fTPCout);
181 //___________________________________________________
182 AliTRDtrackInfo::~AliTRDtrackInfo()
189 if(fTRDtrack) delete fTRDtrack;
192 //___________________________________________________
193 AliTRDtrackInfo::AliMCinfo::~AliMCinfo()
200 for(Int_t ien = 0; ien < 12; ien++){
201 if(fTrackRefs[ien]) delete fTrackRefs[ien];
202 fTrackRefs[ien] = NULL;
206 //___________________________________________________
207 AliTRDtrackInfo::AliESDinfo::~AliESDinfo()
214 delete [] fTRDslices;
218 if(fOP) delete fOP; fOP = NULL;
219 if(fTPCout) delete fTPCout; fTPCout = NULL;
222 //___________________________________________________
223 void AliTRDtrackInfo::AliESDinfo::Delete(const Option_t *){
225 // Delete Pointer members
228 delete [] fTRDslices;
232 if(fOP) delete fOP; fOP = NULL;
233 if(fTPCout) delete fTPCout; fTPCout = NULL;
237 //___________________________________________________
238 AliTRDtrackInfo& AliTRDtrackInfo::operator=(const AliTRDtrackInfo &trdInfo)
243 if(this == &trdInfo) return *this;
245 fNClusters = trdInfo.fNClusters;
249 if(!fMC) fMC = new AliMCinfo(*trdInfo.fMC);
252 new(fMC) AliMCinfo(*trdInfo.fMC);
254 } else {if(fMC) delete fMC; fMC = NULL;}
256 if(trdInfo.fTRDtrack){
257 if(!fTRDtrack) fTRDtrack = new AliTRDtrackV1(*trdInfo.fTRDtrack);
259 fTRDtrack->~AliTRDtrackV1();
260 new(fTRDtrack) AliTRDtrackV1(*trdInfo.fTRDtrack);
262 if(trdInfo.fTRDtrack->IsOwner()) fTRDtrack->SetOwner();
263 } else {if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;}
268 //___________________________________________________
269 AliTRDtrackInfo::AliMCinfo& AliTRDtrackInfo::AliMCinfo::operator=(const AliMCinfo &mc)
272 // Assignment operator
275 if(this == &mc) return *this;
278 fNTrackRefs = mc.fNTrackRefs;
280 AliTrackReference **itr = &fTrackRefs[0];
281 AliTrackReference* const *jtr = &mc.fTrackRefs[0];
282 for(Int_t ien = 0; ien < 12; ien++, itr++, jtr++){
284 if(!(*itr)) (*itr) = new AliTrackReference(*(*jtr));
286 (*itr)->~AliTrackReference();
287 new(&(*itr)) AliTrackReference(*(*jtr));
289 } else {if((*itr)) delete (*itr); (*itr) = NULL;}
294 //___________________________________________________
295 AliTRDtrackInfo::AliESDinfo& AliTRDtrackInfo::AliESDinfo::operator=(const AliESDinfo &esd)
298 // Assignment operator
301 if(this == &esd) return *this;
304 fStatus = esd.fStatus;
305 fKinkIndex = esd.fKinkIndex;
306 fTPCncls = esd.fTPCncls;
308 fTRDpidQuality= esd.fTRDpidQuality;
309 fTRDnSlices = esd.fTRDnSlices;
311 memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
312 memcpy(fTRDv0pid, esd.fTRDv0pid, AliPID::kSPECIES*sizeof(Int_t));
315 if(!fTRDslices) fTRDslices = new Double32_t[fTRDnSlices];
316 memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
320 fOP->~AliExternalTrackParam();
321 // RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
322 new(fOP) AliExternalTrackParam(*esd.fOP);
323 } else fOP = new AliExternalTrackParam(*esd.fOP);
324 } else {if(fOP) delete fOP; fOP = NULL;}
327 fTPCout->~AliExternalTrackParam();
328 // RS: Constructor from VTrack was used instead of Constructor from AliExternalTrackParam
329 new(fTPCout) AliExternalTrackParam(*esd.fTPCout);
330 } else fTPCout = new AliExternalTrackParam(*esd.fTPCout);
331 } else { if(fTPCout) delete fTPCout; fTPCout = NULL;}
336 //___________________________________________________
337 void AliTRDtrackInfo::Delete(const Option_t *)
343 AliDebug(2, Form("track[%p] mc[%p]", (void*)fTRDtrack, (void*)fMC));
345 if(fMC) delete fMC; fMC = NULL;
347 if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
350 //___________________________________________________
351 void AliTRDtrackInfo::SetTrack(const AliTRDtrackV1 *track)
357 if(!fTRDtrack) fTRDtrack = new AliTRDtrackV1(*track);
359 fTRDtrack->~AliTRDtrackV1();
360 new(fTRDtrack) AliTRDtrackV1(*track);
362 fTRDtrack->SetOwner();
363 // Make a copy for the object in order to avoid ownership problems
366 //___________________________________________________
367 void AliTRDtrackInfo::AddTrackRef(const AliTrackReference *tref)
370 // Add track reference
373 if(fMC->fNTrackRefs >= 12){
377 // Make a copy for the object in order to avoid ownership problems
378 fMC->fTrackRefs[fMC->fNTrackRefs++] = new AliTrackReference(*tref);
381 //___________________________________________________
382 AliTrackReference* AliTRDtrackInfo::GetTrackRef(Int_t idx) const
385 // Returns a track reference
387 if(!fMC) return NULL;
388 return (idx>=0 && idx < 12) ? fMC->fTrackRefs[idx] : NULL;
391 //___________________________________________________
392 AliTrackReference* AliTRDtrackInfo::GetTrackRef(const AliTRDseedV1* const tracklet) const
395 // Returns a track reference
397 if(!fMC) return NULL;
398 Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
399 AliTrackReference * const* jtr = &(fMC->fTrackRefs[0]);
400 for(Int_t itr = 0; itr < fMC->fNTrackRefs; itr++, ++jtr){
402 if(TMath::Abs(tracklet->GetX0() - (*jtr)->LocalX()) < cw) return (*jtr);
407 //___________________________________________________
408 Int_t AliTRDtrackInfo::GetNumberOfClusters() const
411 // Returns the number of clusters
415 if(!fTRDtrack) return 0;
416 if(fTRDtrack->GetNumberOfTracklets() == 0) return n;
417 AliTRDseedV1 *tracklet = NULL;
418 for(Int_t ip=0; ip<6; ip++){
419 if(!(tracklet = const_cast<AliTRDseedV1 *>(fTRDtrack->GetTracklet(ip)))) continue;
426 //___________________________________________________
427 void AliTRDtrackInfo::SetOuterParam(const AliExternalTrackParam *op)
430 // Set outer track parameters
435 fESD.fOP->~AliExternalTrackParam();
436 new(fESD.fOP) AliExternalTrackParam(*op);
437 } else fESD.fOP = new AliExternalTrackParam(*op);
440 //___________________________________________________
441 void AliTRDtrackInfo::SetTPCoutParam(const AliExternalTrackParam *op)
444 // Set TPCout track parameters
449 fESD.fTPCout->~AliExternalTrackParam();
450 new(fESD.fTPCout) AliExternalTrackParam(*op);
451 } else fESD.fTPCout = new AliExternalTrackParam(*op);
454 //___________________________________________________
455 Int_t AliTRDtrackInfo::GetNTracklets() const
458 // Return the number of tracklets
461 if(!fTRDtrack) return 0;
462 return fTRDtrack->GetNumberOfTracklets();
465 //___________________________________________________
466 void AliTRDtrackInfo::SetSlices(Int_t n, Double32_t *s)
471 if(fESD.fTRDnSlices != n){
472 fESD.fTRDnSlices = 0;
473 delete [] fESD.fTRDslices;
474 fESD.fTRDslices = NULL;
477 if(!fESD.fTRDnSlices){
478 fESD.fTRDnSlices = n;
479 fESD.fTRDslices = new Double32_t[fESD.fTRDnSlices];
482 memcpy(fESD.fTRDslices, s, n*sizeof(Double32_t));
485 //___________________________________________________
486 Bool_t AliTRDtrackInfo::AliMCinfo::GetDirections(Float_t &x0, Float_t &y0, Float_t &z0, Float_t &dydx, Float_t &dzdx, Float_t &pt, Float_t &eta, Float_t &phi, UChar_t &status) const
488 // Check for 2 track ref for the tracklet defined bythe radial position x0
489 // The "status" is a bit map and gives a more informative output in case of failure:
490 // - 0 : everything is OK
491 // - BIT(0) : 0 track refs found
492 // - BIT(1) : 1 track reference found
493 // - BIT(2) : dx <= 0 between track references
494 // - BIT(3) : dx > 0 && dx < 3.7 - tangent tracks
497 Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
499 AliTrackReference *tr[2] = {NULL, NULL};
500 AliTrackReference * const* jtr = &fTrackRefs[0];
501 for(Int_t itr = 0; itr < fNTrackRefs; itr++, ++jtr){
504 if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, (*jtr)->LocalX());*/
505 if(TMath::Abs(x0 - (*jtr)->LocalX()) > cw) continue;
506 tr[nFound++] = (*jtr);
507 if(nFound == 2) break;
510 AliDebugGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", 1, Form("Missing track ref x0[%6.3f] nref[%d]", x0, nFound));
511 if(!nFound) SETBIT(status, 0);
512 else SETBIT(status, 1);
516 if(pt < 1.e-3) return kFALSE;
518 Double_t dx = tr[1]->LocalX() - tr[0]->LocalX();
520 AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Track ref with wrong radial distances refX0[%6.3f] refX1[%6.3f]", tr[0]->LocalX(), tr[1]->LocalX()));
524 if(TMath::Abs(dx-AliTRDgeometry::CamHght()-AliTRDgeometry::CdrHght())>1.E-3) SETBIT(status, 3);
526 dydx = (tr[1]->LocalY() - tr[0]->LocalY()) / dx;
527 dzdx = (tr[1]->Z() - tr[0]->Z()) / dx;
528 //Float_t dx0 = tr[1]->LocalX() - x0;
529 y0 = tr[1]->LocalY()/* - dydx*dx0*/;
530 z0 = tr[1]->Z()/* - dzdx*dx0*/;
531 x0 = tr[1]->LocalX();
532 eta = -TMath::Log(TMath::Tan(0.5 * tr[1]->Theta()));
533 phi = TMath::ATan2(tr[1]->Y(), tr[1]->X());
537 //___________________________________________________
538 Bool_t AliTRDtrackInfo::AliMCinfo::PropagateKalman(
539 TVectorD *x, TVectorD *y, TVectorD *z,
540 TVectorD *dx, TVectorD *dy, TVectorD *dz,
541 TVectorD *pt, TVectorD *dpt, TVectorD *budget, TVectorD *c, Double_t mass) const
543 // Propagate Kalman from the first TRD track reference to
544 // last one and save residuals in the y, z and pt.
546 // This is to calibrate the dEdx and MS corrections
548 if(!fNTrackRefs) return kFALSE;
549 for(Int_t itr=kNTrackRefs; itr--;){
550 (*x)[itr] = 0.;(*y)[itr] = 0.;(*z)[itr] = 0.;
551 (*dx)[itr] = -1.; (*dy)[itr] = 100.; (*dz)[itr] = 100.; (*dpt)[itr] = 100.;
554 // Initialize TRD track to the first track reference
555 AliTrackReference *tr(NULL);
556 Int_t itr(0); while(!(tr = fTrackRefs[itr])) itr++;
557 if(tr->Pt()<1.e-3) return kFALSE;
560 Double_t xyz[3]={tr->X(),tr->Y(),tr->Z()};
561 Double_t pxyz[3]={tr->Px(),tr->Py(),tr->Pz()};
562 Double_t var[6] = {1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4};
564 var[0], 0., 0., 0., 0., 0.,
565 var[1], 0., 0., 0., 0.,
572 const TParticlePDG *pdg=db.GetParticle(fPDG);
574 AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::PropagateKalman()", Form("PDG entry missing for code %d. References for track %d", fPDG, fNTrackRefs));
577 tt.Set(xyz, pxyz, cov, Short_t(pdg->Charge()));
578 if(mass<0){ // mass 0 use PDG
579 tt.SetMass(pdg->Mass());
580 } else { // use rec value
584 // Double_t bg(tr->P()/pdg->Mass());
585 // printf("\n\nNEW track PDG[%d] bg[%f] x[%f]\n", fPDG, bg, TMath::Log(bg));
586 Double_t x0(tr->LocalX());
587 const Double_t *cc(NULL);
588 for(Int_t ip=0; itr<fNTrackRefs; itr++){
589 if(!(tr = fTrackRefs[itr])) continue;
590 // printf("ip[%d] det[%d]\n", ip, tr->DetectorId());
591 if(!AliTRDtrackerV1::PropagateToX(tt, tr->LocalX(), fgKalmanStep)) continue;
598 (*dx)[ip] = tt.GetX() - x0;
599 (*dy)[ip] = tt.GetY() - tr->LocalY();
600 (*dz)[ip] = tt.GetZ() - tr->Z();
601 (*pt)[ip] = tr->Pt();
602 (*dpt)[ip] = tt.Pt()- tr->Pt();
603 // printf("pt : kalman[%e] MC[%e]\n", tt.Pt(), tr->Pt());
604 (*budget)[ip] = tt.GetBudget(0);
605 cc = tt.GetCovariance();
606 // printf("dx[%5.2f] sy[%f]\n", (*dx)[ip], TMath::Sqrt(cc[0]));
607 for(Int_t ic(0), jp(ip*15); ic<15; ic++, jp++) (*c)[jp]=cc[ic];