]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/global/AliHLTGlobalTrackMerger.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing obsolete mapping macro
[u/mrichter/AliRoot.git] / HLT / global / AliHLTGlobalTrackMerger.cxx
CommitLineData
ec6160d5 1//$Id$
2
3/**************************************************************************
4 * This file is property of and copyright by the ALICE HLT Project *
5 * ALICE Experiment at CERN, All rights reserved. *
6 * *
7 * Primary Authors: Jacek Otwinowski (Jacek.Otwinowski@gsi.de) *
8 * for The ALICE HLT Project. *
9 * *
10 * Permission to use, copy, modify and distribute this software and its *
11 * documentation strictly for non-commercial purposes is hereby granted *
12 * without fee, provided that the above copyright notice appears in all *
13 * copies and that both the copyright notice and this permission notice *
14 * appear in the supporting documentation. The authors make no claims *
15 * about the suitability of this software for any purpose. It is *
16 * provided "as is" without express or implied warranty. *
17 **************************************************************************/
18
19/** @file AliHLTGlobalTrackMerger.cxx
20 @author Jacek Otwinowski
21 @date
22 @brief The HLT global merger base class
23*/
24
25#include "AliHLTTPCTransform.h"
26#include "AliHLTTPCTrack.h"
27#include "AliHLTTPCTrackSegmentData.h"
28#include "AliHLTTPCTransform.h"
29#include "AliHLTTPCTrackArray.h"
30
a9f060e5 31#include "AliTPCReconstructor.h"
32
ec6160d5 33#include "AliTRDtrackV1.h"
34#include "AliESDEvent.h"
35#include "AliESDVertex.h"
36#include "AliTracker.h"
a9f060e5 37#include "TTreeStream.h"
ec6160d5 38
39#include <TClonesArray.h>
40
41#include "AliHLTGlobalTrackMerger.h"
42
43#if __GNUC__>= 3
44using namespace std;
45#endif
46
47ClassImp(AliHLTGlobalTrackMerger)
48
49AliHLTGlobalTrackMerger::AliHLTGlobalTrackMerger()
50 :
51 fMaxY(0.0),
52 fMaxZ(0.0),
53 fMaxSnp(0.0),
54 fMaxTgl(0.0),
55 fMaxSigned1Pt(0.0),
a9f060e5 56 fVertex(0),
57 fDebugStreamer(0)
ec6160d5 58{
59 //Default constructor
60
61 // standard vertex settings at the moment
62 // V(0.,0.,0.), sigmaVx=sigmaVy=5.e-3 [cm], sigmaVz=5.3 [cm]
63 fVertex = new AliESDVertex;
a9f060e5 64
65 if (AliTPCReconstructor::StreamLevel()>0) {
66 fDebugStreamer = new TTreeSRedirector("debugGlobalMerger.root");
67 }
68
ec6160d5 69}
70
71//_____________________________________________________________________________
72AliHLTGlobalTrackMerger::~AliHLTGlobalTrackMerger()
73{
74 //Destructor
75 if(fVertex) delete fVertex; fVertex =0;
a9f060e5 76 if(fDebugStreamer) delete fDebugStreamer; fDebugStreamer =0;
ec6160d5 77}
78
79//_____________________________________________________________________________
80Bool_t AliHLTGlobalTrackMerger::LoadTracks(TClonesArray *aTRDTracks, AliESDEvent *esdEvent)
81{
82 // load TRD tracks
83 if(!aTRDTracks) return kFALSE;
84
85 Int_t entries = aTRDTracks->GetEntriesFast();
86 for(Int_t i=0; i<entries; ++i) {
87 AliTRDtrackV1 *track = (AliTRDtrackV1*)aTRDTracks->At(i);
88 if(!track) continue;
89
90 FillTRDESD(track,AliESDtrack::kTRDin,esdEvent);
91 }
92
93return kTRUE;
94}
95
96//_____________________________________________________________________________
97Bool_t AliHLTGlobalTrackMerger::LoadTracks(AliHLTTPCTrackArray *aTPCTracks, AliESDEvent *esdEvent)
98{
99 // load TPC tracks
100 if(!aTPCTracks) return kFALSE;
101
a9f060e5 102 for (Int_t i=0; i<aTPCTracks->GetNTracks();++i) {
ec6160d5 103 AliHLTTPCTrack* track=(*aTPCTracks)[i];
104 if(!track) continue;
105
106 // convert to the AliKalmanTrack
107 Int_t local=track->Convert2AliKalmanTrack();
108 if(local<0) continue;
109
110 FillTPCESD(track,AliESDtrack::kTPCin,esdEvent);
111 }
112
113return kTRUE;
114}
115
116//_____________________________________________________________________________
117void AliHLTGlobalTrackMerger::FillTPCESD(AliHLTTPCTrack *tpcTrack, ULong_t flags, AliESDEvent* esdEvent)
118{
119 // create AliESDtracks from AliHLTTPCTracks
120 // and add them to AliESDEvent
121
122 if(!tpcTrack) return;
123 if(!esdEvent) return;
124
125 AliESDtrack iotrack;
126 iotrack.UpdateTrackParams(tpcTrack,flags);
127
128 Float_t points[4]={tpcTrack->GetFirstPointX(),tpcTrack->GetFirstPointY(),tpcTrack->GetLastPointX(),tpcTrack->GetLastPointY()};
129 if(tpcTrack->GetSector() == -1){ // Set first and last points for global tracks
130 Double_t s = TMath::Sin( tpcTrack->GetAlpha() );
131 Double_t c = TMath::Cos( tpcTrack->GetAlpha() );
132 points[0] = tpcTrack->GetFirstPointX()*c + tpcTrack->GetFirstPointY()*s;
133 points[1] = -tpcTrack->GetFirstPointX()*s + tpcTrack->GetFirstPointY()*c;
134 points[2] = tpcTrack->GetLastPointX() *c + tpcTrack->GetLastPointY() *s;
135 points[3] = -tpcTrack->GetLastPointX() *s + tpcTrack->GetLastPointY() *c;
136 }
137 iotrack.SetTPCPoints(points);
138
139 esdEvent->AddTrack(&iotrack);
140}
141
142//_____________________________________________________________________________
143void AliHLTGlobalTrackMerger::FillTRDESD(AliTRDtrackV1* trdTrack, ULong_t flags, AliESDEvent* esdEvent)
144{
145 // create AliESDtracks from AliTRDtrackV1
146 // and add them to AliESDEvent
147
148 if(!trdTrack) return;
149 if(!esdEvent) return;
150
151 AliESDtrack iotrack;
152 iotrack.UpdateTrackParams(trdTrack,flags);
153 esdEvent->AddTrack(&iotrack);
154}
155
156//_____________________________________________________________________________
157Bool_t AliHLTGlobalTrackMerger::Merge(AliESDEvent* esdEvent)
158{
159 // merge TPC and TRD tracks
160 // 1. propagate TPC track to the radius between TPC and TRD
161 // 2. propagate TRD track to the same radius between TPC and TRD
162 // 3. matches TPC and TRD tracks at the radius
4b6d0cc6 163 // 4. propagate matched TRD track to the merging radius (first measured TPC point - x coordinate)
ec6160d5 164 // 5. merge TPC and TRD track parameters at the merging radius
165 // 6. create AliESDtrack from merged tracks
166 // 7. add AliESDtrack to AliESDEvent
167
168 if(!esdEvent) return kFALSE;
169
170 const Double_t kMaxStep = 10.0; // [cm] track propagation step
171 const Double_t kMatchRadius = 285.0; // [cm] matching at radius between TPC and TRD
172 Double_t kMergeRadius = 0.0;
173 Bool_t isOk = kFALSE;
174 Bool_t isMatched = kFALSE;
a9f060e5 175 AliESDtrack *track=0;
4b6d0cc6 176 AliExternalTrackParam *extTPCTrack = 0;
ec6160d5 177
178 Int_t nTracks = esdEvent->GetNumberOfTracks();
ec6160d5 179 HLTWarning("nTracks %d",nTracks);
180
a9f060e5 181 Int_t nTracksTPC =0;
ec6160d5 182 for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack)
183 {
a9f060e5 184 track = esdEvent->GetTrack(iTrack);
ec6160d5 185 if(!track) continue;
186
187 // TPC tracks
188 if((track->GetStatus()&AliESDtrack::kTPCin)==0) continue;
189 AliESDtrack *tpcTrack = track;
a9f060e5 190 nTracksTPC++;
ec6160d5 191
4b6d0cc6 192 // create external tpc track param (needed to propagate to matching radius)
193 if ((extTPCTrack = new AliExternalTrackParam(*tpcTrack)) == 0) continue;
194
a9f060e5 195 kMergeRadius = tpcTrack->GetTPCPoints(0); // [cm] merging at first measured TPC point
ec6160d5 196
a9f060e5 197 HLTInfo("-------------------------------------------------------------------------------------");
4b6d0cc6 198 //HLTInfo("-----tpc track before matching: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
199 HLTInfo("-----ext tpc track before matching: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
ec6160d5 200
201 // propagate tracks to the matching radius
4b6d0cc6 202 //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
203 isOk = AliTracker::PropagateTrackTo(extTPCTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
ec6160d5 204 if(!isOk) continue;
205
4b6d0cc6 206 HLTInfo("-----ext tpc track at matching radius: alpha %f, x %f, y, %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
ec6160d5 207
a9f060e5 208 Int_t nTracksTRD =0;
ec6160d5 209 for(Int_t jTrack = 0; jTrack<nTracks; ++jTrack)
210 {
a9f060e5 211 track = esdEvent->GetTrack(jTrack);
ec6160d5 212 if(!track) continue;
213
214 // TRD tracks
215 if((track->GetStatus()&AliESDtrack::kTRDin)==0) continue;
216 AliESDtrack *trdTrack = track;
a9f060e5 217 nTracksTRD++;
218
4b6d0cc6 219 HLTInfo("-----trd track before matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5 220
221 isOk = AliTracker::PropagateTrackTo(trdTrack,kMatchRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
222 if(!isOk) continue;
223
4b6d0cc6 224 HLTInfo("-----trd track at matching radius: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5 225
226 // match TPC and TRD tracks
4b6d0cc6 227 //isMatched = MatchTracks(tpcTrack,trdTrack);
228 isMatched = MatchTracks(extTPCTrack,trdTrack);
ec6160d5 229 if(!isMatched) continue;
230
4b6d0cc6 231 //HLTInfo("-----tpc track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
232 HLTInfo("-----ext tpc track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",extTPCTrack->GetAlpha(),extTPCTrack->GetX(),extTPCTrack->GetY(),extTPCTrack->GetZ(),extTPCTrack->GetSnp(),extTPCTrack->GetTgl(),extTPCTrack->GetSigned1Pt());
233 HLTInfo("-----trd track after matching: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5 234
4b6d0cc6 235 // propagate TRD tracks to the merging radius
236 //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMergeRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
237 //if(!isOk) continue;
ec6160d5 238 isOk = AliTracker::PropagateTrackTo(trdTrack,kMergeRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
239 if(!isOk) continue;
240
4b6d0cc6 241 HLTInfo("-----tpc before merging: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
242 HLTInfo("-----trd before merging: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",trdTrack->GetAlpha(),trdTrack->GetX(),trdTrack->GetY(),trdTrack->GetZ(),trdTrack->GetSnp(),trdTrack->GetTgl(),trdTrack->GetSigned1Pt());
ec6160d5 243
244 // merge TPC and TRD tracks
245 // create AliESDtrack and add it to AliESDevent
246 Bool_t isMerged = MergeTracks(tpcTrack,trdTrack,esdEvent);
a9f060e5 247 if(!isMerged)
248 HLTInfo("No merged tracks");
ec6160d5 249 }
a9f060e5 250 HLTInfo("nTracksTRD %d",nTracksTRD);
4b6d0cc6 251
252 // delete external TPC track
253 if(extTPCTrack) delete extTPCTrack; extTPCTrack=0;
ec6160d5 254 }
a9f060e5 255 HLTInfo("nTracksTPC %d",nTracksTPC);
ec6160d5 256
257return kTRUE;
258}
259
260//_____________________________________________________________________________
261Bool_t AliHLTGlobalTrackMerger::MergeTracks(AliESDtrack *tpcTrack, AliESDtrack* trdTrack, AliESDEvent *esdEvent)
262{
263 // merge TPC and TRD track parameters
264 // create new AliESDtrack with TPC+TRD merged track parameters
265 // add AliESDtrack to AliESDEvent
266
267 if(!tpcTrack) return kFALSE;
268 if(!trdTrack) return kFALSE;
269
ec6160d5 270 Double_t trackParam[5], trackCovar[15];
271 Double_t trackChi2;
272 Int_t trackNDF;
273
274 // calculate merged track parameters
275 Bool_t isNotOK = SmoothTracks(tpcTrack->GetParameter(), tpcTrack->GetCovariance(), tpcTrack->GetTPCchi2(), 5,
276 trdTrack->GetParameter(), trdTrack->GetCovariance(), trdTrack->GetTRDchi2(), 5,
277 trackParam, trackCovar, trackChi2, trackNDF,5);
278
279 if(isNotOK)
280 return kFALSE;
281
282 //
283 // create AliESDtrack
284 // merged TPC+TRD information
285 //
286
287 AliESDtrack track;
288 //track.UpdateTrackParams(tpcTrack, AliESDtrack::kTPCrefit);
289 track.SetStatus(AliESDtrack::kGlobalMerge);
290 track.SetLabel(tpcTrack->GetLabel());
291 track.Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),trackParam,trackCovar);
292 track.SetGlobalChi2(trackChi2);
293
294 //track.SetTPCLabel(tpcTrack->GetLabel());
295 Double32_t tpcPID[AliPID::kSPECIES];
296 tpcTrack->GetTPCpid(tpcPID);
297 track.SetTPCpid(tpcPID);
298 //fTPCncls=t->GetNumberOfClusters(); // no cluster on HLT
299 //fTPCchi2=t->GetChi2();
300
301 //track.SetTRDLabel(trdTrack->GetLabel());
302 Double32_t trdPID[AliPID::kSPECIES];
303 trdTrack->GetTRDpid(trdPID);
304 track.SetTRDpid(trdPID);
305 //fTRDchi2 = t->GetChi2();
306 //fTRDncls = t->GetNumberOfClusters();
307 //for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
308
309 // add track to AliESDEvent
310 esdEvent->AddTrack(&track);
311
312return kTRUE;
313}
314
315//_____________________________________________________________________________
316void AliHLTGlobalTrackMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxsnp, Double_t maxtgl, Double_t signed1Pt)
317{
318 //set parameters for merger
319 fMaxY = maxy;
320 fMaxZ = maxz;
321 fMaxSnp = maxsnp;
322 fMaxTgl = maxtgl;
323 fMaxSigned1Pt = signed1Pt;
324}
325
326//_____________________________________________________________________________
4b6d0cc6 327Bool_t AliHLTGlobalTrackMerger::MatchTracks(AliExternalTrackParam *trackTPC, AliESDtrack *trackTRD)
ec6160d5 328{
329 // match TPC and TRD tracks
330 // return kTRUE in case of matching
331
332 if(!trackTPC) return kFALSE;
333 if(!trackTRD) return kFALSE;
334
a9f060e5 335 Double_t x_tpc=trackTPC->GetX();
336 Double_t y_tpc=trackTPC->GetY();
337 Double_t z_tpc=trackTPC->GetZ();
338 Double_t snp_tpc=trackTPC->GetSnp();
339 Double_t tgl_tpc=trackTPC->GetTgl();
340 Double_t signed1Pt_tpc=trackTPC->GetSigned1Pt();
341
342 Double_t x_trd=trackTRD->GetX();
343 Double_t y_trd=trackTRD->GetY();
344 Double_t z_trd=trackTRD->GetZ();
345 Double_t snp_trd=trackTRD->GetSnp();
346 Double_t tgl_trd=trackTRD->GetTgl();
347 Double_t signed1Pt_trd=trackTRD->GetSigned1Pt();
348
349 // debug stream
350 if (AliTPCReconstructor::StreamLevel()>0) {
351 //TTreeSRedirector &cstream = *fDebugStreamer;
352 *fDebugStreamer<<"match"<<
353 "x_tpc="<<x_tpc<<
354 "y_tpc="<<y_tpc<<
355 "z_tpc="<<z_tpc<<
356 "snp_tpc="<<snp_tpc<<
357 "tgl_tpc="<<tgl_tpc<<
358 "signed1Pt_tpc="<<signed1Pt_tpc<<
359 "x_trd="<<x_trd<<
360 "y_trd="<<y_trd<<
361 "z_trd="<<z_trd<<
362 "snp_trd="<<snp_trd<<
363 "tgl_trd="<<tgl_trd<<
364 "signed1Pt_trd="<<signed1Pt_trd<<
365 "\n";
366 }
367
368 if (TMath::Abs(y_tpc-y_trd) > fMaxY) return kFALSE;
369 if (TMath::Abs(z_tpc-z_trd) > fMaxZ) return kFALSE;
370 if (TMath::Abs(snp_tpc-snp_trd) > fMaxSnp) return kFALSE;
371 if (TMath::Abs(tgl_tpc-tgl_trd) > fMaxTgl) return kFALSE;
372 if (TMath::Abs(signed1Pt_tpc-signed1Pt_trd) > fMaxSigned1Pt) return kFALSE;
ec6160d5 373
374return kTRUE;
375}
376
377//_____________________________________________________________________________
378Bool_t AliHLTGlobalTrackMerger::SmoothTracks( const Double_t T1[], const Double_t C1[], Double_t Chi21, Int_t NDF1,
379 const Double_t T2[], const Double_t C2[], Double_t Chi22, Int_t NDF2,
380 Double_t T [], Double_t C [], Double_t &Chi2, Int_t &NDF,
381 Int_t N )
382{
383 //* Smooth two tracks with parameter vectors of size N
384 //*
385 //* Input:
386 //*
387 //* T1[N], T2[N] - tracks
388 //* C1[N*(N+1)/2], C2[N*(N+1)/2] - covariance matrices in low-diagonal form:
389 //* C = { c00,
390 //* c10, c11,
391 //* c20, c21, c22,
392 //* ... };
393 //* Chi2{1,2}, NDF{1,2} - \Chi^2 and "Number of Degrees of Freedom" values for both tracks
394 //* Output:
395 //*
396 //* T[N], C[N] ( can be aqual to {T1,C1}, or {T2,C2} )
397 //* Chi2, NDF
398 //*
399 //* returns error flag (0 means OK, 1 not OK )
400
401 Int_t M = N*(N+1)/2;
402
403 Double_t A[M];
404 Double_t K[N*N];
405
406 for(Int_t k=0; k<M; k++) A[k] = C1[k] + C2[k];
407 Bool_t err = InvertS(A,N);
408 if( err ) return 1;
409
410 Chi2 = Chi21 + Chi22;
411 NDF = NDF1 + NDF2;
412
413 MultSSQ( C1, A, K, N);
414 Double_t r[N];
415 for( Int_t k=0; k<N;k++) r[k] = T1[k] - T2[k];
416 for( Int_t k=0; k<N;k++ )
417 for( Int_t l=0;l<N;l++) T[k] = T1[k] - K[k*N+l]*r[l];
418
419 for( Int_t ind=0,i=0; i<N; i++ ){
420 for( Int_t j=0; j<i; j++ ) Chi2+= 2*r[i]*r[j]*A[ind++];
421 Chi2+= r[i]*r[i]*A[ind++];
422 }
423 NDF+=N;
424
425 for( Int_t l=0; l<N; l++ ) K[ (N+1)*l ] -= 1;
426
427 for( Int_t ind = 0, l=0; l<N; ++l ){
428 for( Int_t j=0; j<=l; ++j, ind++ ){
429 A[ind] = 0;
430 for( Int_t k=0; k<N; ++k ) A[ind] -= K[l*N+k] * C1[IndexS(j,k)];
431 }
432 }
433 for( Int_t l=0; l<N; l++ ) C[l] = A[l];
434 return 0;
435}
436
437//_____________________________________________________________________________
438void AliHLTGlobalTrackMerger::MultSSQ( const Double_t *A, const Double_t *B, Double_t *C, Int_t N )
439{
440 for( Int_t ind=0, i=0; i<N; ++i ){
441 for( Int_t j=0; j<N; ++j, ++ind ){
442 C[ind] = 0;
443 for( Int_t k=0; k<N; ++k ) C[ind] += A[IndexS(i,k)] * B[IndexS(k,j)];
444 }
445 }
446}
447
448//_____________________________________________________________________________
449Bool_t AliHLTGlobalTrackMerger::InvertS( Double_t A[], Int_t N )
450{
451 //* input: simmetric > 0 NxN matrix A = {a11,a21,a22,a31..a33,..}
452 //* output: inverse A, in case of problems fill zero and return 1
453 //*
454 //* A->low triangular Anew : A = Anew x Anew^T
455 //* method:
456 //* for(j=1,N) for(i=j,N) Aij=(Aii-sum_{k=1}^{j-1}Aik*Ajk )/Ajj
457 //*
458
459 Bool_t ret = 0;
460
461 const Double_t ZERO = 1.E-20;
462
463 {
464 Double_t *j1 = A, *jj = A;
465 for( Int_t j=1; j<=N; j1+=j++, jj+=j ){
466 Double_t *ik = j1, x = 0;
467 while( ik!=jj ){
468 x -= (*ik) * (*ik);
469 ik++;
470 }
471 x += *ik;
472 if( x > ZERO ){
473 x = sqrt(x);
474 *ik = x;
475 ik++;
476 x = 1 / x;
477 for( Int_t step=1; step<=N-j; ik+=++step ){ // ik==Ai1
478 Double_t sum = 0;
bb7774d6 479 for( Double_t *jk=j1; jk!=jj; sum += (*(jk++)) * (*(ik++)) ) {}
ec6160d5 480 *ik = (*ik - sum) * x; // ik == Aij
481 }
482 }else{
483 Double_t *ji=jj;
484 for( Int_t i=j; i<N; i++ ) *(ji+=i) = 0.;
485 ret = 1;
486 }
487 }
488 }
489
490 //* A -> Ainv
491 //* method :
492 //* for(i=1,N){
493 //* Aii = 1/Aii;
494 //* for(j=1,i-1) Aij=-(sum_{k=j}^{i-1} Aik * Akj) / Aii ;
495 //* }
496
497 {
498 Double_t *ii=A,*ij=A;
499 for( Int_t i = 1; i<=N; ij=ii+1, ii+=++i ){
500 if( *ii > ZERO ){
501 Double_t x = -(*ii = 1./ *ii);
502 {
503 Double_t *jj = A;
504 for( Int_t j=1; j<i; jj+=++j, ij++ ){
505 Double_t *ik = ij, *kj = jj, sum = 0.;
506 for( Int_t k=j; ik!=ii; kj+=k++, ik++ ){
507 sum += *ik * *kj;
508 }
509 *kj = sum * x;
510 }
511 }
512 }else{
513 for( Double_t *ik = ij; ik!=ii+1; ik++ ){
514 *ik = 0.;
515 }
516 ret = 1;
517 }
518 }
519 }
520
521 //* A -> A^T x A
522 //* method:
523 //* Aij = sum_{k=i}^N Aki * Akj
524
525 {
526 Double_t *ii=A, *ij=A;
527 for( Int_t i=1; i<=N; ii+=++i ){
528 do{
529 Double_t *ki = ii, *kj = ij, sum = 0.;
530 for( Int_t k=i; k<=N; ki+=k, kj+=k++ ) sum += (*ki) * (*kj);
531 *ij = sum;
532 }while( (ij++)!=ii );
533 }
534 }
535 return ret;
536}
537
538//_____________________________________________________________________________
539void AliHLTGlobalTrackMerger::PropagateTracksToDCA(AliESDEvent *esdEvent)
540{
541 // try to propagate all tracks to DCA to primary vertex
542 if(!esdEvent) return;
543
544 const Double_t kBz = esdEvent->GetMagneticField();
545 const Double_t kSmallRadius = 2.8; // [cm] something less than the beam pipe radius
546 const Double_t kMaxStep = 10.0; // [cm] track propagation step
547 Bool_t isOK = kFALSE;
548
549 Int_t nTracks = esdEvent->GetNumberOfTracks();
550 for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
551 AliESDtrack *track = esdEvent->GetTrack(iTrack);
552 if(!track) continue;
553
ec6160d5 554 // propagate to small radius (material budget included)
555 isOK = AliTracker::PropagateTrackTo(track,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);
556
ec6160d5 557 // relate tracks to DCA to primary vertex
558 if(isOK)
559 {
560 track->RelateToVertex(fVertex, kBz, kVeryBig);
a9f060e5 561 HLTInfo("1-------: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",track->GetAlpha(),track->GetX(),track->GetY(),track->GetZ(),track->GetSnp(),track->GetTgl(),track->GetSigned1Pt());
562 }
563
564 //
565 // the same procedure must be repeated for TPCinner (TPC only) tracks
566 //
567 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
568 if(!tpcTrack) continue;
569
570 // propagate tpcTracks to small radius (material budget included)
571 isOK = AliTracker::PropagateTrackTo(tpcTrack,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);
572
573 // relate tracks to DCA to primary vertex
574 if(isOK)
575 {
576 Double_t par[2], cov[3];
577 tpcTrack->PropagateToDCA(fVertex, kBz, kVeryBig,par,cov);
578 HLTInfo("2-------: alpha %f, x %f, y %f, z %f, snp %f, tgl %f, 1pt %f",tpcTrack->GetAlpha(),tpcTrack->GetX(),tpcTrack->GetY(),tpcTrack->GetZ(),tpcTrack->GetSnp(),tpcTrack->GetTgl(),tpcTrack->GetSigned1Pt());
ec6160d5 579 }
580 }
581}