]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/info/AliTRDtrackInfo.cxx
save MC pt for Kalman propagator test
[u/mrichter/AliRoot.git] / PWG1 / TRD / info / AliTRDtrackInfo.cxx
CommitLineData
1ee39b3a 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: AliTRDtrackInfo.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// Reconstruction QA //
21// //
22// Authors: //
23// Alex Bercuci <A.Bercuci@gsi.de> //
24// Markus Fasel <M.Fasel@gsi.de> //
25// //
26////////////////////////////////////////////////////////////////////////////
27
28#include "TDatabasePDG.h"
4226db3e 29#include "TVectorT.h"
1ee39b3a 30
31#include "AliTrackReference.h"
32#include "AliExternalTrackParam.h"
33
34#include "AliTRDseedV1.h"
35#include "AliTRDtrackV1.h"
36#include "AliTRDgeometry.h"
37#include "AliTRDtrackerV1.h"
38
39#include "AliTRDtrackInfo.h"
40
41ClassImp(AliTRDtrackInfo)
42ClassImp(AliTRDtrackInfo::AliMCinfo)
43ClassImp(AliTRDtrackInfo::AliESDinfo)
5066aa9a 44Double_t AliTRDtrackInfo::AliMCinfo::fgKalmanStep = 2.;
1ee39b3a 45
46//___________________________________________________
47AliTRDtrackInfo::AliTRDtrackInfo():
48 TObject()
49 ,fNClusters(0)
f232df0d 50 ,fTRDtrack(NULL)
51 ,fMC(NULL)
1ee39b3a 52 ,fESD()
53{
54 //
55 // Default constructor
56 //
57}
58
59
60//___________________________________________________
61AliTRDtrackInfo::AliTRDtrackInfo(const AliTRDtrackInfo &trdInfo):
62 TObject((const TObject&)trdInfo)
63 ,fNClusters(trdInfo.fNClusters)
f232df0d 64 ,fTRDtrack(NULL)
1ee39b3a 65 ,fMC(new AliMCinfo(*trdInfo.fMC))
66 ,fESD(trdInfo.fESD)
67{
68 //
69 // copy Entries
70 //
71
72 if(trdInfo.fMC) fMC = new AliMCinfo(*trdInfo.fMC);
73
74 if(trdInfo.fTRDtrack){
75 fTRDtrack = new AliTRDtrackV1(*trdInfo.fTRDtrack);
76 if(trdInfo.fTRDtrack->IsOwner()) fTRDtrack->SetOwner();
77 }
78}
79
80//___________________________________________________
81AliTRDtrackInfo::AliMCinfo::AliMCinfo()
82 :fLabel(0)
83 ,fPDG(0)
84 ,fNTrackRefs(0)
85{
86 // Set 0-Pointers
87 memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
88}
89
90//___________________________________________________
91AliTRDtrackInfo::AliMCinfo::AliMCinfo(const AliMCinfo &mc)
92 :fLabel(mc.fLabel)
93 ,fPDG(mc.fPDG)
94 ,fNTrackRefs(mc.fNTrackRefs)
95{
96 //
97 // Constructor
98 //
99
100 memset(fTrackRefs, 0, sizeof(AliTrackReference *) * 12);
101 for(Int_t ien = 0; ien < 12; ien++){
102 if(mc.fTrackRefs[ien])
103 fTrackRefs[ien] = new AliTrackReference(*(mc.fTrackRefs[ien]));
104 }
105}
106
107//___________________________________________________
108AliTRDtrackInfo::AliESDinfo::AliESDinfo()
109 :fId(-1)
110 ,fStatus(0)
111 ,fKinkIndex(0)
112 ,fTPCncls(0)
113 ,fTRDpidQuality(0)
114 ,fTRDnSlices(0)
f232df0d 115 ,fTRDslices(NULL)
116 ,fOP(NULL)
1ee39b3a 117{
118 //
119 // Constructor
120 //
121
122 memset(fTRDr, 0, AliPID::kSPECIES*sizeof(Double32_t));
123}
124
125//___________________________________________________
126AliTRDtrackInfo::AliESDinfo::AliESDinfo(const AliESDinfo &esd)
127 :fId(esd.fId)
128 ,fStatus(esd.fStatus)
129 ,fKinkIndex(esd.fKinkIndex)
130 ,fTPCncls(esd.fTPCncls)
131 ,fTRDpidQuality(esd.fTRDpidQuality)
132 ,fTRDnSlices(esd.fTRDnSlices)
f232df0d 133 ,fTRDslices(NULL)
134 ,fOP(NULL)
1ee39b3a 135{
136 //
137 // Constructor
138 //
139
140 memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
141
142 if(fTRDnSlices){
143 fTRDslices = new Double32_t[fTRDnSlices];
144 memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
145 }
146 if(esd.fOP) fOP = new AliExternalTrackParam(*esd.fOP);
147}
148
149
150//___________________________________________________
151AliTRDtrackInfo::~AliTRDtrackInfo()
152{
153 //
154 // Destructor
155 //
156
157 if(fMC) delete fMC;
158 if(fTRDtrack) delete fTRDtrack;
159}
160
161//___________________________________________________
162AliTRDtrackInfo::AliMCinfo::~AliMCinfo()
163{
164 //
165 // Destructor
166 //
167
168 fNTrackRefs = 0;
169 for(Int_t ien = 0; ien < 12; ien++){
170 if(fTrackRefs[ien]) delete fTrackRefs[ien];
f232df0d 171 fTrackRefs[ien] = NULL;
1ee39b3a 172 }
173}
174
175//___________________________________________________
176AliTRDtrackInfo::AliESDinfo::~AliESDinfo()
177{
178 //
179 // Destructor
180 //
181
182 if(fTRDnSlices){
183 delete [] fTRDslices;
f232df0d 184 fTRDslices = NULL;
1ee39b3a 185 fTRDnSlices = 0;
186 }
f232df0d 187 if(fOP) delete fOP; fOP = NULL;
1ee39b3a 188}
189
190
191//___________________________________________________
192AliTRDtrackInfo& AliTRDtrackInfo::operator=(const AliTRDtrackInfo &trdInfo)
193{
194 //
195 // = Operator
196 //
197
198 fNClusters = trdInfo.fNClusters;
199 fESD = trdInfo.fESD;
200
201 if(trdInfo.fMC){
202 if(!fMC)
203 fMC = new AliMCinfo(*trdInfo.fMC);
204 else
205 new(fMC) AliMCinfo(*trdInfo.fMC);
206 }
207
208 if(trdInfo.fTRDtrack){
209 if(!fTRDtrack)
210 fTRDtrack = new AliTRDtrackV1(*trdInfo.fTRDtrack);
211 else
212 new(fTRDtrack) AliTRDtrackV1(*trdInfo.fTRDtrack);
213 if(trdInfo.fTRDtrack->IsOwner()) fTRDtrack->SetOwner();
214 }
215
216 return *this;
217}
218
219//___________________________________________________
220AliTRDtrackInfo::AliMCinfo& AliTRDtrackInfo::AliMCinfo::operator=(const AliMCinfo &mc)
221{
222 //
223 // Assignment operator
224 //
225
226 fLabel = mc.fLabel;
227 fPDG = mc.fPDG;
228 fNTrackRefs = mc.fNTrackRefs;
229
230 AliTrackReference **itr = &fTrackRefs[0];
231 AliTrackReference* const *jtr = &mc.fTrackRefs[0];
232 for(Int_t ien = 0; ien < 12; ien++, itr++, jtr++){
233 if((*jtr)){
234 if(!(*itr)) (*itr) = new AliTrackReference(*(*jtr));
235 else new(&(*itr)) AliTrackReference(*(*jtr));
f232df0d 236 } else (*itr) = NULL;
1ee39b3a 237 }
238 return *this;
239}
240
241//___________________________________________________
242AliTRDtrackInfo::AliESDinfo& AliTRDtrackInfo::AliESDinfo::operator=(const AliESDinfo &esd)
243{
244 //
245 // Assignment operator
246 //
247
248 fId = esd.fId;
249 fStatus = esd.fStatus;
250 fTRDpidQuality = esd.fTRDpidQuality;
251 fTRDnSlices = esd.fTRDnSlices;
f232df0d 252 fTRDslices = NULL;
1ee39b3a 253
254 memcpy(fTRDr, esd.fTRDr, AliPID::kSPECIES*sizeof(Double32_t));
255
256 if(fTRDnSlices){
257 fTRDslices = new Double32_t[fTRDnSlices];
258 memcpy(fTRDslices, esd.fTRDslices, fTRDnSlices*sizeof(Double32_t));
259 }
260 if(esd.fOP){
261 if(fOP) new(fOP) AliExternalTrackParam(esd.fOP);
262 else fOP = new AliExternalTrackParam(esd.fOP);
f232df0d 263 } else fOP = NULL;
1ee39b3a 264
265 return *this;
266}
267
268//___________________________________________________
269void AliTRDtrackInfo::Delete(const Option_t *)
270{
271 //
272 // Delete
273 //
274
275 fNClusters = 0;
f232df0d 276 if(fMC) delete fMC; fMC = NULL;
277 if(fTRDtrack) delete fTRDtrack; fTRDtrack = NULL;
1ee39b3a 278}
279
280//___________________________________________________
281void AliTRDtrackInfo::SetTrack(const AliTRDtrackV1 *track)
282{
283 //
284 // Set the TRD track
285 //
286
287 if(!fTRDtrack)
288 fTRDtrack = new AliTRDtrackV1(*track);
289 else
290 new(fTRDtrack)AliTRDtrackV1(*track);
291 fTRDtrack->SetOwner();
292 // Make a copy for the object in order to avoid ownership problems
293}
294
295//___________________________________________________
296void AliTRDtrackInfo::AddTrackRef(const AliTrackReference *tref)
297{
298 //
299 // Add track reference
300 //
301
302 SetMC();
303 if(fMC->fNTrackRefs >= 12){
304 SetCurved();
305 return;
306 }
307 // Make a copy for the object in order to avoid ownership problems
308 fMC->fTrackRefs[fMC->fNTrackRefs++] = new AliTrackReference(*tref);
309}
310
311//___________________________________________________
312AliTrackReference* AliTRDtrackInfo::GetTrackRef(Int_t idx) const
313{
314//
315// Returns a track reference
316//
f232df0d 317 if(!fMC) return NULL;
318 return (idx>=0 && idx < 12) ? fMC->fTrackRefs[idx] : NULL;
1ee39b3a 319}
320
321//___________________________________________________
322AliTrackReference* AliTRDtrackInfo::GetTrackRef(AliTRDseedV1* const tracklet) const
323{
324//
325// Returns a track reference
326//
f232df0d 327 if(!fMC) return NULL;
1ee39b3a 328 Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
329 AliTrackReference * const* jtr = &(fMC->fTrackRefs[0]);
330 for(Int_t itr = 0; itr < fMC->fNTrackRefs; itr++, ++jtr){
331 if(!(*jtr)) break;
332 if(TMath::Abs(tracklet->GetX0() - (*jtr)->LocalX()) < cw) return (*jtr);
333 }
f232df0d 334 return NULL;
1ee39b3a 335}
336
337//___________________________________________________
338Int_t AliTRDtrackInfo::GetNumberOfClusters() const
339{
340 //
341 // Returns the number of clusters
342 //
343
344 Int_t n = 0;
345 if(!fTRDtrack) return 0;
346 if(fTRDtrack->GetNumberOfTracklets() == 0) return n;
f232df0d 347 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 348 for(Int_t ip=0; ip<6; ip++){
349 if(!(tracklet = const_cast<AliTRDseedV1 *>(fTRDtrack->GetTracklet(ip)))) continue;
350 n+=tracklet->GetN();
351 }
352 return n;
353}
354
355
356//___________________________________________________
357void AliTRDtrackInfo::SetOuterParam(const AliExternalTrackParam *op)
358{
359 //
360 // Set outer track parameters
361 //
362
363 if(!op) return;
364 if(fESD.fOP) new(fESD.fOP) AliExternalTrackParam(*op);
365 else fESD.fOP = new AliExternalTrackParam(*op);
366}
367
368//___________________________________________________
369Int_t AliTRDtrackInfo::GetNTracklets() const
370{
371 //
372 // Return the number of tracklets
373 //
374
f232df0d 375 if(!fTRDtrack) return NULL;
1ee39b3a 376 return fTRDtrack->GetNumberOfTracklets();
377}
378
379//___________________________________________________
380void AliTRDtrackInfo::SetSlices(Int_t n, Double32_t *s)
381{
382 //
383 // Set the slices
384 //
1ee39b3a 385 if(fESD.fTRDnSlices != n){
386 fESD.fTRDnSlices = 0;
387 delete [] fESD.fTRDslices;
f232df0d 388 fESD.fTRDslices = NULL;
1ee39b3a 389 }
390
391 if(!fESD.fTRDnSlices){
392 fESD.fTRDnSlices = n;
f232df0d 393 fESD.fTRDslices = new Double32_t[fESD.fTRDnSlices];
1ee39b3a 394 }
395
396 memcpy(fESD.fTRDslices, s, n*sizeof(Double32_t));
397}
398
399//___________________________________________________
400Bool_t AliTRDtrackInfo::AliMCinfo::GetDirections(Float_t &x0, Float_t &y0, Float_t &z0, Float_t &dydx, Float_t &dzdx, Float_t &pt, UChar_t &status) const
401{
402// Check for 2 track ref for the tracklet defined bythe radial position x0
403// The "status" is a bit map and gives a more informative output in case of failure:
404// - 0 : everything is OK
405// - BIT(0) : 0 track refs found
406// - BIT(1) : 1 track reference found
407// - BIT(2) : dx <= 0 between track references
408// - BIT(3) : dx > 0 && dx < 3.7 - tangent tracks
409
410 status = 0;
411 Double_t cw = AliTRDgeometry::CamHght() + AliTRDgeometry::CdrHght();
412 Int_t nFound = 0;
f232df0d 413 AliTrackReference *tr[2] = {NULL, NULL};
1ee39b3a 414 AliTrackReference * const* jtr = &fTrackRefs[0];
415 for(Int_t itr = 0; itr < fNTrackRefs; itr++, ++jtr){
416 if(!(*jtr)) break;
417/*
418 if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, (*jtr)->LocalX());*/
419 if(TMath::Abs(x0 - (*jtr)->LocalX()) > cw) continue;
420 tr[nFound++] = (*jtr);
421 if(nFound == 2) break;
422 }
423 if(nFound < 2){
f232df0d 424 AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Missing track ref x0[%6.3f] nref[%d]", x0, nFound));
1ee39b3a 425 if(!nFound) SETBIT(status, 0);
426 else SETBIT(status, 1);
427 return kFALSE;
428 }
429 if((pt=tr[1]->Pt()) < 1.e-3) return kFALSE;
430
431 Double_t dx = tr[1]->LocalX() - tr[0]->LocalX();
432 if(dx <= 0.){
433 AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::GetDirections()", Form("Track ref with wrong radial distances refX0[%6.3f] refX1[%6.3f]", tr[0]->LocalX(), tr[1]->LocalX()));
434 SETBIT(status, 2);
435 return kFALSE;
436 }
437 if(TMath::Abs(dx-AliTRDgeometry::CamHght()-AliTRDgeometry::CdrHght())>1.E-3) SETBIT(status, 3);
438
439 dydx = (tr[1]->LocalY() - tr[0]->LocalY()) / dx;
440 dzdx = (tr[1]->Z() - tr[0]->Z()) / dx;
441 //Float_t dx0 = tr[1]->LocalX() - x0;
442 y0 = tr[1]->LocalY()/* - dydx*dx0*/;
443 z0 = tr[1]->Z()/* - dzdx*dx0*/;
444 x0 = tr[1]->LocalX();
445 return kTRUE;
446}
447
448//___________________________________________________
5066aa9a 449void AliTRDtrackInfo::AliMCinfo::PropagateKalman(TVectorD *dx, TVectorD *dy, TVectorD *dz, TVectorD *pt, TVectorD *dpt, TVectorD *c) const
1ee39b3a 450{
451// Propagate Kalman from the first TRD track reference to
452// last one and save residuals in the y, z and pt.
453//
454// This is to calibrate the dEdx and MS corrections
455
4226db3e 456 if(!fNTrackRefs) return;
1ee39b3a 457 for(Int_t itr=kNTrackRefs; itr--;){
4226db3e 458 (*dx)[itr] = -1.; (*dy)[itr] = 100.; (*dz)[itr] = 100.; (*dpt)[itr] = 100.;
1ee39b3a 459 }
1ee39b3a 460
461 // Initialize TRD track to the first track reference
4226db3e 462 AliTrackReference *tr(NULL);
463 Int_t itr(0); while(!(tr = fTrackRefs[itr])) itr++;
1ee39b3a 464 if(tr->Pt()<1.e-3) return;
465
466 AliTRDtrackV1 tt;
467 Double_t xyz[3]={tr->X(),tr->Y(),tr->Z()};
468 Double_t pxyz[3]={tr->Px(),tr->Py(),tr->Pz()};
469 Double_t var[6] = {1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4};
470 Double_t cov[21]={
471 var[0], 0., 0., 0., 0., 0.,
472 var[1], 0., 0., 0., 0.,
473 var[2], 0., 0., 0.,
474 var[3], 0., 0.,
475 var[4], 0.,
476 var[5]
477 };
478 TDatabasePDG db;
479 const TParticlePDG *pdg=db.GetParticle(fPDG);
480 if(!pdg){
481 AliWarningGeneral("AliTRDtrackInfo::AliMCinfo::PropagateKalman()", Form("PDG entry missing for code %d. References for track %d", fPDG, fNTrackRefs));
482 return;
483 }
484 tt.Set(xyz, pxyz, cov, Short_t(pdg->Charge()));
485 tt.SetMass(pdg->Mass());
486
4226db3e 487 Double_t x0(tr->LocalX());
488 const Double_t *cc(NULL);
489 for(Int_t ip=0; itr<fNTrackRefs; itr++){
490 if(!(tr = fTrackRefs[itr])) continue;
5066aa9a 491 if(!AliTRDtrackerV1::PropagateToX(tt, tr->LocalX(), fgKalmanStep)) continue;
1ee39b3a 492
493 //if(update) ...
4226db3e 494 (*dx)[ip] = tt.GetX() - x0;
495 (*dy)[ip] = tt.GetY() - tr->LocalY();
496 (*dz)[ip] = tt.GetZ() - tr->Z();
5066aa9a 497 (*pt)[ip] = tr->Pt();
4226db3e 498 (*dpt)[ip] = tt.Pt()- tr->Pt();
1ee39b3a 499 cc = tt.GetCovariance();
4226db3e 500 c->Use(ip*15, (ip+1)*15, cc);
1ee39b3a 501 ip++;
502 }
503}