]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtrackV1.cxx
Copy arrays in assignment instead of the pointer; avoid double delete.
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackV1.cxx
CommitLineData
d9950a5a 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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/* $Id$ */
17
eb38ed55 18#include "AliESDtrack.h"
41880fac 19#include "AliTracker.h"
eb38ed55 20
d9950a5a 21#include "AliTRDtrackV1.h"
22#include "AliTRDcluster.h"
23#include "AliTRDcalibDB.h"
24#include "AliTRDrecoParam.h"
25
d9950a5a 26ClassImp(AliTRDtrackV1)
27
0906e73e 28///////////////////////////////////////////////////////////////////////////////
29// //
30// Represents a reconstructed TRD track //
31// Local TRD Kalman track //
32// //
33// Authors: //
34// Alex Bercuci <A.Bercuci@gsi.de> //
35// Markus Fasel <M.Fasel@gsi.de> //
36// //
37///////////////////////////////////////////////////////////////////////////////
d9950a5a 38
39//_______________________________________________________________
0906e73e 40AliTRDtrackV1::AliTRDtrackV1()
41 :AliTRDtrack()
0906e73e 42{
d9950a5a 43 //
44 // Default constructor
45 //
0906e73e 46
d9950a5a 47 for(int ip=0; ip<6; ip++){
48 fTrackletIndex[ip] = -1;
49 fTracklet[ip].Reset();
50 }
51}
52
53//_______________________________________________________________
0906e73e 54AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t)
55 :AliTRDtrack(t)
d9950a5a 56{
57 //
0906e73e 58 // Constructor from AliESDtrack
d9950a5a 59 //
60
61 //AliInfo(Form("alpha %f", GetAlpha()));
62 t.GetTRDtracklets(&fTrackletIndex[0]);
63 for(int ip=0; ip<6; ip++) fTracklet[ip].Reset();
64}
65
66//_______________________________________________________________
0906e73e 67AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref)
68 :AliTRDtrack(ref)
d9950a5a 69{
70 //
71 // Copy constructor
72 //
73
0906e73e 74 for(int ip=0; ip<6; ip++){
75 fTrackletIndex[ip] = ref.fTrackletIndex[ip];
76 fTracklet[ip] = ref.fTracklet[ip];
77 }
d9950a5a 78}
79
80//_______________________________________________________________
81// AliTRDtrackV1::~AliTRDtrackV1()
82// {
83//
84// }
d9950a5a 85
86//_______________________________________________________________
eb38ed55 87AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
88 , Double_t x, Double_t alpha)
0906e73e 89 :AliTRDtrack()
d9950a5a 90{
0906e73e 91 //
92 // The stand alone tracking constructor
93 // TEMPORARY !!!!!!!!!!!
94 // to check :
95 // 1. covariance matrix
96 // 2. dQdl calculation
97 //
d9950a5a 98
0906e73e 99 Double_t cnv = 1.0 / (GetBz() * kB2C);
d9950a5a 100
101 Double_t pp[5] = { p[0]
102 , p[1]
103 , x*p[4] - p[2]
104 , p[3]
105 , p[4]*cnv };
106
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];
112
113 Double_t cc[15] = { cov[ 0]
114 , cov[ 1], cov[ 2]
115 , c20, c21, c22
116 , cov[ 6], cov[ 7], c32, cov[ 9]
117 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
118
119 Set(x,alpha,pp,cc);
120 Double_t s = GetSnp();
121 Double_t t = GetTgl();
122
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];
127 nclPlane = 0;
128 for(int ic = 0; ic<AliTRDseed::knTimebins; ic++){
129 if(!fTracklet[iplane].IsUsable(ic)) continue;
130 if(!(c = fTracklet[iplane].GetClusters(ic))) continue;
131
132 fIndex[ncl] = fTracklet[iplane].GetIndexes(ic);
133 Double_t q = TMath::Abs(c->GetQ());
134 fClusters[ncl] = c;
135 // temporary !!!
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.;
138 ncl++;
139 nclPlane++;
140 }
141 //printf("%d N clusters plane %d [%d %d].\n", iplane, nclPlane, fTracklet[iplane].GetN2(), trklts[iplane].GetN());
142 }
143 //printf("N clusters in AliTRDtrackV1 %d.\n", ncl);
eb38ed55 144 SetNumberOfClusters(/*ncl*/);
d9950a5a 145}
146
bb56afff 147//_______________________________________________________________
148Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
149{
150 // set MC label for this tracklet
151
152 Int_t s[kMAXCLUSTERSPERTRACK][2];
153 for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
154 s[i][0] = -1;
155 s[i][1] = 0;
156 }
157
158 Bool_t labelAdded;
159 Int_t label;
160 AliTRDcluster *c = 0x0;
6984f7c1 161 for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
bb56afff 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);
167 labelAdded = kFALSE;
168 Int_t j = 0;
169 if (label >= 0) {
170 while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
171 if ((s[j][0] == label) ||
172 (s[j][1] == 0)) {
173 s[j][0] = label;
174 s[j][1]++;
175 labelAdded = kTRUE;
176 }
177 j++;
178 }
179 }
180 }
181 }
182 }
183
184 Int_t max = 0;
185 label = -123456789;
186 for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
187 if (s[i][1] <= max) continue;
188 max = s[i][1];
189 label = s[i][0];
190 }
191
192 if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
193
194 SetLabel(label);
195
196 return kTRUE;
197}
198
d9950a5a 199//_______________________________________________________________
200Bool_t AliTRDtrackV1::CookPID()
201{
0906e73e 202 //
203 // Cook the PID information
204 //
205
d9950a5a 206// CookdEdx(); // truncated mean ... do we still need it ?
207
208// CookdEdxTimBin(seed->GetID());
0906e73e 209
210 // Sets the a priori probabilities
211 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
212 fPID[ispec] = 1.0 / AliPID::kSPECIES;
213 }
214
215 // steer PID calculation @ tracklet level
216 Double_t *prob = 0x0;
217 fPIDquality = 0;
6984f7c1 218 for(int itrklt=0; itrklt<AliESDtrack::kTRDnPlanes; itrklt++){
0906e73e 219 //for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) fdEdxPlane[itrklt][iSlice] = -1.;
d9950a5a 220
221 if(fTrackletIndex[itrklt]<0) continue;
6888a49b 222 if(!fTracklet[itrklt].IsOK()) continue;
0906e73e 223 if(!(prob = fTracklet[itrklt].GetProbability())) return kFALSE;
224
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];
229 nspec++;
230 }
231 if(!nspec) continue;
232
233 fPIDquality++;
d9950a5a 234 }
0906e73e 235
236 // no tracklet found for PID calculations
237 if(!fPIDquality) return kTRUE;
238
239 // slot for PID calculation @ track level
240
241
242 // normalize probabilities
243 Double_t probTotal = 0.0;
244 for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
d9950a5a 245
d9950a5a 246
0906e73e 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.");
d9950a5a 249 return kFALSE;
250 }
251
0906e73e 252 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
d9950a5a 253
0906e73e 254 return kTRUE;
d9950a5a 255}
256
0906e73e 257//_______________________________________________________________
258Float_t AliTRDtrackV1::GetMomentum(Int_t plane) const
259{
260 //
261 // Get the momentum at a given plane
262 //
263
264 return plane >=0 && plane < 6 && fTrackletIndex[plane] >= 0 ? fTracklet[plane].GetMomentum() : -1.;
265}
d9950a5a 266
267//_______________________________________________________________
268Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
269{
270 //
0906e73e 271 // Get the predicted chi2
d9950a5a 272 //
273
274 Double_t x = trklt->GetX0();
275 Double_t p[2] = { trklt->GetYat(x)
276 , trklt->GetZat(x) };
277 Double_t cov[3];
278 trklt->GetCovAt(x, cov);
279
280 return AliExternalTrackParam::GetPredictedChi2(p, cov);
281
282}
283
0906e73e 284//_______________________________________________________________
285Bool_t AliTRDtrackV1::IsOwner() const
286{
287 //
288 // Check whether track owns the tracklets
289 //
290
6984f7c1 291 for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
0906e73e 292 if(fTrackletIndex[ip] < 0) continue;
293 if(!fTracklet[ip].IsOwner()) return kFALSE;
294 }
295 return kTRUE;
296}
297
87a7fa94 298//_____________________________________________________________________________
299void AliTRDtrackV1::MakeBackupTrack()
300{
301 //
302 // Creates a backup track
303 //
304
305 if (fBackupTrack) {
306 fBackupTrack->~AliTRDtrack();
307 new(fBackupTrack) AliTRDtrack((AliTRDtrack&)(*this));
308 }
309 fBackupTrack = new AliTRDtrack((AliTRDtrack&)(*this));
310}
311
312
eb38ed55 313//___________________________________________________________
314void AliTRDtrackV1::SetNumberOfClusters()
315{
316// Calculate the number of clusters attached to this track
317
318 Int_t ncls = 0;
319 for(int ip=0; ip<6; ip++){
320 if(fTrackletIndex[ip] >= 0) ncls += fTracklet[ip].GetN();
321 }
322 AliKalmanTrack::SetNumberOfClusters(ncls);
323}
324
325
0906e73e 326//_______________________________________________________________
327void AliTRDtrackV1::SetOwner(Bool_t own)
328{
329 //
330 // Toggle ownership of tracklets
331 //
332
6984f7c1 333 for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
0906e73e 334 if(fTrackletIndex[ip] < 0) continue;
335 //AliInfo(Form("p[%d] index[%d]", ip, fTrackletIndex[ip]));
336 fTracklet[ip].SetOwner(own);
337 }
338}
339
d9950a5a 340//_______________________________________________________________
341void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
342{
343 //
0906e73e 344 // Set the tracklets
d9950a5a 345 //
346
6984f7c1 347 if(plane < 0 || plane >= AliESDtrack::kTRDnPlanes) return;
d9950a5a 348 fTracklet[plane] = (*trklt);
349 fTrackletIndex[plane] = index;
350}
351
352//_______________________________________________________________
0906e73e 353Bool_t AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
d9950a5a 354{
355 //
0906e73e 356 // Update track parameters
d9950a5a 357 //
358
0906e73e 359 Double_t x = trklt->GetX0();
d9950a5a 360 Double_t p[2] = { trklt->GetYat(x)
361 , trklt->GetZat(x) };
362 Double_t cov[3];
363 trklt->GetCovAt(x, cov);
364
365// Print();
366// AliInfo(Form("cov[%f %f %f]", cov[0], cov[1], cov[2]));
367
368 if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
369 //Print();
370
41880fac 371 AliTRDcluster *c = 0x0;
372 Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
373 AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
374
375
d9950a5a 376 // Register info to track
377// Int_t n = GetNumberOfClusters();
378// fIndex[n] = index;
379// fClusters[n] = c;
eb38ed55 380 SetNumberOfClusters(/*GetNumberOfClusters()+trklt->GetN()*/);
d9950a5a 381 SetChi2(GetChi2() + chisq);
0906e73e 382
383 // update tracklet
384 trklt->SetMomentum(GetP());
bcb6fb78 385 trklt->SetSnp(GetSnp());
386 trklt->SetTgl(GetTgl());
d9950a5a 387 return kTRUE;
388}
389
390//_______________________________________________________________
0906e73e 391void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
d9950a5a 392{
393 //
6984f7c1 394 // Update the TRD PID information in the ESD track
d9950a5a 395 //
6984f7c1 396
397 track->SetNumberOfTRDslices(kNslice);
0906e73e 398
6984f7c1 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);
405 }
406 }
0906e73e 407
6984f7c1 408 // copy PID to ESD
409 track->SetTRDpid(fPID);
410 track->SetTRDpidQuality(fPIDquality);
d9950a5a 411}