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 **************************************************************************/
18 #include "AliESDtrack.h"
19 #include "AliTracker.h"
21 #include "AliTRDtrackV1.h"
22 #include "AliTRDcluster.h"
23 #include "AliTRDcalibDB.h"
24 #include "AliTRDrecoParam.h"
26 ClassImp(AliTRDtrackV1)
28 ///////////////////////////////////////////////////////////////////////////////
30 // Represents a reconstructed TRD track //
31 // Local TRD Kalman track //
34 // Alex Bercuci <A.Bercuci@gsi.de> //
35 // Markus Fasel <M.Fasel@gsi.de> //
37 ///////////////////////////////////////////////////////////////////////////////
39 //_______________________________________________________________
40 AliTRDtrackV1::AliTRDtrackV1()
44 // Default constructor
47 for(int ip=0; ip<6; ip++){
48 fTrackletIndex[ip] = -1;
49 fTracklet[ip].Reset();
53 //_______________________________________________________________
54 AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t)
58 // Constructor from AliESDtrack
61 //AliInfo(Form("alpha %f", GetAlpha()));
62 t.GetTRDtracklets(&fTrackletIndex[0]);
63 for(int ip=0; ip<6; ip++) fTracklet[ip].Reset();
66 //_______________________________________________________________
67 AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref)
74 for(int ip=0; ip<6; ip++){
75 fTrackletIndex[ip] = ref.fTrackletIndex[ip];
76 fTracklet[ip] = ref.fTracklet[ip];
80 //_______________________________________________________________
81 // AliTRDtrackV1::~AliTRDtrackV1()
86 //_______________________________________________________________
87 AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
88 , Double_t x, Double_t alpha)
92 // The stand alone tracking constructor
93 // TEMPORARY !!!!!!!!!!!
95 // 1. covariance matrix
96 // 2. dQdl calculation
99 Double_t cnv = 1.0 / (GetBz() * kB2C);
101 Double_t pp[5] = { p[0]
107 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
108 Double_t c32 = x*cov[13] - cov[ 8];
109 Double_t c20 = x*cov[10] - cov[ 3];
110 Double_t c21 = x*cov[11] - cov[ 4];
111 Double_t c42 = x*cov[14] - cov[12];
113 Double_t cc[15] = { cov[ 0]
116 , cov[ 6], cov[ 7], c32, cov[ 9]
117 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
120 Double_t s = GetSnp();
121 Double_t t = GetTgl();
123 Int_t ncl = 0, nclPlane; AliTRDcluster *c = 0x0;
124 for(int iplane=0; iplane<6; iplane++){
125 fTrackletIndex[iplane] = -1;
126 fTracklet[iplane] = trklts[iplane];
128 for(int ic = 0; ic<AliTRDseed::knTimebins; ic++){
129 if(!fTracklet[iplane].IsUsable(ic)) continue;
130 if(!(c = fTracklet[iplane].GetClusters(ic))) continue;
132 fIndex[ncl] = fTracklet[iplane].GetIndexes(ic);
133 Double_t q = TMath::Abs(c->GetQ());
136 // This is not correctly. Has to be updated in the AliTRDtrackerV1::FollowBackProlonagation()
137 fdQdl[ncl] = q * (s*s < 1.) ? TMath::Sqrt((1-s*s)/(1+t*t)) : 1.;
141 //printf("%d N clusters plane %d [%d %d].\n", iplane, nclPlane, fTracklet[iplane].GetN2(), trklts[iplane].GetN());
143 //printf("N clusters in AliTRDtrackV1 %d.\n", ncl);
144 SetNumberOfClusters(/*ncl*/);
147 //_______________________________________________________________
148 Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
150 // set MC label for this tracklet
152 Int_t s[kMAXCLUSTERSPERTRACK][2];
153 for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
160 AliTRDcluster *c = 0x0;
161 for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
162 if(fTrackletIndex[ip] < 0) continue;
163 for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
164 if(!(c = fTracklet[ip].GetClusters(ic))) continue;
165 for (Int_t k = 0; k < 3; k++) {
166 label = c->GetLabel(k);
170 while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
171 if ((s[j][0] == label) ||
186 for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
187 if (s[i][1] <= max) continue;
192 if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
199 //_______________________________________________________________
200 Bool_t AliTRDtrackV1::CookPID()
203 // Cook the PID information
206 // CookdEdx(); // truncated mean ... do we still need it ?
208 // CookdEdxTimBin(seed->GetID());
210 // Sets the a priori probabilities
211 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
212 fPID[ispec] = 1.0 / AliPID::kSPECIES;
215 // steer PID calculation @ tracklet level
216 Double_t *prob = 0x0;
218 for(int itrklt=0; itrklt<AliESDtrack::kTRDnPlanes; itrklt++){
219 //for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) fdEdxPlane[itrklt][iSlice] = -1.;
221 if(fTrackletIndex[itrklt]<0) continue;
222 if(!fTracklet[itrklt].IsOK()) continue;
223 if(!(prob = fTracklet[itrklt].GetProbability())) return kFALSE;
225 Int_t nspec = 0; // quality check of tracklet dEdx
226 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
227 if(prob[ispec] < 0.) continue;
228 fPID[ispec] *= prob[ispec];
236 // no tracklet found for PID calculations
237 if(!fPIDquality) return kTRUE;
239 // slot for PID calculation @ track level
242 // normalize probabilities
243 Double_t probTotal = 0.0;
244 for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
247 if (probTotal <= 0.0) {
248 AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
252 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
257 //_______________________________________________________________
258 Float_t AliTRDtrackV1::GetMomentum(Int_t plane) const
261 // Get the momentum at a given plane
264 return plane >=0 && plane < 6 && fTrackletIndex[plane] >= 0 ? fTracklet[plane].GetMomentum() : -1.;
267 //_______________________________________________________________
268 Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
271 // Get the predicted chi2
274 Double_t x = trklt->GetX0();
275 Double_t p[2] = { trklt->GetYat(x)
276 , trklt->GetZat(x) };
278 trklt->GetCovAt(x, cov);
280 return AliExternalTrackParam::GetPredictedChi2(p, cov);
284 //_______________________________________________________________
285 Bool_t AliTRDtrackV1::IsOwner() const
288 // Check whether track owns the tracklets
291 for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
292 if(fTrackletIndex[ip] < 0) continue;
293 if(!fTracklet[ip].IsOwner()) return kFALSE;
298 //_____________________________________________________________________________
299 void AliTRDtrackV1::MakeBackupTrack()
302 // Creates a backup track
306 fBackupTrack->~AliTRDtrack();
307 new(fBackupTrack) AliTRDtrack((AliTRDtrack&)(*this));
309 fBackupTrack = new AliTRDtrack((AliTRDtrack&)(*this));
313 //___________________________________________________________
314 void AliTRDtrackV1::SetNumberOfClusters()
316 // Calculate the number of clusters attached to this track
319 for(int ip=0; ip<6; ip++){
320 if(fTrackletIndex[ip] >= 0) ncls += fTracklet[ip].GetN();
322 AliKalmanTrack::SetNumberOfClusters(ncls);
326 //_______________________________________________________________
327 void AliTRDtrackV1::SetOwner(Bool_t own)
330 // Toggle ownership of tracklets
333 for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
334 if(fTrackletIndex[ip] < 0) continue;
335 //AliInfo(Form("p[%d] index[%d]", ip, fTrackletIndex[ip]));
336 fTracklet[ip].SetOwner(own);
340 //_______________________________________________________________
341 void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
347 if(plane < 0 || plane >= AliESDtrack::kTRDnPlanes) return;
348 fTracklet[plane] = (*trklt);
349 fTrackletIndex[plane] = index;
352 //_______________________________________________________________
353 Bool_t AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
356 // Update track parameters
359 Double_t x = trklt->GetX0();
360 Double_t p[2] = { trklt->GetYat(x)
361 , trklt->GetZat(x) };
363 trklt->GetCovAt(x, cov);
366 // AliInfo(Form("cov[%f %f %f]", cov[0], cov[1], cov[2]));
368 if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
371 AliTRDcluster *c = 0x0;
372 Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
373 AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
376 // Register info to track
377 // Int_t n = GetNumberOfClusters();
378 // fIndex[n] = index;
380 SetNumberOfClusters(/*GetNumberOfClusters()+trklt->GetN()*/);
381 SetChi2(GetChi2() + chisq);
384 trklt->SetMomentum(GetP());
385 trklt->SetSnp(GetSnp());
386 trklt->SetTgl(GetTgl());
390 //_______________________________________________________________
391 void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
394 // Update the TRD PID information in the ESD track
397 track->SetNumberOfTRDslices(kNslice);
399 for (Int_t ip = 0; ip < kNplane; ip++) {
400 if(fTrackletIndex[ip] < 0) continue;
401 fTracklet[ip].CookdEdx(kNslice);
402 Float_t *dedx = fTracklet[ip].GetdEdx();
403 for (Int_t js = 0; js < kNslice; js++) {
404 track->SetTRDslice(dedx[js], ip, js);
409 track->SetTRDpid(fPID);
410 track->SetTRDpidQuality(fPIDquality);