]>
Commit | Line | Data |
---|---|---|
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 | ||
18ada816 | 25 | //#include "AliTPCReconstructor.h" |
ec6160d5 | 26 | |
ec6160d5 | 27 | #include "AliESDEvent.h" |
28 | #include "AliESDVertex.h" | |
29 | #include "AliTracker.h" | |
a9f060e5 | 30 | #include "TTreeStream.h" |
ec6160d5 | 31 | |
32 | #include <TClonesArray.h> | |
33 | ||
34 | #include "AliHLTGlobalTrackMerger.h" | |
35 | ||
36 | #if __GNUC__>= 3 | |
37 | using namespace std; | |
38 | #endif | |
39 | ||
40 | ClassImp(AliHLTGlobalTrackMerger) | |
41 | ||
42 | AliHLTGlobalTrackMerger::AliHLTGlobalTrackMerger() | |
43 | : | |
44 | fMaxY(0.0), | |
45 | fMaxZ(0.0), | |
46 | fMaxSnp(0.0), | |
47 | fMaxTgl(0.0), | |
48 | fMaxSigned1Pt(0.0), | |
a9f060e5 | 49 | fVertex(0), |
50 | fDebugStreamer(0) | |
ec6160d5 | 51 | { |
52 | //Default constructor | |
53 | ||
54 | // standard vertex settings at the moment | |
55 | // V(0.,0.,0.), sigmaVx=sigmaVy=5.e-3 [cm], sigmaVz=5.3 [cm] | |
56 | fVertex = new AliESDVertex; | |
a9f060e5 | 57 | |
18ada816 | 58 | // if (AliTPCReconstructor::StreamLevel()>0) { |
59 | // fDebugStreamer = new TTreeSRedirector("debugGlobalMerger.root"); | |
60 | // } | |
a9f060e5 | 61 | |
ec6160d5 | 62 | } |
63 | ||
64 | //_____________________________________________________________________________ | |
65 | AliHLTGlobalTrackMerger::~AliHLTGlobalTrackMerger() | |
66 | { | |
67 | //Destructor | |
68 | if(fVertex) delete fVertex; fVertex =0; | |
a9f060e5 | 69 | if(fDebugStreamer) delete fDebugStreamer; fDebugStreamer =0; |
ec6160d5 | 70 | } |
71 | ||
ec6160d5 | 72 | //_____________________________________________________________________________ |
73 | Bool_t AliHLTGlobalTrackMerger::Merge(AliESDEvent* esdEvent) | |
74 | { | |
75 | // merge TPC and TRD tracks | |
76 | // 1. propagate TPC track to the radius between TPC and TRD | |
77 | // 2. propagate TRD track to the same radius between TPC and TRD | |
78 | // 3. matches TPC and TRD tracks at the radius | |
4b6d0cc6 | 79 | // 4. propagate matched TRD track to the merging radius (first measured TPC point - x coordinate) |
ec6160d5 | 80 | // 5. merge TPC and TRD track parameters at the merging radius |
81 | // 6. create AliESDtrack from merged tracks | |
82 | // 7. add AliESDtrack to AliESDEvent | |
83 | ||
84 | if(!esdEvent) return kFALSE; | |
85 | ||
86 | const Double_t kMaxStep = 10.0; // [cm] track propagation step | |
87 | const Double_t kMatchRadius = 285.0; // [cm] matching at radius between TPC and TRD | |
88 | Double_t kMergeRadius = 0.0; | |
89 | Bool_t isOk = kFALSE; | |
90 | Bool_t isMatched = kFALSE; | |
a9f060e5 | 91 | AliESDtrack *track=0; |
4b6d0cc6 | 92 | AliExternalTrackParam *extTPCTrack = 0; |
ec6160d5 | 93 | |
94 | Int_t nTracks = esdEvent->GetNumberOfTracks(); | |
ec6160d5 | 95 | HLTWarning("nTracks %d",nTracks); |
96 | ||
a9f060e5 | 97 | Int_t nTracksTPC =0; |
ec6160d5 | 98 | for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) |
99 | { | |
a9f060e5 | 100 | track = esdEvent->GetTrack(iTrack); |
ec6160d5 | 101 | if(!track) continue; |
102 | ||
103 | // TPC tracks | |
104 | if((track->GetStatus()&AliESDtrack::kTPCin)==0) continue; | |
105 | AliESDtrack *tpcTrack = track; | |
a9f060e5 | 106 | nTracksTPC++; |
ec6160d5 | 107 | |
4b6d0cc6 | 108 | // create external tpc track param (needed to propagate to matching radius) |
109 | if ((extTPCTrack = new AliExternalTrackParam(*tpcTrack)) == 0) continue; | |
110 | ||
a9f060e5 | 111 | kMergeRadius = tpcTrack->GetTPCPoints(0); // [cm] merging at first measured TPC point |
ec6160d5 | 112 | |
a9f060e5 | 113 | HLTInfo("-------------------------------------------------------------------------------------"); |
4b6d0cc6 | 114 | //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()); |
115 | 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 | 116 | |
117 | // propagate tracks to the matching radius | |
4b6d0cc6 | 118 | //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE); |
119 | isOk = AliTracker::PropagateTrackTo(extTPCTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE); | |
ec6160d5 | 120 | if(!isOk) continue; |
121 | ||
4b6d0cc6 | 122 | 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 | 123 | |
a9f060e5 | 124 | Int_t nTracksTRD =0; |
ec6160d5 | 125 | for(Int_t jTrack = 0; jTrack<nTracks; ++jTrack) |
126 | { | |
a9f060e5 | 127 | track = esdEvent->GetTrack(jTrack); |
ec6160d5 | 128 | if(!track) continue; |
129 | ||
130 | // TRD tracks | |
131 | if((track->GetStatus()&AliESDtrack::kTRDin)==0) continue; | |
132 | AliESDtrack *trdTrack = track; | |
a9f060e5 | 133 | nTracksTRD++; |
134 | ||
4b6d0cc6 | 135 | 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 | 136 | |
137 | isOk = AliTracker::PropagateTrackTo(trdTrack,kMatchRadius,trdTrack->GetMass(),kMaxStep,kFALSE); | |
138 | if(!isOk) continue; | |
139 | ||
4b6d0cc6 | 140 | 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 | 141 | |
142 | // match TPC and TRD tracks | |
4b6d0cc6 | 143 | //isMatched = MatchTracks(tpcTrack,trdTrack); |
144 | isMatched = MatchTracks(extTPCTrack,trdTrack); | |
ec6160d5 | 145 | if(!isMatched) continue; |
146 | ||
4b6d0cc6 | 147 | //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()); |
148 | 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()); | |
149 | 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 | 150 | |
4b6d0cc6 | 151 | // propagate TRD tracks to the merging radius |
152 | //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMergeRadius,tpcTrack->GetMass(),kMaxStep,kFALSE); | |
153 | //if(!isOk) continue; | |
ec6160d5 | 154 | isOk = AliTracker::PropagateTrackTo(trdTrack,kMergeRadius,trdTrack->GetMass(),kMaxStep,kFALSE); |
155 | if(!isOk) continue; | |
156 | ||
4b6d0cc6 | 157 | 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()); |
158 | 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 | 159 | |
160 | // merge TPC and TRD tracks | |
161 | // create AliESDtrack and add it to AliESDevent | |
162 | Bool_t isMerged = MergeTracks(tpcTrack,trdTrack,esdEvent); | |
a9f060e5 | 163 | if(!isMerged) |
164 | HLTInfo("No merged tracks"); | |
ec6160d5 | 165 | } |
a9f060e5 | 166 | HLTInfo("nTracksTRD %d",nTracksTRD); |
4b6d0cc6 | 167 | |
168 | // delete external TPC track | |
169 | if(extTPCTrack) delete extTPCTrack; extTPCTrack=0; | |
ec6160d5 | 170 | } |
a9f060e5 | 171 | HLTInfo("nTracksTPC %d",nTracksTPC); |
ec6160d5 | 172 | |
173 | return kTRUE; | |
174 | } | |
175 | ||
176 | //_____________________________________________________________________________ | |
177 | Bool_t AliHLTGlobalTrackMerger::MergeTracks(AliESDtrack *tpcTrack, AliESDtrack* trdTrack, AliESDEvent *esdEvent) | |
178 | { | |
179 | // merge TPC and TRD track parameters | |
180 | // create new AliESDtrack with TPC+TRD merged track parameters | |
181 | // add AliESDtrack to AliESDEvent | |
182 | ||
183 | if(!tpcTrack) return kFALSE; | |
184 | if(!trdTrack) return kFALSE; | |
185 | ||
ec6160d5 | 186 | Double_t trackParam[5], trackCovar[15]; |
187 | Double_t trackChi2; | |
188 | Int_t trackNDF; | |
189 | ||
190 | // calculate merged track parameters | |
191 | Bool_t isNotOK = SmoothTracks(tpcTrack->GetParameter(), tpcTrack->GetCovariance(), tpcTrack->GetTPCchi2(), 5, | |
192 | trdTrack->GetParameter(), trdTrack->GetCovariance(), trdTrack->GetTRDchi2(), 5, | |
193 | trackParam, trackCovar, trackChi2, trackNDF,5); | |
194 | ||
195 | if(isNotOK) | |
196 | return kFALSE; | |
197 | ||
198 | // | |
199 | // create AliESDtrack | |
200 | // merged TPC+TRD information | |
201 | // | |
202 | ||
203 | AliESDtrack track; | |
204 | //track.UpdateTrackParams(tpcTrack, AliESDtrack::kTPCrefit); | |
205 | track.SetStatus(AliESDtrack::kGlobalMerge); | |
206 | track.SetLabel(tpcTrack->GetLabel()); | |
207 | track.Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),trackParam,trackCovar); | |
208 | track.SetGlobalChi2(trackChi2); | |
209 | ||
210 | //track.SetTPCLabel(tpcTrack->GetLabel()); | |
211 | Double32_t tpcPID[AliPID::kSPECIES]; | |
212 | tpcTrack->GetTPCpid(tpcPID); | |
213 | track.SetTPCpid(tpcPID); | |
214 | //fTPCncls=t->GetNumberOfClusters(); // no cluster on HLT | |
215 | //fTPCchi2=t->GetChi2(); | |
216 | ||
217 | //track.SetTRDLabel(trdTrack->GetLabel()); | |
218 | Double32_t trdPID[AliPID::kSPECIES]; | |
219 | trdTrack->GetTRDpid(trdPID); | |
220 | track.SetTRDpid(trdPID); | |
221 | //fTRDchi2 = t->GetChi2(); | |
222 | //fTRDncls = t->GetNumberOfClusters(); | |
223 | //for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i); | |
224 | ||
225 | // add track to AliESDEvent | |
226 | esdEvent->AddTrack(&track); | |
227 | ||
228 | return kTRUE; | |
229 | } | |
230 | ||
231 | //_____________________________________________________________________________ | |
232 | void AliHLTGlobalTrackMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxsnp, Double_t maxtgl, Double_t signed1Pt) | |
233 | { | |
234 | //set parameters for merger | |
235 | fMaxY = maxy; | |
236 | fMaxZ = maxz; | |
237 | fMaxSnp = maxsnp; | |
238 | fMaxTgl = maxtgl; | |
239 | fMaxSigned1Pt = signed1Pt; | |
240 | } | |
241 | ||
242 | //_____________________________________________________________________________ | |
4b6d0cc6 | 243 | Bool_t AliHLTGlobalTrackMerger::MatchTracks(AliExternalTrackParam *trackTPC, AliESDtrack *trackTRD) |
ec6160d5 | 244 | { |
245 | // match TPC and TRD tracks | |
246 | // return kTRUE in case of matching | |
247 | ||
248 | if(!trackTPC) return kFALSE; | |
249 | if(!trackTRD) return kFALSE; | |
250 | ||
a9f060e5 | 251 | Double_t x_tpc=trackTPC->GetX(); |
252 | Double_t y_tpc=trackTPC->GetY(); | |
253 | Double_t z_tpc=trackTPC->GetZ(); | |
254 | Double_t snp_tpc=trackTPC->GetSnp(); | |
255 | Double_t tgl_tpc=trackTPC->GetTgl(); | |
256 | Double_t signed1Pt_tpc=trackTPC->GetSigned1Pt(); | |
257 | ||
258 | Double_t x_trd=trackTRD->GetX(); | |
259 | Double_t y_trd=trackTRD->GetY(); | |
260 | Double_t z_trd=trackTRD->GetZ(); | |
261 | Double_t snp_trd=trackTRD->GetSnp(); | |
262 | Double_t tgl_trd=trackTRD->GetTgl(); | |
263 | Double_t signed1Pt_trd=trackTRD->GetSigned1Pt(); | |
264 | ||
265 | // debug stream | |
18ada816 | 266 | // if (AliTPCReconstructor::StreamLevel()>0) { |
267 | // //TTreeSRedirector &cstream = *fDebugStreamer; | |
268 | // *fDebugStreamer<<"match"<< | |
269 | // "x_tpc="<<x_tpc<< | |
270 | // "y_tpc="<<y_tpc<< | |
271 | // "z_tpc="<<z_tpc<< | |
272 | // "snp_tpc="<<snp_tpc<< | |
273 | // "tgl_tpc="<<tgl_tpc<< | |
274 | // "signed1Pt_tpc="<<signed1Pt_tpc<< | |
275 | // "x_trd="<<x_trd<< | |
276 | // "y_trd="<<y_trd<< | |
277 | // "z_trd="<<z_trd<< | |
278 | // "snp_trd="<<snp_trd<< | |
279 | // "tgl_trd="<<tgl_trd<< | |
280 | // "signed1Pt_trd="<<signed1Pt_trd<< | |
281 | // "\n"; | |
282 | // } | |
a9f060e5 | 283 | |
8b7a0f3c | 284 | if (TMath::Abs(x_tpc-x_trd) > 0) {/* get rid of warning*/;} |
a9f060e5 | 285 | if (TMath::Abs(y_tpc-y_trd) > fMaxY) return kFALSE; |
286 | if (TMath::Abs(z_tpc-z_trd) > fMaxZ) return kFALSE; | |
287 | if (TMath::Abs(snp_tpc-snp_trd) > fMaxSnp) return kFALSE; | |
288 | if (TMath::Abs(tgl_tpc-tgl_trd) > fMaxTgl) return kFALSE; | |
289 | if (TMath::Abs(signed1Pt_tpc-signed1Pt_trd) > fMaxSigned1Pt) return kFALSE; | |
ec6160d5 | 290 | |
291 | return kTRUE; | |
292 | } | |
293 | ||
294 | //_____________________________________________________________________________ | |
295 | Bool_t AliHLTGlobalTrackMerger::SmoothTracks( const Double_t T1[], const Double_t C1[], Double_t Chi21, Int_t NDF1, | |
296 | const Double_t T2[], const Double_t C2[], Double_t Chi22, Int_t NDF2, | |
297 | Double_t T [], Double_t C [], Double_t &Chi2, Int_t &NDF, | |
298 | Int_t N ) | |
299 | { | |
300 | //* Smooth two tracks with parameter vectors of size N | |
301 | //* | |
302 | //* Input: | |
303 | //* | |
304 | //* T1[N], T2[N] - tracks | |
305 | //* C1[N*(N+1)/2], C2[N*(N+1)/2] - covariance matrices in low-diagonal form: | |
306 | //* C = { c00, | |
307 | //* c10, c11, | |
308 | //* c20, c21, c22, | |
309 | //* ... }; | |
310 | //* Chi2{1,2}, NDF{1,2} - \Chi^2 and "Number of Degrees of Freedom" values for both tracks | |
311 | //* Output: | |
312 | //* | |
313 | //* T[N], C[N] ( can be aqual to {T1,C1}, or {T2,C2} ) | |
314 | //* Chi2, NDF | |
315 | //* | |
316 | //* returns error flag (0 means OK, 1 not OK ) | |
317 | ||
318 | Int_t M = N*(N+1)/2; | |
319 | ||
320 | Double_t A[M]; | |
321 | Double_t K[N*N]; | |
322 | ||
323 | for(Int_t k=0; k<M; k++) A[k] = C1[k] + C2[k]; | |
324 | Bool_t err = InvertS(A,N); | |
325 | if( err ) return 1; | |
326 | ||
327 | Chi2 = Chi21 + Chi22; | |
328 | NDF = NDF1 + NDF2; | |
329 | ||
330 | MultSSQ( C1, A, K, N); | |
331 | Double_t r[N]; | |
332 | for( Int_t k=0; k<N;k++) r[k] = T1[k] - T2[k]; | |
333 | for( Int_t k=0; k<N;k++ ) | |
334 | for( Int_t l=0;l<N;l++) T[k] = T1[k] - K[k*N+l]*r[l]; | |
335 | ||
336 | for( Int_t ind=0,i=0; i<N; i++ ){ | |
337 | for( Int_t j=0; j<i; j++ ) Chi2+= 2*r[i]*r[j]*A[ind++]; | |
338 | Chi2+= r[i]*r[i]*A[ind++]; | |
339 | } | |
340 | NDF+=N; | |
341 | ||
342 | for( Int_t l=0; l<N; l++ ) K[ (N+1)*l ] -= 1; | |
343 | ||
344 | for( Int_t ind = 0, l=0; l<N; ++l ){ | |
345 | for( Int_t j=0; j<=l; ++j, ind++ ){ | |
346 | A[ind] = 0; | |
347 | for( Int_t k=0; k<N; ++k ) A[ind] -= K[l*N+k] * C1[IndexS(j,k)]; | |
348 | } | |
349 | } | |
350 | for( Int_t l=0; l<N; l++ ) C[l] = A[l]; | |
351 | return 0; | |
352 | } | |
353 | ||
354 | //_____________________________________________________________________________ | |
355 | void AliHLTGlobalTrackMerger::MultSSQ( const Double_t *A, const Double_t *B, Double_t *C, Int_t N ) | |
356 | { | |
357 | for( Int_t ind=0, i=0; i<N; ++i ){ | |
358 | for( Int_t j=0; j<N; ++j, ++ind ){ | |
359 | C[ind] = 0; | |
360 | for( Int_t k=0; k<N; ++k ) C[ind] += A[IndexS(i,k)] * B[IndexS(k,j)]; | |
361 | } | |
362 | } | |
363 | } | |
364 | ||
365 | //_____________________________________________________________________________ | |
366 | Bool_t AliHLTGlobalTrackMerger::InvertS( Double_t A[], Int_t N ) | |
367 | { | |
368 | //* input: simmetric > 0 NxN matrix A = {a11,a21,a22,a31..a33,..} | |
369 | //* output: inverse A, in case of problems fill zero and return 1 | |
370 | //* | |
371 | //* A->low triangular Anew : A = Anew x Anew^T | |
372 | //* method: | |
373 | //* for(j=1,N) for(i=j,N) Aij=(Aii-sum_{k=1}^{j-1}Aik*Ajk )/Ajj | |
374 | //* | |
375 | ||
376 | Bool_t ret = 0; | |
377 | ||
378 | const Double_t ZERO = 1.E-20; | |
379 | ||
380 | { | |
381 | Double_t *j1 = A, *jj = A; | |
382 | for( Int_t j=1; j<=N; j1+=j++, jj+=j ){ | |
383 | Double_t *ik = j1, x = 0; | |
384 | while( ik!=jj ){ | |
385 | x -= (*ik) * (*ik); | |
386 | ik++; | |
387 | } | |
388 | x += *ik; | |
389 | if( x > ZERO ){ | |
390 | x = sqrt(x); | |
391 | *ik = x; | |
392 | ik++; | |
393 | x = 1 / x; | |
394 | for( Int_t step=1; step<=N-j; ik+=++step ){ // ik==Ai1 | |
395 | Double_t sum = 0; | |
bb7774d6 | 396 | for( Double_t *jk=j1; jk!=jj; sum += (*(jk++)) * (*(ik++)) ) {} |
ec6160d5 | 397 | *ik = (*ik - sum) * x; // ik == Aij |
398 | } | |
399 | }else{ | |
400 | Double_t *ji=jj; | |
401 | for( Int_t i=j; i<N; i++ ) *(ji+=i) = 0.; | |
402 | ret = 1; | |
403 | } | |
404 | } | |
405 | } | |
406 | ||
407 | //* A -> Ainv | |
408 | //* method : | |
409 | //* for(i=1,N){ | |
410 | //* Aii = 1/Aii; | |
411 | //* for(j=1,i-1) Aij=-(sum_{k=j}^{i-1} Aik * Akj) / Aii ; | |
412 | //* } | |
413 | ||
414 | { | |
415 | Double_t *ii=A,*ij=A; | |
416 | for( Int_t i = 1; i<=N; ij=ii+1, ii+=++i ){ | |
417 | if( *ii > ZERO ){ | |
418 | Double_t x = -(*ii = 1./ *ii); | |
419 | { | |
420 | Double_t *jj = A; | |
421 | for( Int_t j=1; j<i; jj+=++j, ij++ ){ | |
422 | Double_t *ik = ij, *kj = jj, sum = 0.; | |
423 | for( Int_t k=j; ik!=ii; kj+=k++, ik++ ){ | |
424 | sum += *ik * *kj; | |
425 | } | |
426 | *kj = sum * x; | |
427 | } | |
428 | } | |
429 | }else{ | |
430 | for( Double_t *ik = ij; ik!=ii+1; ik++ ){ | |
431 | *ik = 0.; | |
432 | } | |
433 | ret = 1; | |
434 | } | |
435 | } | |
436 | } | |
437 | ||
438 | //* A -> A^T x A | |
439 | //* method: | |
440 | //* Aij = sum_{k=i}^N Aki * Akj | |
441 | ||
442 | { | |
443 | Double_t *ii=A, *ij=A; | |
444 | for( Int_t i=1; i<=N; ii+=++i ){ | |
445 | do{ | |
446 | Double_t *ki = ii, *kj = ij, sum = 0.; | |
447 | for( Int_t k=i; k<=N; ki+=k, kj+=k++ ) sum += (*ki) * (*kj); | |
448 | *ij = sum; | |
449 | }while( (ij++)!=ii ); | |
450 | } | |
451 | } | |
452 | return ret; | |
453 | } | |
454 | ||
455 | //_____________________________________________________________________________ | |
456 | void AliHLTGlobalTrackMerger::PropagateTracksToDCA(AliESDEvent *esdEvent) | |
457 | { | |
458 | // try to propagate all tracks to DCA to primary vertex | |
459 | if(!esdEvent) return; | |
460 | ||
461 | const Double_t kBz = esdEvent->GetMagneticField(); | |
462 | const Double_t kSmallRadius = 2.8; // [cm] something less than the beam pipe radius | |
463 | const Double_t kMaxStep = 10.0; // [cm] track propagation step | |
464 | Bool_t isOK = kFALSE; | |
465 | ||
466 | Int_t nTracks = esdEvent->GetNumberOfTracks(); | |
467 | for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) { | |
468 | AliESDtrack *track = esdEvent->GetTrack(iTrack); | |
469 | if(!track) continue; | |
470 | ||
ec6160d5 | 471 | // propagate to small radius (material budget included) |
472 | isOK = AliTracker::PropagateTrackTo(track,kSmallRadius,track->GetMass(),kMaxStep,kFALSE); | |
473 | ||
ec6160d5 | 474 | // relate tracks to DCA to primary vertex |
475 | if(isOK) | |
476 | { | |
477 | track->RelateToVertex(fVertex, kBz, kVeryBig); | |
a9f060e5 | 478 | 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()); |
479 | } | |
480 | ||
481 | // | |
482 | // the same procedure must be repeated for TPCinner (TPC only) tracks | |
483 | // | |
484 | AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam(); | |
485 | if(!tpcTrack) continue; | |
486 | ||
487 | // propagate tpcTracks to small radius (material budget included) | |
488 | isOK = AliTracker::PropagateTrackTo(tpcTrack,kSmallRadius,track->GetMass(),kMaxStep,kFALSE); | |
489 | ||
490 | // relate tracks to DCA to primary vertex | |
491 | if(isOK) | |
492 | { | |
493 | Double_t par[2], cov[3]; | |
494 | tpcTrack->PropagateToDCA(fVertex, kBz, kVeryBig,par,cov); | |
495 | 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 | 496 | } |
497 | } | |
498 | } |