3 /**************************************************************************
4 * This file is property of and copyright by the ALICE HLT Project *
5 * ALICE Experiment at CERN, All rights reserved. *
7 * Primary Authors: Jacek Otwinowski (Jacek.Otwinowski@gsi.de) *
8 * for The ALICE HLT Project. *
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 **************************************************************************/
19 // @file AliHLTGlobalTrackMerger.cxx
20 // @author Jacek Otwinowski
22 // @brief The HLT global merger base class
25 //#include "AliTPCReconstructor.h"
27 #include "AliESDEvent.h"
28 #include "AliESDVertex.h"
29 #include "AliESDtrack.h"
30 #include "AliTracker.h"
31 #include "TTreeStream.h"
33 #include <TClonesArray.h>
35 #include "AliHLTGlobalTrackMerger.h"
41 ClassImp(AliHLTGlobalTrackMerger)
43 AliHLTGlobalTrackMerger::AliHLTGlobalTrackMerger()
55 // standard vertex settings at the moment
56 // V(0.,0.,0.), sigmaVx=sigmaVy=5.e-3 [cm], sigmaVz=5.3 [cm]
57 fVertex = new AliESDVertex;
59 // if (AliTPCReconstructor::StreamLevel()>0) {
60 // fDebugStreamer = new TTreeSRedirector("debugGlobalMerger.root");
65 //_____________________________________________________________________________
66 AliHLTGlobalTrackMerger::~AliHLTGlobalTrackMerger()
69 if(fVertex) delete fVertex; fVertex =0;
70 if(fDebugStreamer) delete fDebugStreamer; fDebugStreamer =0;
73 //_____________________________________________________________________________
74 Bool_t AliHLTGlobalTrackMerger::Merge(AliESDEvent* esdEvent)
76 // merge TPC and TRD tracks
77 // 1. propagate TPC track to the radius between TPC and TRD
78 // 2. propagate TRD track to the same radius between TPC and TRD
79 // 3. matches TPC and TRD tracks at the radius
80 // 4. propagate matched TRD track to the merging radius (first measured TPC point - x coordinate)
81 // 5. merge TPC and TRD track parameters at the merging radius
82 // 6. create AliESDtrack from merged tracks
83 // 7. add AliESDtrack to AliESDEvent
85 if(!esdEvent) return kFALSE;
87 const Double_t kMaxStep = 10.0; // [cm] track propagation step
88 const Double_t kMatchRadius = 285.0; // [cm] matching at radius between TPC and TRD
89 Double_t kMergeRadius = 0.0;
91 Bool_t isMatched = kFALSE;
93 AliExternalTrackParam *extTPCTrack = 0;
95 Int_t nTracks = esdEvent->GetNumberOfTracks();
96 HLTWarning("nTracks %d",nTracks);
99 for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack)
101 track = esdEvent->GetTrack(iTrack);
105 if((track->GetStatus()&AliESDtrack::kTPCin)==0) continue;
106 AliESDtrack *tpcTrack = track;
109 // create external tpc track param (needed to propagate to matching radius)
110 if ((extTPCTrack = new AliExternalTrackParam(*tpcTrack)) == 0) continue;
112 kMergeRadius = tpcTrack->GetTPCPoints(0); // [cm] merging at first measured TPC point
114 HLTInfo("-------------------------------------------------------------------------------------");
115 //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());
116 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());
118 // propagate tracks to the matching radius
119 //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
120 isOk = AliTracker::PropagateTrackTo(extTPCTrack,kMatchRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
123 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());
126 for(Int_t jTrack = 0; jTrack<nTracks; ++jTrack)
128 track = esdEvent->GetTrack(jTrack);
132 if((track->GetStatus()&AliESDtrack::kTRDin)==0) continue;
133 AliESDtrack *trdTrack = track;
136 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());
138 isOk = AliTracker::PropagateTrackTo(trdTrack,kMatchRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
141 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());
143 // match TPC and TRD tracks
144 //isMatched = MatchTracks(tpcTrack,trdTrack);
145 isMatched = MatchTracks(extTPCTrack,trdTrack);
146 if(!isMatched) continue;
148 //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());
149 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());
150 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());
152 // propagate TRD tracks to the merging radius
153 //isOk = AliTracker::PropagateTrackTo(tpcTrack,kMergeRadius,tpcTrack->GetMass(),kMaxStep,kFALSE);
154 //if(!isOk) continue;
155 isOk = AliTracker::PropagateTrackTo(trdTrack,kMergeRadius,trdTrack->GetMass(),kMaxStep,kFALSE);
158 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());
159 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());
161 // merge TPC and TRD tracks
162 // create AliESDtrack and add it to AliESDevent
163 Bool_t isMerged = MergeTracks(tpcTrack,trdTrack,esdEvent);
165 HLTInfo("No merged tracks");
167 HLTInfo("nTracksTRD %d",nTracksTRD);
169 // delete external TPC track
170 if(extTPCTrack) delete extTPCTrack; extTPCTrack=0;
172 HLTInfo("nTracksTPC %d",nTracksTPC);
177 //_____________________________________________________________________________
178 Bool_t AliHLTGlobalTrackMerger::MergeTracks(AliESDtrack *tpcTrack, AliESDtrack* trdTrack, AliESDEvent *esdEvent)
180 // merge TPC and TRD track parameters
181 // create new AliESDtrack with TPC+TRD merged track parameters
182 // add AliESDtrack to AliESDEvent
184 if(!tpcTrack) return kFALSE;
185 if(!trdTrack) return kFALSE;
187 Double_t trackParam[5], trackCovar[15];
191 // calculate merged track parameters
192 Bool_t isNotOK = SmoothTracks(tpcTrack->GetParameter(), tpcTrack->GetCovariance(), tpcTrack->GetTPCchi2(), 5,
193 trdTrack->GetParameter(), trdTrack->GetCovariance(), trdTrack->GetTRDchi2(), 5,
194 trackParam, trackCovar, trackChi2, trackNDF,5);
200 // create AliESDtrack
201 // merged TPC+TRD information
205 //track.UpdateTrackParams(tpcTrack, AliESDtrack::kTPCrefit);
206 track.SetStatus(AliESDtrack::kGlobalMerge);
207 track.SetLabel(tpcTrack->GetLabel());
208 track.Set(tpcTrack->GetX(),tpcTrack->GetAlpha(),trackParam,trackCovar);
209 track.SetGlobalChi2(trackChi2);
211 //track.SetTPCLabel(tpcTrack->GetLabel());
212 Double32_t tpcPID[AliPID::kSPECIES];
213 tpcTrack->GetTPCpid(tpcPID);
214 track.SetTPCpid(tpcPID);
215 //fTPCncls=t->GetNumberOfClusters(); // no cluster on HLT
216 //fTPCchi2=t->GetChi2();
218 //track.SetTRDLabel(trdTrack->GetLabel());
219 Double32_t trdPID[AliPID::kSPECIES];
220 trdTrack->GetTRDpid(trdPID);
221 track.SetTRDpid(trdPID);
222 //fTRDchi2 = t->GetChi2();
223 //fTRDncls = t->GetNumberOfClusters();
224 //for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
226 // add track to AliESDEvent
227 esdEvent->AddTrack(&track);
232 //_____________________________________________________________________________
233 void AliHLTGlobalTrackMerger::SetParameter(Double_t maxy, Double_t maxz, Double_t maxsnp, Double_t maxtgl, Double_t signed1Pt)
235 //set parameters for merger
240 fMaxSigned1Pt = signed1Pt;
243 //_____________________________________________________________________________
244 Bool_t AliHLTGlobalTrackMerger::MatchTracks(AliExternalTrackParam *trackTPC, const AliESDtrack *trackTRD)
246 // match TPC and TRD tracks
247 // return kTRUE in case of matching
249 if(!trackTPC) return kFALSE;
250 if(!trackTRD) return kFALSE;
252 Double_t xTpc=trackTPC->GetX();
253 Double_t yTpc=trackTPC->GetY();
254 Double_t zTpc=trackTPC->GetZ();
255 Double_t snpTpc=trackTPC->GetSnp();
256 Double_t tglTpc=trackTPC->GetTgl();
257 Double_t signed1PtTpc=trackTPC->GetSigned1Pt();
259 Double_t xTrd=trackTRD->GetX();
260 Double_t yTrd=trackTRD->GetY();
261 Double_t zTrd=trackTRD->GetZ();
262 Double_t snpTrd=trackTRD->GetSnp();
263 Double_t tglTrd=trackTRD->GetTgl();
264 Double_t signed1PtTrd=trackTRD->GetSigned1Pt();
267 // if (AliTPCReconstructor::StreamLevel()>0) {
268 // //TTreeSRedirector &cstream = *fDebugStreamer;
269 // *fDebugStreamer<<"match"<<
273 // "snpTpc="<<snpTpc<<
274 // "tglTpc="<<tglTpc<<
275 // "signed1PtTpc="<<signed1PtTpc<<
279 // "snpTrd="<<snpTrd<<
280 // "tglTrd="<<tglTrd<<
281 // "signed1PtTrd="<<signed1PtTrd<<
285 if (TMath::Abs(xTpc-xTrd) > 0) {/* get rid of warning*/;}
286 if (TMath::Abs(yTpc-yTrd) > fMaxY) return kFALSE;
287 if (TMath::Abs(zTpc-zTrd) > fMaxZ) return kFALSE;
288 if (TMath::Abs(snpTpc-snpTrd) > fMaxSnp) return kFALSE;
289 if (TMath::Abs(tglTpc-tglTrd) > fMaxTgl) return kFALSE;
290 if (TMath::Abs(signed1PtTpc-signed1PtTrd) > fMaxSigned1Pt) return kFALSE;
295 //_____________________________________________________________________________
296 Bool_t AliHLTGlobalTrackMerger::SmoothTracks( const Double_t T1[], const Double_t C1[], Double_t Chi21, Int_t NDF1,
297 const Double_t T2[], const Double_t C2[], Double_t Chi22, Int_t NDF2,
298 Double_t T [], Double_t C [], Double_t &Chi2, Int_t &NDF,
301 //* Smooth two tracks with parameter vectors of size N
305 //* T1[N], T2[N] - tracks
306 //* C1[N*(N+1)/2], C2[N*(N+1)/2] - covariance matrices in low-diagonal form:
311 //* Chi2{1,2}, NDF{1,2} - \Chi^2 and "Number of Degrees of Freedom" values for both tracks
314 //* T[N], C[N] ( can be aqual to {T1,C1}, or {T2,C2} )
317 //* returns error flag (0 means OK, 1 not OK )
324 for(Int_t k=0; k<M; k++) A[k] = C1[k] + C2[k];
325 Bool_t err = InvertS(A,N);
328 Chi2 = Chi21 + Chi22;
331 MultSSQ( C1, A, K, N);
333 for( Int_t k=0; k<N;k++) r[k] = T1[k] - T2[k];
334 for( Int_t k=0; k<N;k++ )
335 for( Int_t l=0;l<N;l++) T[k] = T1[k] - K[k*N+l]*r[l];
337 for( Int_t ind=0,i=0; i<N; i++ ){
338 for( Int_t j=0; j<i; j++ ) Chi2+= 2*r[i]*r[j]*A[ind++];
339 Chi2+= r[i]*r[i]*A[ind++];
343 for( Int_t l=0; l<N; l++ ) K[ (N+1)*l ] -= 1;
345 for( Int_t ind = 0, l=0; l<N; ++l ){
346 for( Int_t j=0; j<=l; ++j, ind++ ){
348 for( Int_t k=0; k<N; ++k ) A[ind] -= K[l*N+k] * C1[IndexS(j,k)];
351 for( Int_t l=0; l<N; l++ ) C[l] = A[l];
355 //_____________________________________________________________________________
356 void AliHLTGlobalTrackMerger::MultSSQ( const Double_t *A, const Double_t *B, Double_t *C, Int_t N )
359 for( Int_t ind=0, i=0; i<N; ++i ){
360 for( Int_t j=0; j<N; ++j, ++ind ){
362 for( Int_t k=0; k<N; ++k ) C[ind] += A[IndexS(i,k)] * B[IndexS(k,j)];
367 //_____________________________________________________________________________
368 Bool_t AliHLTGlobalTrackMerger::InvertS( Double_t A[], Int_t N )
370 //* input: simmetric > 0 NxN matrix A = {a11,a21,a22,a31..a33,..}
371 //* output: inverse A, in case of problems fill zero and return 1
373 //* A->low triangular Anew : A = Anew x Anew^T
375 //* for(j=1,N) for(i=j,N) Aij=(Aii-sum_{k=1}^{j-1}Aik*Ajk )/Ajj
380 const Double_t ZERO = 1.E-20;
383 Double_t *j1 = A, *jj = A;
384 for( Int_t j=1; j<=N; j1+=j++, jj+=j ){
385 Double_t *ik = j1, x = 0;
396 for( Int_t step=1; step<=N-j; ik+=++step ){ // ik==Ai1
398 for( Double_t *jk=j1; jk!=jj; sum += (*(jk++)) * (*(ik++)) ) {}
399 *ik = (*ik - sum) * x; // ik == Aij
403 for( Int_t i=j; i<N; i++ ) *(ji+=i) = 0.;
413 //* for(j=1,i-1) Aij=-(sum_{k=j}^{i-1} Aik * Akj) / Aii ;
417 Double_t *ii=A,*ij=A;
418 for( Int_t i = 1; i<=N; ij=ii+1, ii+=++i ){
420 Double_t x = -(*ii = 1./ *ii);
423 for( Int_t j=1; j<i; jj+=++j, ij++ ){
424 Double_t *ik = ij, *kj = jj, sum = 0.;
425 for( Int_t k=j; ik!=ii; kj+=k++, ik++ ){
432 for( Double_t *ik = ij; ik!=ii+1; ik++ ){
442 //* Aij = sum_{k=i}^N Aki * Akj
445 Double_t *ii=A, *ij=A;
446 for( Int_t i=1; i<=N; ii+=++i ){
448 Double_t *ki = ii, *kj = ij, sum = 0.;
449 for( Int_t k=i; k<=N; ki+=k, kj+=k++ ) sum += (*ki) * (*kj);
451 }while( (ij++)!=ii );
457 //_____________________________________________________________________________
458 void AliHLTGlobalTrackMerger::PropagateTracksToDCA(const AliESDEvent *esdEvent)
460 // try to propagate all tracks to DCA to primary vertex
461 if(!esdEvent) return;
463 const Double_t kBz = esdEvent->GetMagneticField();
464 const Double_t kSmallRadius = 2.8; // [cm] something less than the beam pipe radius
465 const Double_t kMaxStep = 10.0; // [cm] track propagation step
466 Bool_t isOK = kFALSE;
468 Int_t nTracks = esdEvent->GetNumberOfTracks();
469 for(Int_t iTrack = 0; iTrack<nTracks; ++iTrack) {
470 AliESDtrack *track = esdEvent->GetTrack(iTrack);
473 // propagate to small radius (material budget included)
474 isOK = AliTracker::PropagateTrackTo(track,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);
476 // relate tracks to DCA to primary vertex
479 track->RelateToVertex(fVertex, kBz, kVeryBig);
480 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());
484 // the same procedure must be repeated for TPCinner (TPC only) tracks
486 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
487 if(!tpcTrack) continue;
489 // propagate tpcTracks to small radius (material budget included)
490 isOK = AliTracker::PropagateTrackTo(tpcTrack,kSmallRadius,track->GetMass(),kMaxStep,kFALSE);
492 // relate tracks to DCA to primary vertex
495 Double_t par[2], cov[3];
496 tpcTrack->PropagateToDCA(fVertex, kBz, kVeryBig,par,cov);
497 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());