Line data Source code
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 :
18 : //-----------------------------------------------------------------------------
19 : /// \class AliMUONTrackReconstructorK
20 : ///
21 : /// MUON track reconstructor using the kalman method
22 : ///
23 : /// This class contains as data:
24 : /// - the parameters for the track reconstruction
25 : ///
26 : /// It contains as methods, among others:
27 : /// - MakeTracks to build the tracks
28 : ///
29 : //-----------------------------------------------------------------------------
30 :
31 : #include "AliMUONTrackReconstructorK.h"
32 :
33 : #include "AliMUONConstants.h"
34 : #include "AliMUONVCluster.h"
35 : #include "AliMUONVClusterServer.h"
36 : #include "AliMUONVClusterStore.h"
37 : #include "AliMUONTrack.h"
38 : #include "AliMUONTrackParam.h"
39 : #include "AliMUONTrackExtrap.h"
40 : #include "AliMUONRecoParam.h"
41 : #include "AliMUONGeometryTransformer.h"
42 :
43 : #include "AliMpArea.h"
44 :
45 : #include "AliLog.h"
46 :
47 : #include <Riostream.h>
48 : #include <TMath.h>
49 : #include <TMatrixD.h>
50 : #include <TClonesArray.h>
51 :
52 : using std::endl;
53 : using std::cout;
54 : /// \cond CLASSIMP
55 18 : ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
56 : /// \endcond
57 :
58 : //__________________________________________________________________________
59 : AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(const AliMUONRecoParam* recoParam, AliMUONVClusterServer* clusterServer,
60 : const AliMUONGeometryTransformer* transformer)
61 2 : : AliMUONVTrackReconstructor(recoParam, clusterServer, transformer)
62 10 : {
63 : /// Constructor
64 4 : }
65 :
66 : //__________________________________________________________________________
67 : AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
68 8 : {
69 : /// Destructor
70 8 : }
71 :
72 : //__________________________________________________________________________
73 : Bool_t AliMUONTrackReconstructorK::MakeTrackCandidates(AliMUONVClusterStore& clusterStore)
74 : {
75 : /// To make track candidates (assuming linear propagation if AliMUONRecoParam::MakeTrackCandidatesFast() return kTRUE):
76 : /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
77 : /// Good candidates are made of at least three clusters if both stations are requested (two otherwise).
78 : /// Keep only best candidates or all of them according to the flag AliMUONRecoParam::TrackAllTracks().
79 :
80 : TClonesArray *segments;
81 : AliMUONObjectPair *segment;
82 : AliMUONTrack *track;
83 : Int_t iCandidate = 0;
84 : Bool_t clusterFound;
85 :
86 32 : AliDebug(1,"Enter MakeTrackCandidates");
87 :
88 : // Unless we're doing combined tracking, we'll clusterize all stations at once
89 : Int_t firstChamber(0);
90 : Int_t lastChamber(9);
91 :
92 8 : if (GetRecoParam()->CombineClusterTrackReco()) {
93 : // ... Here's the exception : ask the clustering to reconstruct
94 : // clusters *only* in station 4 and 5 for combined tracking
95 : firstChamber = 6;
96 0 : }
97 :
98 176 : for (Int_t i = firstChamber; i <= lastChamber; ++i )
99 : {
100 320 : if (fClusterServer && GetRecoParam()->UseChamber(i)) fClusterServer->Clusterize(i, clusterStore, AliMpArea(), GetRecoParam());
101 : }
102 :
103 : // Loop over stations(1..) 5 and 4 and make track candidates
104 56 : for (Int_t istat=4; istat>=3; istat--) {
105 :
106 : // Make segments in the station
107 16 : segments = MakeSegmentsBetweenChambers(clusterStore, 2*istat, 2*istat+1);
108 :
109 : // Loop over segments
110 116 : for (Int_t iSegment=0; iSegment<segments->GetEntriesFast(); iSegment++) {
111 102 : AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
112 34 : segment = (AliMUONObjectPair*) segments->UncheckedAt(iSegment);
113 :
114 : // Transform segments to tracks and put them at the end of fRecTracksPtr
115 34 : track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(segment,GetRecoParam()->GetBendingVertexDispersion());
116 34 : fNRecTracks++;
117 :
118 : // Look for compatible cluster(s) in the other station
119 34 : if (GetRecoParam()->MakeTrackCandidatesFast()) clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
120 34 : else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
121 :
122 : // Remove track if no cluster found on a requested station
123 : // or abort tracking if there are too many candidates
124 34 : if (GetRecoParam()->RequestStation(7-istat)) {
125 34 : if (!clusterFound) {
126 0 : fRecTracksPtr->Remove(track);
127 0 : fNRecTracks--;
128 34 : } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
129 0 : AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks));
130 0 : return kFALSE;
131 : }
132 : } else {
133 0 : if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
134 0 : AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
135 0 : return kFALSE;
136 : }
137 : }
138 :
139 : }
140 :
141 : }
142 :
143 24 : AliDebug(1,Form("Number of candidates before cleaning = %d",fNRecTracks));
144 :
145 : // Keep all different tracks if required
146 16 : if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
147 :
148 : // Retrace tracks using Kalman filter and select them if needed
149 8 : Int_t nCurrentTracks = fRecTracksPtr->GetLast()+1;
150 52 : for (Int_t iRecTrack = 0; iRecTrack < nCurrentTracks; iRecTrack++) {
151 18 : track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
152 :
153 : // skip empty slots
154 18 : if(!track) continue;
155 :
156 : // retrace tracks using Kalman filter and remove the ones for which extrap failed or that are out of limits
157 36 : if (!RetraceTrack(*track,kTRUE) || !IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
158 0 : fRecTracksPtr->Remove(track);
159 0 : fNRecTracks--;
160 0 : }
161 :
162 : }
163 :
164 : // Keep only the best tracks if required
165 8 : if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
166 8 : else fRecTracksPtr->Compress();
167 :
168 24 : AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
169 :
170 : return kTRUE;
171 :
172 8 : }
173 :
174 : //__________________________________________________________________________
175 : Bool_t AliMUONTrackReconstructorK::MakeMoreTrackCandidates(AliMUONVClusterStore& clusterStore)
176 : {
177 : /// To make extra track candidates assuming linear propagation:
178 : /// clustering is supposed to be already done
179 : /// Start with segments made of 1 cluster in each of the stations 4 and 5 then follow track in remaining chambers.
180 : /// Good candidates are made of at least three clusters if both stations are requested (two otherwise).
181 : /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
182 :
183 : TClonesArray *segments;
184 : AliMUONObjectPair *segment;
185 : AliMUONTrack *track;
186 : Int_t iCandidate = 0, iCurrentTrack, nCurrentTracks;
187 0 : Int_t initialNRecTracks = fNRecTracks;
188 : Bool_t clusterFound;
189 :
190 0 : AliDebug(1,"Enter MakeMoreTrackCandidates");
191 :
192 : // Double loop over chambers in stations(1..) 4 and 5 to make track candidates
193 0 : for (Int_t ich1 = 6; ich1 <= 7; ich1++) {
194 0 : for (Int_t ich2 = 8; ich2 <= 9; ich2++) {
195 :
196 : // Make segments between ch1 and ch2
197 0 : segments = MakeSegmentsBetweenChambers(clusterStore, ich1, ich2);
198 :
199 : /// Remove segments already attached to a track
200 0 : RemoveUsedSegments(*segments);
201 :
202 : // Loop over segments
203 0 : for (Int_t iSegment=0; iSegment<segments->GetEntriesFast(); iSegment++) {
204 0 : AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
205 0 : segment = (AliMUONObjectPair*) segments->UncheckedAt(iSegment);
206 :
207 : // Transform segments to tracks and put them at the end of fRecTracksPtr
208 0 : iCurrentTrack = fRecTracksPtr->GetLast()+1;
209 0 : track = new ((*fRecTracksPtr)[iCurrentTrack]) AliMUONTrack(segment,GetRecoParam()->GetBendingVertexDispersion());
210 0 : fNRecTracks++;
211 :
212 : // Look for compatible cluster(s) in the second chamber of station 5
213 0 : clusterFound = FollowLinearTrackInChamber(*track, clusterStore, 17-ich2);
214 :
215 : // skip the original track in case it has been removed
216 0 : if (GetRecoParam()->TrackAllTracks() && clusterFound) iCurrentTrack++;
217 :
218 : // loop over every new tracks
219 0 : nCurrentTracks = fRecTracksPtr->GetLast()+1;
220 0 : while (iCurrentTrack < nCurrentTracks) {
221 0 : track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iCurrentTrack);
222 :
223 : // Look for compatible cluster(s) in the second chamber of station 4
224 0 : FollowLinearTrackInChamber(*track, clusterStore, 13-ich1);
225 :
226 0 : iCurrentTrack++;
227 : }
228 :
229 : // abort tracking if there are too many candidates
230 0 : if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
231 0 : AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
232 0 : return kFALSE;
233 : }
234 :
235 : }
236 :
237 : }
238 : }
239 :
240 0 : AliDebug(1,Form("Number of candidates before cleaning = %d",fNRecTracks));
241 :
242 : // Retrace tracks using Kalman filter (also compute track chi2) and select them
243 0 : nCurrentTracks = fRecTracksPtr->GetLast()+1;
244 0 : for (Int_t iRecTrack = initialNRecTracks; iRecTrack < nCurrentTracks; iRecTrack++) {
245 0 : track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
246 :
247 : // skip empty slots
248 0 : if(!track) continue;
249 :
250 : // retrace tracks using Kalman filter and remove the ones for which extrap failed or that are out of limits
251 0 : if (!RetraceTrack(*track,kTRUE) || !IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
252 0 : fRecTracksPtr->Remove(track);
253 0 : fNRecTracks--;
254 0 : }
255 :
256 : }
257 :
258 : // Keep only the best tracks if required
259 0 : if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
260 0 : else fRecTracksPtr->Compress();
261 :
262 0 : AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
263 :
264 0 : return kTRUE;
265 :
266 0 : }
267 :
268 : //__________________________________________________________________________
269 : Bool_t AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed)
270 : {
271 : /// Re-run the kalman filter from the most downstream cluster to the most uptream one
272 : /// Return kFALSE in case of failure (i.e. extrapolation problem)
273 66 : AliDebug(1,"Enter RetraceTrack");
274 :
275 22 : AliMUONTrackParam* lastTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Last();
276 :
277 : // Reset the "seed" (= track parameters and their covariances at last cluster) if required
278 22 : if (resetSeed) {
279 :
280 : // parameters at last cluster
281 22 : AliMUONVCluster* cluster2 = lastTrackParam->GetClusterPtr();
282 22 : Double_t x2 = cluster2->GetX();
283 22 : Double_t y2 = cluster2->GetY();
284 22 : Double_t z2 = cluster2->GetZ();
285 :
286 : // parameters at last but one cluster
287 22 : AliMUONTrackParam* previousTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(lastTrackParam);
288 22 : AliMUONVCluster* cluster1 = previousTrackParam->GetClusterPtr();
289 : // make sure it is on the previous chamber (can have 2 clusters in the same chamber after "ComplementTrack")
290 22 : if (cluster2->GetChamberId() == cluster1->GetChamberId()) {
291 0 : previousTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(previousTrackParam);
292 0 : cluster1 = previousTrackParam->GetClusterPtr();
293 0 : }
294 22 : Double_t x1 = cluster1->GetX();
295 22 : Double_t y1 = cluster1->GetY();
296 22 : Double_t z1 = cluster1->GetZ();
297 :
298 : // reset track parameters
299 22 : Double_t dZ = z1 - z2;
300 22 : lastTrackParam->SetNonBendingCoor(x2);
301 22 : lastTrackParam->SetBendingCoor(y2);
302 22 : lastTrackParam->SetZ(z2);
303 22 : lastTrackParam->SetNonBendingSlope((x1 - x2) / dZ);
304 22 : lastTrackParam->SetBendingSlope((y1 - y2) / dZ);
305 22 : Double_t bendingImpact = y2 - z2 * lastTrackParam->GetBendingSlope();
306 22 : Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
307 22 : lastTrackParam->SetInverseBendingMomentum(inverseBendingMomentum);
308 :
309 : // => Reset track parameter covariances at last cluster (as if the other clusters did not exist)
310 22 : TMatrixD lastParamCov(5,5);
311 22 : lastParamCov.Zero();
312 : // Non bending plane
313 66 : lastParamCov(0,0) = cluster2->GetErrX2();
314 66 : lastParamCov(0,1) = - cluster2->GetErrX2() / dZ;
315 66 : lastParamCov(1,0) = lastParamCov(0,1);
316 88 : lastParamCov(1,1) = ( 1000. * cluster1->GetErrX2() + cluster2->GetErrX2() ) / dZ / dZ;
317 : // Bending plane
318 66 : lastParamCov(2,2) = cluster2->GetErrY2();
319 66 : lastParamCov(2,3) = - cluster2->GetErrY2() / dZ;
320 66 : lastParamCov(3,2) = lastParamCov(2,3);
321 88 : lastParamCov(3,3) = ( 1000. * cluster1->GetErrY2() + cluster2->GetErrY2() ) / dZ / dZ;
322 : // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
323 22 : if (AliMUONTrackExtrap::IsFieldON()) {
324 88 : lastParamCov(4,4) = ((GetRecoParam()->GetBendingVertexDispersion() *
325 44 : GetRecoParam()->GetBendingVertexDispersion() +
326 88 : (z1 * z1 * cluster2->GetErrY2() + z2 * z2 * 1000. * cluster1->GetErrY2()) / dZ / dZ) /
327 22 : bendingImpact / bendingImpact + 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum;
328 66 : lastParamCov(2,4) = z1 * cluster2->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
329 66 : lastParamCov(4,2) = lastParamCov(2,4);
330 110 : lastParamCov(3,4) = - (z1 * cluster2->GetErrY2() + z2 * 1000. * cluster1->GetErrY2()) *
331 22 : inverseBendingMomentum / bendingImpact / dZ / dZ;
332 66 : lastParamCov(4,3) = lastParamCov(3,4);
333 22 : } else lastParamCov(4,4) = inverseBendingMomentum*inverseBendingMomentum;
334 22 : lastTrackParam->SetCovariances(lastParamCov);
335 :
336 : // Reset the track chi2
337 22 : lastTrackParam->SetTrackChi2(0.);
338 :
339 22 : }
340 :
341 : // Redo the tracking
342 22 : return RetracePartialTrack(trackCandidate, lastTrackParam);
343 :
344 0 : }
345 :
346 : //__________________________________________________________________________
347 : Bool_t AliMUONTrackReconstructorK::RetracePartialTrack(AliMUONTrack &trackCandidate, const AliMUONTrackParam* startingTrackParam)
348 : {
349 : /// Re-run the kalman filter from the cluster attached to startingTrackParam to the most uptream cluster
350 : /// Return kFALSE in case of failure (i.e. extrapolation problem)
351 88 : AliDebug(1,"Enter RetracePartialTrack");
352 :
353 : // Printout for debuging
354 44 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
355 0 : cout << "RetracePartialTrack: track chi2 before re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
356 0 : }
357 :
358 : // Reset the track chi2
359 22 : trackCandidate.SetGlobalChi2(startingTrackParam->GetTrackChi2());
360 :
361 : // loop over attached clusters until the first one and recompute track parameters and covariances using kalman filter
362 : Bool_t extrapStatus = kTRUE;
363 22 : Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() - 1;
364 : Int_t currentChamber;
365 : Double_t addChi2TrackAtCluster;
366 22 : AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam);
367 232 : while (trackParamAtCluster) {
368 :
369 : // reset track parameters and their covariances
370 94 : trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
371 94 : trackParamAtCluster->SetZ(startingTrackParam->GetZ());
372 94 : trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
373 :
374 : // add MCS effect
375 94 : AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber+1),-1.);
376 :
377 : // reset propagator for smoother
378 188 : if (GetRecoParam()->UseSmoother()) trackParamAtCluster->ResetPropagator();
379 :
380 : // add MCS in missing chambers if any
381 94 : currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
382 188 : while (currentChamber < expectedChamber) {
383 : // extrapolation to the missing chamber (update the propagator)
384 0 : if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber),
385 0 : GetRecoParam()->UseSmoother())) extrapStatus = kFALSE;
386 : // add MCS effect
387 0 : AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber),-1.);
388 0 : expectedChamber--;
389 : }
390 :
391 : // extrapolation to the plane of the cluster attached to the current trackParamAtCluster (update the propagator)
392 188 : if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ(),
393 94 : GetRecoParam()->UseSmoother())) extrapStatus = kFALSE;
394 :
395 94 : if (GetRecoParam()->UseSmoother()) {
396 : // save extrapolated parameters for smoother
397 94 : trackParamAtCluster->SetExtrapParameters(trackParamAtCluster->GetParameters());
398 :
399 : // save extrapolated covariance matrix for smoother
400 94 : trackParamAtCluster->SetExtrapCovariances(trackParamAtCluster->GetCovariances());
401 94 : }
402 :
403 : // Compute new track parameters using kalman filter
404 94 : addChi2TrackAtCluster = RunKalmanFilter(*trackParamAtCluster);
405 :
406 : // Update the track chi2
407 94 : trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2() + addChi2TrackAtCluster);
408 94 : trackParamAtCluster->SetTrackChi2(trackCandidate.GetGlobalChi2());
409 :
410 : // prepare next step
411 94 : expectedChamber = currentChamber - 1;
412 : startingTrackParam = trackParamAtCluster;
413 94 : trackParamAtCluster = (AliMUONTrackParam*) (trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam));
414 : }
415 :
416 : // Printout for debuging
417 44 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
418 0 : cout << "RetracePartialTrack: track chi2 after re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
419 0 : }
420 :
421 : // set global chi2 to max value in case of problem during track extrapolation
422 22 : if (!extrapStatus) trackCandidate.SetGlobalChi2(2.*AliMUONTrack::MaxChi2());
423 22 : return extrapStatus;
424 :
425 : }
426 :
427 : //__________________________________________________________________________
428 : Bool_t AliMUONTrackReconstructorK::FollowTracks(AliMUONVClusterStore& clusterStore)
429 : {
430 : /// Follow tracks in stations(1..) 3, 2 and 1
431 32 : AliDebug(1,"Enter FollowTracks");
432 :
433 : AliMUONTrack *track;
434 : Int_t currentNRecTracks;
435 :
436 72 : for (Int_t station = 2; station >= 0; station--) {
437 :
438 : // Save the actual number of reconstructed track in case of
439 : // tracks are added or suppressed during the tracking procedure
440 : // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
441 24 : currentNRecTracks = fNRecTracks;
442 :
443 188 : for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
444 174 : AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
445 :
446 58 : track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
447 :
448 : // Look for compatible cluster(s) in station(0..) "station"
449 58 : if (!FollowTrackInStation(*track, clusterStore, station)) {
450 :
451 : // Try to recover track if required
452 0 : if (GetRecoParam()->RecoverTracks()) {
453 :
454 : // work on a copy of the track if this station is not required
455 : // to keep the case where no cluster is reconstructed as a possible candidate
456 0 : if (!GetRecoParam()->RequestStation(station)) {
457 0 : track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*track);
458 0 : fNRecTracks++;
459 0 : }
460 :
461 : // try to recover
462 0 : if (!RecoverTrack(*track, clusterStore, station)) {
463 : // remove track if no cluster found
464 0 : fRecTracksPtr->Remove(track);
465 0 : fNRecTracks--;
466 0 : }
467 :
468 0 : } else if (GetRecoParam()->RequestStation(station)) {
469 : // remove track if no cluster found
470 0 : fRecTracksPtr->Remove(track);
471 0 : fNRecTracks--;
472 0 : }
473 :
474 : }
475 :
476 : // abort tracking if there are too many candidates
477 116 : if (GetRecoParam()->RequestStation(station)) {
478 116 : if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
479 0 : AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks));
480 0 : return kFALSE;
481 : }
482 : } else {
483 0 : if ((fNRecTracks + currentNRecTracks - iRecTrack - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
484 0 : AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + currentNRecTracks - iRecTrack - 1));
485 0 : return kFALSE;
486 : }
487 : }
488 :
489 : }
490 :
491 24 : fRecTracksPtr->Compress(); // this is essential before checking tracks
492 :
493 72 : AliDebug(1,Form("In stations(1..) %d: Number of candidates before cleaning = %d",station+1,fNRecTracks));
494 :
495 : // Keep only the best tracks if required
496 24 : if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
497 :
498 72 : AliDebug(1,Form("In stations(1..) %d: Number of good candidates = %d",station+1,fNRecTracks));
499 :
500 : }
501 :
502 8 : return kTRUE;
503 :
504 8 : }
505 :
506 : //__________________________________________________________________________
507 : Bool_t AliMUONTrackReconstructorK::FollowTrackInChamber(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextChamber)
508 : {
509 : /// Follow trackCandidate in chamber(0..) nextChamber and search for compatible cluster(s)
510 : /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
511 : /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
512 : /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
513 : /// Remove the obsolete "trackCandidate" at the end.
514 : /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
515 : /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
516 0 : AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1));
517 :
518 : Double_t chi2OfCluster;
519 0 : Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
520 0 : GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
521 : Double_t addChi2TrackAtCluster;
522 0 : Double_t bestAddChi2TrackAtCluster = AliMUONTrack::MaxChi2();
523 : Bool_t foundOneCluster = kFALSE;
524 : AliMUONTrack *newTrack = 0x0;
525 : AliMUONVCluster *cluster;
526 0 : AliMUONTrackParam extrapTrackParamAtCh;
527 0 : AliMUONTrackParam extrapTrackParamAtCluster;
528 0 : AliMUONTrackParam bestTrackParamAtCluster;
529 :
530 : // Get track parameters according to the propagation direction
531 0 : if (nextChamber > 7) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
532 0 : else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
533 :
534 : // Printout for debuging
535 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
536 0 : cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
537 0 : extrapTrackParamAtCh.GetParameters().Print();
538 0 : extrapTrackParamAtCh.GetCovariances().Print();
539 : }
540 :
541 : // Add MCS effect
542 0 : Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
543 0 : AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
544 :
545 : // reset propagator for smoother
546 0 : if (GetRecoParam()->UseSmoother()) extrapTrackParamAtCh.ResetPropagator();
547 :
548 : // Add MCS in the missing chamber(s) if any
549 0 : while (currentChamber > nextChamber + 1) {
550 : // extrapolation to the missing chamber
551 0 : currentChamber--;
552 0 : if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber),
553 0 : GetRecoParam()->UseSmoother())) return kFALSE;
554 : // add MCS effect
555 0 : AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
556 : }
557 :
558 : //Extrapolate trackCandidate to chamber
559 0 : if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(nextChamber),
560 0 : GetRecoParam()->UseSmoother())) return kFALSE;
561 :
562 : // Printout for debuging
563 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
564 0 : cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(nextChamber)<<":"<<endl;
565 0 : extrapTrackParamAtCh.GetParameters().Print();
566 0 : extrapTrackParamAtCh.GetCovariances().Print();
567 : }
568 :
569 : // Printout for debuging
570 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
571 0 : cout << "FollowTrackInChamber: look for clusters in chamber(1..): " << nextChamber+1 << endl;
572 : }
573 :
574 : // Ask the clustering to reconstruct new clusters around the track position in the current chamber
575 : // except for station 4 and 5 that are already entirely clusterized
576 0 : if (GetRecoParam()->CombineClusterTrackReco()) {
577 0 : if (nextChamber < 6) AskForNewClustersInChamber(extrapTrackParamAtCh, clusterStore, nextChamber);
578 : }
579 :
580 : // Create iterators to loop over clusters in both chambers
581 0 : TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber));
582 :
583 : // look for candidates in chamber
584 0 : while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) ) {
585 :
586 : // try to add the current cluster fast
587 0 : if (!TryOneClusterFast(extrapTrackParamAtCh, cluster)) continue;
588 :
589 : // try to add the current cluster accuratly
590 0 : chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, cluster, extrapTrackParamAtCluster,
591 0 : GetRecoParam()->UseSmoother());
592 :
593 : // if good chi2 then consider to add cluster
594 0 : if (chi2OfCluster < maxChi2OfCluster) {
595 :
596 0 : if (GetRecoParam()->UseSmoother()) {
597 : // save extrapolated parameters for smoother
598 0 : extrapTrackParamAtCluster.SetExtrapParameters(extrapTrackParamAtCluster.GetParameters());
599 :
600 : // save extrapolated covariance matrix for smoother
601 0 : extrapTrackParamAtCluster.SetExtrapCovariances(extrapTrackParamAtCluster.GetCovariances());
602 : }
603 :
604 : // Compute new track parameters including new cluster using kalman filter
605 0 : addChi2TrackAtCluster = RunKalmanFilter(extrapTrackParamAtCluster);
606 :
607 : // skip track out of limits
608 0 : if (!IsAcceptable(extrapTrackParamAtCluster)) continue;
609 :
610 : // remember a cluster was found
611 : foundOneCluster = kTRUE;
612 :
613 : // Printout for debuging
614 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
615 0 : cout << "FollowTrackInChamber: found one cluster in chamber(1..): " << nextChamber+1
616 0 : << " (Chi2 = " << chi2OfCluster << ")" << endl;
617 0 : cluster->Print();
618 : }
619 :
620 0 : if (GetRecoParam()->TrackAllTracks()) {
621 : // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
622 0 : newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
623 0 : UpdateTrack(*newTrack,extrapTrackParamAtCluster,addChi2TrackAtCluster);
624 0 : fNRecTracks++;
625 :
626 : // Printout for debuging
627 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
628 0 : cout << "FollowTrackInChamber: added one cluster in chamber(1..): " << nextChamber+1 << endl;
629 0 : if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
630 : }
631 :
632 0 : } else if (addChi2TrackAtCluster < bestAddChi2TrackAtCluster) {
633 : // keep track of the best cluster
634 : bestAddChi2TrackAtCluster = addChi2TrackAtCluster;
635 0 : bestTrackParamAtCluster = extrapTrackParamAtCluster;
636 : }
637 :
638 : }
639 :
640 : }
641 :
642 : // fill out the best track if required else clean up the fRecTracksPtr array
643 0 : if (!GetRecoParam()->TrackAllTracks()) {
644 0 : if (foundOneCluster) {
645 0 : UpdateTrack(trackCandidate,bestTrackParamAtCluster,bestAddChi2TrackAtCluster);
646 :
647 : // Printout for debuging
648 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
649 0 : cout << "FollowTrackInChamber: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl;
650 0 : if (AliLog::GetGlobalDebugLevel() >= 3) trackCandidate.RecursiveDump();
651 : }
652 :
653 0 : } else return kFALSE;
654 :
655 0 : } else if (foundOneCluster) {
656 :
657 : // remove obsolete track
658 0 : fRecTracksPtr->Remove(&trackCandidate);
659 0 : fNRecTracks--;
660 :
661 0 : } else return kFALSE;
662 :
663 0 : return kTRUE;
664 :
665 0 : }
666 :
667 : //__________________________________________________________________________
668 : Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
669 : {
670 : /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
671 : /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
672 : /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
673 : /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
674 : /// Remove the obsolete "trackCandidate" at the end.
675 : /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
676 : /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
677 368 : AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
678 :
679 : // Order the chamber according to the propagation direction (tracking starts with chamber 2):
680 : // - nextStation == station(1...) 5 => forward propagation
681 : // - nextStation < station(1...) 5 => backward propagation
682 : Int_t ch1, ch2;
683 184 : if (nextStation==4) {
684 110 : ch1 = 2*nextStation+1;
685 : ch2 = 2*nextStation;
686 18 : } else {
687 : ch1 = 2*nextStation;
688 74 : ch2 = 2*nextStation+1;
689 : }
690 :
691 : Double_t chi2OfCluster;
692 184 : Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
693 92 : GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
694 : Double_t addChi2TrackAtCluster1;
695 : Double_t addChi2TrackAtCluster2;
696 92 : Double_t bestAddChi2TrackAtCluster1 = AliMUONTrack::MaxChi2();
697 92 : Double_t bestAddChi2TrackAtCluster2 = AliMUONTrack::MaxChi2();
698 : Bool_t foundOneCluster = kFALSE;
699 : Bool_t foundTwoClusters = kFALSE;
700 : AliMUONTrack *newTrack = 0x0;
701 : AliMUONVCluster *clusterCh1, *clusterCh2;
702 92 : AliMUONTrackParam extrapTrackParam;
703 92 : AliMUONTrackParam extrapTrackParamAtCh;
704 92 : AliMUONTrackParam extrapTrackParamAtCluster1;
705 92 : AliMUONTrackParam extrapTrackParamAtCluster2;
706 92 : AliMUONTrackParam bestTrackParamAtCluster1;
707 92 : AliMUONTrackParam bestTrackParamAtCluster2;
708 :
709 : // Get track parameters according to the propagation direction
710 146 : if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
711 222 : else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
712 :
713 : // Printout for debuging
714 368 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
715 0 : cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
716 0 : extrapTrackParamAtCh.GetParameters().Print();
717 0 : extrapTrackParamAtCh.GetCovariances().Print();
718 : }
719 :
720 : // Add MCS effect
721 92 : Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
722 92 : AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
723 :
724 : // reset propagator for smoother
725 184 : if (GetRecoParam()->UseSmoother()) extrapTrackParamAtCh.ResetPropagator();
726 :
727 : // Add MCS in the missing chamber(s) if any
728 170 : while (ch1 < ch2 && currentChamber > ch2 + 1) {
729 : // extrapolation to the missing chamber
730 2 : currentChamber--;
731 6 : if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber),
732 2 : GetRecoParam()->UseSmoother())) return kFALSE;
733 : // add MCS effect
734 2 : AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.);
735 : }
736 :
737 : //Extrapolate trackCandidate to chamber "ch2"
738 276 : if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2),
739 92 : GetRecoParam()->UseSmoother())) return kFALSE;
740 :
741 : // Printout for debuging
742 368 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
743 0 : cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
744 0 : extrapTrackParamAtCh.GetParameters().Print();
745 0 : extrapTrackParamAtCh.GetCovariances().Print();
746 : }
747 :
748 : // Printout for debuging
749 368 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
750 0 : cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
751 : }
752 :
753 : // Ask the clustering to reconstruct new clusters around the track position in the current station
754 : // except for station 4 and 5 that are already entirely clusterized
755 184 : if (GetRecoParam()->CombineClusterTrackReco()) {
756 92 : if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation);
757 : }
758 :
759 92 : Int_t nClusters = clusterStore.GetSize();
760 92 : Bool_t *clusterCh1Used = new Bool_t[nClusters];
761 3968 : for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
762 : Int_t iCluster1;
763 :
764 : // Create iterators to loop over clusters in both chambers
765 184 : TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
766 184 : TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
767 :
768 : // look for candidates in chamber 2
769 560 : while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
770 :
771 : // try to add the current cluster fast
772 376 : if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
773 :
774 : // try to add the current cluster accuratly
775 94 : chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2,
776 94 : GetRecoParam()->UseSmoother());
777 :
778 : // if good chi2 then try to attach a cluster in the other chamber too
779 94 : if (chi2OfCluster < maxChi2OfCluster) {
780 :
781 94 : if (GetRecoParam()->UseSmoother()) {
782 : // save extrapolated parameters for smoother
783 94 : extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
784 :
785 : // save extrapolated covariance matrix for smoother
786 188 : extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
787 : }
788 :
789 : // Compute new track parameters including "clusterCh2" using kalman filter
790 94 : addChi2TrackAtCluster2 = RunKalmanFilter(extrapTrackParamAtCluster2);
791 :
792 : // skip track out of limits
793 188 : if (!IsAcceptable(extrapTrackParamAtCluster2)) continue;
794 :
795 : // Printout for debuging
796 376 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
797 0 : cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
798 0 : << " (Chi2 = " << chi2OfCluster << ")" << endl;
799 0 : clusterCh2->Print();
800 0 : cout << " look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
801 : }
802 :
803 : // copy new track parameters for next step
804 94 : extrapTrackParam = extrapTrackParamAtCluster2;
805 :
806 : // add MCS effect
807 94 : AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
808 :
809 : // reset propagator for smoother
810 188 : if (GetRecoParam()->UseSmoother()) extrapTrackParam.ResetPropagator();
811 :
812 : //Extrapolate track parameters to chamber "ch1"
813 282 : Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1),
814 94 : GetRecoParam()->UseSmoother());
815 :
816 : // reset cluster iterator of chamber 1
817 94 : nextInCh1.Reset();
818 : iCluster1 = -1;
819 :
820 : // look for second candidates in chamber 1
821 : Bool_t foundSecondCluster = kFALSE;
822 666 : if (normalExtrap) while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
823 192 : iCluster1++;
824 :
825 : // try to add the current cluster fast
826 384 : if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
827 :
828 : // try to add the current cluster accuratly
829 96 : chi2OfCluster = TryOneCluster(extrapTrackParam, clusterCh1, extrapTrackParamAtCluster1,
830 96 : GetRecoParam()->UseSmoother());
831 :
832 : // if good chi2 then consider to add the 2 clusters to the "trackCandidate"
833 96 : if (chi2OfCluster < maxChi2OfCluster) {
834 :
835 94 : if (GetRecoParam()->UseSmoother()) {
836 : // save extrapolated parameters for smoother
837 94 : extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
838 :
839 : // save extrapolated covariance matrix for smoother
840 188 : extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
841 : }
842 :
843 : // Compute new track parameters including "clusterCh1" using kalman filter
844 94 : addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
845 :
846 : // skip track out of limits
847 188 : if (!IsAcceptable(extrapTrackParamAtCluster1)) continue;
848 :
849 : // remember a second cluster was found
850 : foundSecondCluster = kTRUE;
851 : foundTwoClusters = kTRUE;
852 :
853 : // Printout for debuging
854 376 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
855 0 : cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
856 0 : << " (Chi2 = " << chi2OfCluster << ")" << endl;
857 0 : clusterCh1->Print();
858 : }
859 :
860 94 : if (GetRecoParam()->TrackAllTracks()) {
861 : // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
862 376 : newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
863 94 : UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2,addChi2TrackAtCluster1,addChi2TrackAtCluster2);
864 94 : fNRecTracks++;
865 :
866 : // Tag clusterCh1 as used
867 94 : clusterCh1Used[iCluster1] = kTRUE;
868 :
869 : // Printout for debuging
870 376 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
871 0 : cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
872 0 : if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
873 : }
874 :
875 0 : } else if (addChi2TrackAtCluster1+addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1+bestAddChi2TrackAtCluster2) {
876 : // keep track of the best couple of clusters
877 : bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
878 : bestAddChi2TrackAtCluster2 = addChi2TrackAtCluster2;
879 0 : bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
880 0 : bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
881 : }
882 :
883 : }
884 :
885 : }
886 :
887 : // if no clusterCh1 found then consider to add clusterCh2 only
888 94 : if (!foundSecondCluster) {
889 :
890 : // remember a cluster was found
891 : foundOneCluster = kTRUE;
892 :
893 2 : if (GetRecoParam()->TrackAllTracks()) {
894 : // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
895 8 : newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
896 2 : UpdateTrack(*newTrack,extrapTrackParamAtCluster2,addChi2TrackAtCluster2);
897 2 : fNRecTracks++;
898 :
899 : // Printout for debuging
900 8 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
901 0 : cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
902 0 : if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
903 : }
904 :
905 0 : } else if (!foundTwoClusters && addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1) {
906 : // keep track of the best single cluster except if a couple of clusters has already been found
907 : bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster2;
908 0 : bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
909 : }
910 :
911 : }
912 :
913 94 : }
914 :
915 : }
916 :
917 : // look for candidates in chamber 1 not already attached to a track
918 : // if we want to keep all possible tracks or if no good couple of clusters has been found
919 92 : if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
920 :
921 : // add MCS effect for next step
922 92 : AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(ch2),-1.);
923 :
924 : //Extrapolate trackCandidate to chamber "ch1"
925 276 : Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1),
926 92 : GetRecoParam()->UseSmoother());
927 :
928 : // Printout for debuging
929 368 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
930 0 : cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch1)<<":"<<endl;
931 0 : extrapTrackParamAtCh.GetParameters().Print();
932 0 : extrapTrackParamAtCh.GetCovariances().Print();
933 : }
934 :
935 : // Printout for debuging
936 368 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
937 0 : cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
938 : }
939 :
940 : // reset cluster iterator of chamber 1
941 92 : nextInCh1.Reset();
942 : iCluster1 = -1;
943 :
944 : // look for second candidates in chamber 1
945 652 : if (normalExtrap) while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
946 188 : iCluster1++;
947 :
948 188 : if (clusterCh1Used[iCluster1]) continue; // Skip clusters already used
949 :
950 : // try to add the current cluster fast
951 188 : if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
952 :
953 : // try to add the current cluster accuratly
954 0 : chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1,
955 0 : GetRecoParam()->UseSmoother());
956 :
957 : // if good chi2 then consider to add clusterCh1
958 : // We do not try to attach a cluster in the other chamber too since it has already been done above
959 0 : if (chi2OfCluster < maxChi2OfCluster) {
960 :
961 0 : if (GetRecoParam()->UseSmoother()) {
962 : // save extrapolated parameters for smoother
963 0 : extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
964 :
965 : // save extrapolated covariance matrix for smoother
966 0 : extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
967 : }
968 :
969 : // Compute new track parameters including "clusterCh1" using kalman filter
970 0 : addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
971 :
972 : // skip track out of limits
973 0 : if (!IsAcceptable(extrapTrackParamAtCluster1)) continue;
974 :
975 : // remember a cluster was found
976 : foundOneCluster = kTRUE;
977 :
978 : // Printout for debuging
979 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
980 0 : cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
981 0 : << " (Chi2 = " << chi2OfCluster << ")" << endl;
982 0 : clusterCh1->Print();
983 : }
984 :
985 0 : if (GetRecoParam()->TrackAllTracks()) {
986 : // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
987 0 : newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
988 0 : UpdateTrack(*newTrack,extrapTrackParamAtCluster1,addChi2TrackAtCluster1);
989 0 : fNRecTracks++;
990 :
991 : // Printout for debuging
992 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
993 0 : cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
994 0 : if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
995 : }
996 :
997 0 : } else if (addChi2TrackAtCluster1 < bestAddChi2TrackAtCluster1) {
998 : // keep track of the best single cluster except if a couple of clusters has already been found
999 : bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
1000 0 : bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
1001 : }
1002 :
1003 : }
1004 :
1005 : }
1006 :
1007 92 : }
1008 :
1009 : // fill out the best track if required else clean up the fRecTracksPtr array
1010 92 : if (!GetRecoParam()->TrackAllTracks()) {
1011 0 : if (foundTwoClusters) {
1012 0 : UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2,bestAddChi2TrackAtCluster1,bestAddChi2TrackAtCluster2);
1013 :
1014 : // Printout for debuging
1015 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1016 0 : cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
1017 0 : if (AliLog::GetGlobalDebugLevel() >= 3) trackCandidate.RecursiveDump();
1018 : }
1019 :
1020 0 : } else if (foundOneCluster) {
1021 0 : UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestAddChi2TrackAtCluster1);
1022 :
1023 : // Printout for debuging
1024 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1025 0 : cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
1026 0 : if (AliLog::GetGlobalDebugLevel() >= 3) trackCandidate.RecursiveDump();
1027 : }
1028 :
1029 : } else {
1030 0 : delete [] clusterCh1Used;
1031 0 : return kFALSE;
1032 : }
1033 :
1034 182 : } else if (foundOneCluster || foundTwoClusters) {
1035 :
1036 : // remove obsolete track
1037 92 : fRecTracksPtr->Remove(&trackCandidate);
1038 92 : fNRecTracks--;
1039 :
1040 : } else {
1041 0 : delete [] clusterCh1Used;
1042 0 : return kFALSE;
1043 : }
1044 :
1045 184 : delete [] clusterCh1Used;
1046 92 : return kTRUE;
1047 :
1048 184 : }
1049 :
1050 : //__________________________________________________________________________
1051 : Double_t AliMUONTrackReconstructorK::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster)
1052 : {
1053 : /// Compute new track parameters and their covariances including new cluster using kalman filter
1054 : /// return the additional track chi2
1055 1144 : AliDebug(1,"Enter RunKalmanFilter");
1056 :
1057 : // Get actual track parameters (p)
1058 286 : TMatrixD param(trackParamAtCluster.GetParameters());
1059 :
1060 : // Get new cluster parameters (m)
1061 286 : AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
1062 286 : TMatrixD clusterParam(5,1);
1063 286 : clusterParam.Zero();
1064 858 : clusterParam(0,0) = cluster->GetX();
1065 858 : clusterParam(2,0) = cluster->GetY();
1066 :
1067 : // Compute the actual parameter weight (W)
1068 572 : TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
1069 572 : if (paramWeight.Determinant() != 0) {
1070 286 : paramWeight.Invert();
1071 : } else {
1072 0 : AliWarning(" Determinant = 0");
1073 0 : return 2.*AliMUONTrack::MaxChi2();
1074 : }
1075 :
1076 : // Compute the new cluster weight (U)
1077 286 : TMatrixD clusterWeight(5,5);
1078 286 : clusterWeight.Zero();
1079 858 : clusterWeight(0,0) = 1. / cluster->GetErrX2();
1080 858 : clusterWeight(2,2) = 1. / cluster->GetErrY2();
1081 :
1082 : // Compute the new parameters covariance matrix ( (W+U)^-1 )
1083 286 : TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
1084 572 : if (newParamCov.Determinant() != 0) {
1085 286 : newParamCov.Invert();
1086 : } else {
1087 0 : AliWarning(" Determinant = 0");
1088 0 : return 2.*AliMUONTrack::MaxChi2();
1089 : }
1090 :
1091 : // Save the new parameters covariance matrix
1092 286 : trackParamAtCluster.SetCovariances(newParamCov);
1093 :
1094 : // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
1095 286 : TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
1096 286 : TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
1097 286 : TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
1098 286 : newParam += param; // ((W+U)^-1)U(m-p) + p
1099 :
1100 : // Save the new parameters
1101 286 : trackParamAtCluster.SetParameters(newParam);
1102 :
1103 : // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
1104 286 : tmp = newParam; // p'
1105 286 : tmp -= param; // (p'-p)
1106 286 : TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
1107 286 : TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
1108 286 : tmp = newParam; // p'
1109 286 : tmp -= clusterParam; // (p'-m)
1110 286 : TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
1111 858 : addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
1112 :
1113 572 : return addChi2Track(0,0);
1114 :
1115 858 : }
1116 :
1117 : //__________________________________________________________________________
1118 : void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster, Double_t addChi2)
1119 : {
1120 : /// Add 1 cluster to the track candidate
1121 : /// Update chi2 of the track
1122 :
1123 : // Flag cluster as being (not) removable
1124 4 : if (GetRecoParam()->RequestStation(trackParamAtCluster.GetClusterPtr()->GetChamberId()/2))
1125 2 : trackParamAtCluster.SetRemovable(kFALSE);
1126 0 : else trackParamAtCluster.SetRemovable(kTRUE);
1127 2 : trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
1128 :
1129 : // Update the track chi2 into trackParamAtCluster
1130 2 : trackParamAtCluster.SetTrackChi2(track.GetGlobalChi2() + addChi2);
1131 :
1132 : // Update the chi2 of the new track
1133 2 : track.SetGlobalChi2(trackParamAtCluster.GetTrackChi2());
1134 :
1135 : // Update array of TrackParamAtCluster
1136 2 : track.AddTrackParamAtCluster(trackParamAtCluster,*(trackParamAtCluster.GetClusterPtr()));
1137 :
1138 2 : }
1139 :
1140 : //__________________________________________________________________________
1141 : void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2,
1142 : Double_t addChi2AtCluster1, Double_t addChi2AtCluster2)
1143 : {
1144 : /// Add 2 clusters to the track candidate (order is important)
1145 : /// Update track and local chi2
1146 :
1147 : // Update local chi2 at first cluster
1148 188 : AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
1149 94 : Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
1150 94 : Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
1151 188 : Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
1152 94 : deltaY*deltaY / cluster1->GetErrY2();
1153 94 : trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
1154 :
1155 : // Flag first cluster as being removable
1156 94 : trackParamAtCluster1.SetRemovable(kTRUE);
1157 :
1158 : // Update local chi2 at second cluster
1159 94 : AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
1160 94 : AliMUONTrackParam extrapTrackParamAtCluster2(trackParamAtCluster1);
1161 94 : AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParamAtCluster2, trackParamAtCluster2.GetZ());
1162 282 : deltaX = extrapTrackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
1163 282 : deltaY = extrapTrackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
1164 282 : Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
1165 188 : deltaY*deltaY / cluster2->GetErrY2();
1166 94 : trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
1167 :
1168 : // Flag second cluster as being removable
1169 94 : trackParamAtCluster2.SetRemovable(kTRUE);
1170 :
1171 : // Update the track chi2 into trackParamAtCluster2
1172 94 : trackParamAtCluster2.SetTrackChi2(track.GetGlobalChi2() + addChi2AtCluster2);
1173 :
1174 : // Update the track chi2 into trackParamAtCluster1
1175 94 : trackParamAtCluster1.SetTrackChi2(trackParamAtCluster2.GetTrackChi2() + addChi2AtCluster1);
1176 :
1177 : // Update the chi2 of the new track
1178 94 : track.SetGlobalChi2(trackParamAtCluster1.GetTrackChi2());
1179 :
1180 : // Update array of trackParamAtCluster
1181 94 : track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
1182 94 : track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
1183 :
1184 94 : }
1185 :
1186 : //__________________________________________________________________________
1187 : Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
1188 : {
1189 : /// Try to recover the track candidate in the next station
1190 : /// by removing the worst of the two clusters attached in the current station
1191 : /// Return kTRUE if recovering succeeds
1192 0 : AliDebug(1,"Enter RecoverTrack");
1193 :
1194 : // Do not try to recover track until we have attached cluster(s) on station(1..) 3
1195 0 : if (nextStation > 1) return kFALSE;
1196 :
1197 : Int_t worstClusterNumber = -1;
1198 : Double_t localChi2, worstLocalChi2 = -1.;
1199 :
1200 : // Look for the cluster to remove
1201 0 : for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
1202 0 : AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
1203 :
1204 : // check if current cluster is in the previous station
1205 0 : if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break;
1206 :
1207 : // check if current cluster is removable
1208 0 : if (!trackParamAtCluster->IsRemovable()) return kFALSE;
1209 :
1210 : // reset the current cluster as being not removable if it is on a required station
1211 0 : if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE);
1212 :
1213 : // Pick up cluster with the worst chi2
1214 0 : localChi2 = trackParamAtCluster->GetLocalChi2();
1215 0 : if (localChi2 > worstLocalChi2) {
1216 : worstLocalChi2 = localChi2;
1217 : worstClusterNumber = clusterNumber;
1218 0 : }
1219 0 : }
1220 :
1221 : // check if worst cluster found
1222 0 : if (worstClusterNumber < 0) return kFALSE;
1223 :
1224 : // Remove the worst cluster
1225 0 : trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
1226 :
1227 : // Re-calculate track parameters at the (new) first cluster
1228 0 : if (!RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(1))) return kFALSE;
1229 :
1230 : // skip track out of limits
1231 0 : if (!IsAcceptable(*((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()))) return kFALSE;
1232 :
1233 : // Look for new cluster(s) in next station
1234 0 : return FollowTrackInStation(trackCandidate, clusterStore, nextStation);
1235 :
1236 0 : }
1237 :
1238 : //__________________________________________________________________________
1239 : Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track)
1240 : {
1241 : /// Compute new track parameters and their covariances using smoother
1242 72 : AliDebug(1,"Enter UseSmoother");
1243 :
1244 18 : AliMUONTrackParam *previousTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1245 :
1246 : // Smoothed parameters and covariances at first cluster = filtered parameters and covariances
1247 18 : previousTrackParam->SetSmoothParameters(previousTrackParam->GetParameters());
1248 18 : previousTrackParam->SetSmoothCovariances(previousTrackParam->GetCovariances());
1249 :
1250 18 : AliMUONTrackParam *currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
1251 :
1252 : // Save local chi2 at first cluster = last additional chi2 provided by Kalman
1253 18 : previousTrackParam->SetLocalChi2(previousTrackParam->GetTrackChi2() - currentTrackParam->GetTrackChi2());
1254 :
1255 : // if the track contains only 2 clusters simply copy the filtered parameters
1256 18 : if (track.GetNClusters() == 2) {
1257 0 : currentTrackParam->SetSmoothParameters(currentTrackParam->GetParameters());
1258 0 : currentTrackParam->SetSmoothCovariances(currentTrackParam->GetCovariances());
1259 0 : currentTrackParam->SetLocalChi2(currentTrackParam->GetTrackChi2());
1260 0 : return kTRUE;
1261 : }
1262 :
1263 180 : while (currentTrackParam) {
1264 :
1265 : // Get variables
1266 162 : const TMatrixD &extrapParameters = previousTrackParam->GetExtrapParameters(); // X(k+1 k)
1267 162 : const TMatrixD &filteredParameters = currentTrackParam->GetParameters(); // X(k k)
1268 162 : const TMatrixD &previousSmoothParameters = previousTrackParam->GetSmoothParameters(); // X(k+1 n)
1269 162 : const TMatrixD &propagator = previousTrackParam->GetPropagator(); // F(k)
1270 162 : const TMatrixD &extrapCovariances = previousTrackParam->GetExtrapCovariances(); // C(k+1 k)
1271 162 : const TMatrixD &filteredCovariances = currentTrackParam->GetCovariances(); // C(k k)
1272 162 : const TMatrixD &previousSmoothCovariances = previousTrackParam->GetSmoothCovariances(); // C(k+1 n)
1273 :
1274 : // Compute smoother gain: A(k) = C(kk) * F(k)^t * (C(k+1 k))^-1
1275 162 : TMatrixD extrapWeight(extrapCovariances);
1276 324 : if (extrapWeight.Determinant() != 0) {
1277 162 : extrapWeight.Invert(); // (C(k+1 k))^-1
1278 : } else {
1279 0 : AliWarning(" Determinant = 0");
1280 0 : return kFALSE;
1281 : }
1282 162 : TMatrixD smootherGain(filteredCovariances,TMatrixD::kMultTranspose,propagator); // C(kk) * F(k)^t
1283 162 : smootherGain *= extrapWeight; // C(kk) * F(k)^t * (C(k+1 k))^-1
1284 :
1285 : // Compute smoothed parameters: X(k n) = X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
1286 162 : TMatrixD tmpParam(previousSmoothParameters,TMatrixD::kMinus,extrapParameters); // X(k+1 n) - X(k+1 k)
1287 162 : TMatrixD smoothParameters(smootherGain,TMatrixD::kMult,tmpParam); // A(k) * (X(k+1 n) - X(k+1 k))
1288 162 : smoothParameters += filteredParameters; // X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
1289 :
1290 : // Save smoothed parameters
1291 162 : currentTrackParam->SetSmoothParameters(smoothParameters);
1292 :
1293 : // Compute smoothed covariances: C(k n) = C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
1294 162 : TMatrixD tmpCov(previousSmoothCovariances,TMatrixD::kMinus,extrapCovariances); // C(k+1 n) - C(k+1 k)
1295 162 : TMatrixD tmpCov2(tmpCov,TMatrixD::kMultTranspose,smootherGain); // (C(k+1 n) - C(k+1 k)) * (A(k))^t
1296 162 : TMatrixD smoothCovariances(smootherGain,TMatrixD::kMult,tmpCov2); // A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
1297 162 : smoothCovariances += filteredCovariances; // C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
1298 :
1299 : // Save smoothed covariances
1300 162 : currentTrackParam->SetSmoothCovariances(smoothCovariances);
1301 :
1302 : // Compute smoothed residual: r(k n) = cluster - X(k n)
1303 162 : AliMUONVCluster* cluster = currentTrackParam->GetClusterPtr();
1304 162 : TMatrixD smoothResidual(2,1);
1305 162 : smoothResidual.Zero();
1306 648 : smoothResidual(0,0) = cluster->GetX() - smoothParameters(0,0);
1307 648 : smoothResidual(1,0) = cluster->GetY() - smoothParameters(2,0);
1308 :
1309 : // Compute weight of smoothed residual: W(k n) = (clusterCov - C(k n))^-1
1310 162 : TMatrixD smoothResidualWeight(2,2);
1311 648 : smoothResidualWeight(0,0) = cluster->GetErrX2() - smoothCovariances(0,0);
1312 486 : smoothResidualWeight(0,1) = - smoothCovariances(0,2);
1313 486 : smoothResidualWeight(1,0) = - smoothCovariances(2,0);
1314 648 : smoothResidualWeight(1,1) = cluster->GetErrY2() - smoothCovariances(2,2);
1315 324 : if (smoothResidualWeight.Determinant() != 0) {
1316 162 : smoothResidualWeight.Invert();
1317 : } else {
1318 0 : AliWarning(" Determinant = 0");
1319 0 : return kFALSE;
1320 : }
1321 :
1322 : // Compute local chi2 = (r(k n))^t * W(k n) * r(k n)
1323 162 : TMatrixD tmpChi2(smoothResidual,TMatrixD::kTransposeMult,smoothResidualWeight); // (r(k n))^t * W(k n)
1324 162 : TMatrixD localChi2(tmpChi2,TMatrixD::kMult,smoothResidual); // (r(k n))^t * W(k n) * r(k n)
1325 :
1326 : // Save local chi2
1327 324 : currentTrackParam->SetLocalChi2(localChi2(0,0));
1328 :
1329 : previousTrackParam = currentTrackParam;
1330 486 : currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
1331 486 : }
1332 :
1333 18 : return kTRUE;
1334 :
1335 18 : }
1336 :
1337 : //__________________________________________________________________________
1338 : Bool_t AliMUONTrackReconstructorK::ComplementTracks(const AliMUONVClusterStore& clusterStore)
1339 : {
1340 : /// Complete tracks by adding missing clusters (if there is an overlap between
1341 : /// two detection elements, the track may have two clusters in the same chamber).
1342 : /// Recompute track parameters and covariances at each clusters.
1343 : /// Remove tracks getting abnormal (i.e. extrapolation failed...) after being complemented.
1344 : /// Return kTRUE if one or more tracks have been complemented or removed.
1345 32 : AliDebug(1,"Enter ComplementTracks");
1346 :
1347 : Int_t chamberId, detElemId;
1348 : Double_t chi2OfCluster, addChi2TrackAtCluster, bestAddChi2TrackAtCluster;
1349 16 : Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
1350 8 : GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
1351 : Bool_t foundOneCluster, trackModified, hasChanged = kFALSE;
1352 : AliMUONVCluster *cluster;
1353 8 : AliMUONTrackParam *trackParam, *previousTrackParam, *nextTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
1354 : AliMUONTrack *nextTrack;
1355 :
1356 16 : AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
1357 56 : while (track) {
1358 : trackModified = kFALSE;
1359 :
1360 60 : trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1361 : previousTrackParam = trackParam;
1362 436 : while (trackParam) {
1363 : foundOneCluster = kFALSE;
1364 198 : bestAddChi2TrackAtCluster = AliMUONTrack::MaxChi2();
1365 198 : chamberId = trackParam->GetClusterPtr()->GetChamberId();
1366 198 : detElemId = trackParam->GetClusterPtr()->GetDetElemId();
1367 :
1368 : // prepare nextTrackParam before adding new cluster because of the sorting
1369 594 : nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
1370 :
1371 : // Create iterators to loop over clusters in current chamber
1372 396 : TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
1373 :
1374 : // look for one second candidate in the same chamber
1375 1212 : while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
1376 :
1377 : // look for a cluster in another detection element
1378 816 : if (cluster->GetDetElemId() == detElemId) continue;
1379 :
1380 : // try to add the current cluster fast
1381 340 : if (!TryOneClusterFast(*trackParam, cluster)) continue;
1382 :
1383 : // try to add the current cluster accurately
1384 : // never use track parameters at last cluster because the covariance matrix is meaningless
1385 24 : if (nextTrackParam) chi2OfCluster = TryOneCluster(*trackParam, cluster, trackParamAtCluster);
1386 0 : else chi2OfCluster = TryOneCluster(*previousTrackParam, cluster, trackParamAtCluster);
1387 :
1388 : // if good chi2 then consider to add this cluster to the track
1389 8 : if (chi2OfCluster < maxChi2OfCluster) {
1390 :
1391 : // Compute local track parameters including current cluster using kalman filter
1392 4 : addChi2TrackAtCluster = RunKalmanFilter(trackParamAtCluster);
1393 :
1394 : // keep track of the best cluster
1395 4 : if (addChi2TrackAtCluster < bestAddChi2TrackAtCluster) {
1396 : bestAddChi2TrackAtCluster = addChi2TrackAtCluster;
1397 4 : bestTrackParamAtCluster = trackParamAtCluster;
1398 : foundOneCluster = kTRUE;
1399 4 : }
1400 :
1401 : }
1402 :
1403 : }
1404 :
1405 : // add new cluster if any
1406 198 : if (foundOneCluster) {
1407 :
1408 : // Printout for debuging
1409 16 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1410 0 : cout << "ComplementTracks: found one cluster in chamber(1..): " << chamberId+1 << endl;
1411 0 : bestTrackParamAtCluster.GetClusterPtr()->Print();
1412 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
1413 0 : cout<<endl<<"Track parameters and covariances at cluster:"<<endl;
1414 0 : bestTrackParamAtCluster.GetParameters().Print();
1415 0 : bestTrackParamAtCluster.GetCovariances().Print();
1416 : }
1417 : }
1418 :
1419 4 : trackParam->SetRemovable(kTRUE);
1420 4 : bestTrackParamAtCluster.SetRemovable(kTRUE);
1421 4 : track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
1422 : trackModified = kTRUE;
1423 : hasChanged = kTRUE;
1424 4 : }
1425 :
1426 : previousTrackParam = trackParam;
1427 : trackParam = nextTrackParam;
1428 198 : }
1429 :
1430 : // prepare next track
1431 40 : nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1432 :
1433 : // re-compute track parameters using kalman filter if needed
1434 28 : if (trackModified && !RetraceTrack(*track,kTRUE)) {
1435 0 : AliWarning("track modified but problem occur during refitting --> remove track");
1436 0 : fRecTracksPtr->Remove(track);
1437 0 : fNRecTracks--;
1438 0 : }
1439 :
1440 : track = nextTrack;
1441 : }
1442 :
1443 8 : return hasChanged;
1444 :
1445 8 : }
1446 :
1447 : //__________________________________________________________________________
1448 : void AliMUONTrackReconstructorK::ImproveTrack(AliMUONTrack &track)
1449 : {
1450 : /// Improve the given track by removing removable clusters with local chi2 highter than the defined cut
1451 : /// Removable clusters are identified by the method AliMUONTrack::TagRemovableClusters()
1452 : /// Recompute track parameters and covariances at the remaining clusters
1453 : /// and if something goes wrong (i.e. extrapolation failed...) set track chi2 to max value
1454 72 : AliDebug(1,"Enter ImproveTrack");
1455 :
1456 : Double_t localChi2, worstLocalChi2;
1457 : AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster, *nextTrackParam, *next2nextTrackParam;
1458 : Int_t nextChamber, next2nextChamber;
1459 : Bool_t smoothed;
1460 36 : Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForImprovement() *
1461 18 : GetRecoParam()->GetSigmaCutForImprovement();
1462 :
1463 36 : while (!track.IsImproved()) {
1464 :
1465 : // identify removable clusters
1466 18 : track.TagRemovableClusters(GetRecoParam()->RequestedStationMask());
1467 :
1468 : // Run smoother if required
1469 : smoothed = kFALSE;
1470 36 : if (GetRecoParam()->UseSmoother()) smoothed = RunSmoother(track);
1471 :
1472 : // Use standard procedure to compute local chi2 if smoother not required or not working
1473 18 : if (!smoothed) {
1474 :
1475 : // Update track parameters and covariances
1476 0 : if (!track.UpdateCovTrackParamAtCluster()) {
1477 0 : AliWarning("unable to update track parameters and covariances --> stop improvement");
1478 : // restore the kalman parameters
1479 0 : RetraceTrack(track,kTRUE);
1480 0 : break;
1481 : }
1482 :
1483 : // Compute local chi2 of each clusters
1484 0 : track.ComputeLocalChi2(kTRUE);
1485 0 : }
1486 :
1487 : // Look for the cluster to remove
1488 : worstTrackParamAtCluster = 0x0;
1489 : worstLocalChi2 = -1.;
1490 18 : trackParamAtCluster = (AliMUONTrackParam*)track.GetTrackParamAtCluster()->First();
1491 396 : while (trackParamAtCluster) {
1492 :
1493 : // save parameters into smooth parameters in case of smoother did not work properly
1494 360 : if (GetRecoParam()->UseSmoother() && !smoothed) {
1495 0 : trackParamAtCluster->SetSmoothParameters(trackParamAtCluster->GetParameters());
1496 0 : trackParamAtCluster->SetSmoothCovariances(trackParamAtCluster->GetCovariances());
1497 0 : }
1498 :
1499 : // Pick up cluster with the worst chi2
1500 180 : localChi2 = trackParamAtCluster->GetLocalChi2();
1501 180 : if (localChi2 > worstLocalChi2) {
1502 : worstLocalChi2 = localChi2;
1503 : worstTrackParamAtCluster = trackParamAtCluster;
1504 60 : }
1505 :
1506 180 : trackParamAtCluster = (AliMUONTrackParam*)track.GetTrackParamAtCluster()->After(trackParamAtCluster);
1507 : }
1508 :
1509 : // Check whether the worst chi2 is under requirement or not
1510 18 : if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1511 18 : track.SetImproved(kTRUE);
1512 18 : break;
1513 : }
1514 :
1515 : // if the worst cluster is not removable then stop improvement
1516 0 : if (!worstTrackParamAtCluster->IsRemovable()) {
1517 : // restore the kalman parameters in case they have been lost
1518 0 : if (!smoothed) RetraceTrack(track,kTRUE);
1519 : break;
1520 : }
1521 :
1522 : // get track parameters at cluster next to the one to be removed
1523 0 : nextTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(worstTrackParamAtCluster);
1524 :
1525 : // Remove the worst cluster
1526 0 : track.RemoveTrackParamAtCluster(worstTrackParamAtCluster);
1527 :
1528 : // Re-calculate track parameters
1529 : // - from the cluster immediately downstream the one suppressed
1530 : // - or from the begining - if parameters have been re-computed using the standard method (kalman parameters have been lost)
1531 : // - or if the removed cluster was used to compute the tracking seed
1532 : Bool_t normalExtrap;
1533 0 : if (smoothed && nextTrackParam) {
1534 :
1535 0 : nextChamber = nextTrackParam->GetClusterPtr()->GetChamberId();
1536 : next2nextTrackParam = nextTrackParam;
1537 0 : do {
1538 :
1539 0 : next2nextChamber = next2nextTrackParam->GetClusterPtr()->GetChamberId();
1540 0 : next2nextTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(next2nextTrackParam);
1541 :
1542 0 : } while (next2nextTrackParam && (next2nextChamber == nextChamber));
1543 :
1544 0 : if (next2nextChamber == nextChamber) normalExtrap = RetraceTrack(track,kTRUE);
1545 0 : else normalExtrap = RetracePartialTrack(track,nextTrackParam);
1546 :
1547 0 : } else normalExtrap = RetraceTrack(track,kTRUE);
1548 :
1549 : // stop in case of extrapolation problem
1550 0 : if (!normalExtrap) {
1551 0 : AliWarning("track partially improved but problem occur during refitting --> stop improvement");
1552 0 : break;
1553 : }
1554 :
1555 : // Printout for debuging
1556 0 : if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1557 0 : cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl;
1558 0 : }
1559 :
1560 0 : }
1561 :
1562 18 : }
1563 :
1564 : //__________________________________________________________________________
1565 : Bool_t AliMUONTrackReconstructorK::FinalizeTrack(AliMUONTrack &track)
1566 : {
1567 : /// Update track parameters and covariances at each attached cluster
1568 : /// using smoother if required, if not already done
1569 : /// return kFALSE if the track cannot be extrapolated uo to the last chamber
1570 :
1571 : AliMUONTrackParam *trackParamAtCluster;
1572 : Bool_t smoothed = kFALSE;
1573 :
1574 : // update track parameters (using smoother if required) if not already done
1575 48 : if (track.IsImproved()) smoothed = GetRecoParam()->UseSmoother();
1576 : else {
1577 0 : if (GetRecoParam()->UseSmoother()) smoothed = RunSmoother(track);
1578 0 : if (!smoothed) {
1579 0 : if (track.UpdateCovTrackParamAtCluster()) track.ComputeLocalChi2(kTRUE);
1580 : else {
1581 0 : AliWarning("finalization failed due to extrapolation problem");
1582 0 : return kFALSE;
1583 : }
1584 0 : }
1585 : }
1586 :
1587 : // copy smoothed parameters and covariances if any
1588 16 : if (smoothed) {
1589 :
1590 16 : trackParamAtCluster = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First());
1591 356 : while (trackParamAtCluster) {
1592 :
1593 162 : trackParamAtCluster->SetParameters(trackParamAtCluster->GetSmoothParameters());
1594 162 : trackParamAtCluster->SetCovariances(trackParamAtCluster->GetSmoothCovariances());
1595 :
1596 162 : trackParamAtCluster = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->After(trackParamAtCluster));
1597 : }
1598 :
1599 : }
1600 :
1601 16 : return kTRUE;
1602 :
1603 16 : }
1604 :
1605 : //__________________________________________________________________________
1606 : Bool_t AliMUONTrackReconstructorK::RefitTrack(AliMUONTrack &track, Bool_t enableImprovement)
1607 : {
1608 : /// re-fit the given track
1609 0 : AliDebug(1,"Enter RefitTrack");
1610 :
1611 : // check validity of the track (i.e. at least 2 chambers hit on stations 4 and 5)
1612 0 : if (!track.IsValid(0)) {
1613 0 : AliWarning("the track is not valid --> unable to refit");
1614 0 : return kFALSE;
1615 : }
1616 :
1617 : // re-compute track parameters and covariances using Kalman filter
1618 0 : if (!RetraceTrack(track,kTRUE)) {
1619 0 : AliWarning("bad track refitting due to extrapolation failure");
1620 0 : return kFALSE;
1621 : }
1622 :
1623 : // Improve the reconstructed tracks if required
1624 0 : track.SetImproved(kFALSE);
1625 0 : if (enableImprovement && GetRecoParam()->ImproveTracks()) ImproveTrack(track);
1626 :
1627 : // Fill AliMUONTrack data members
1628 0 : if (track.GetGlobalChi2() < AliMUONTrack::MaxChi2()) return FinalizeTrack(track);
1629 : else {
1630 0 : AliWarning("track not finalized due to extrapolation failure");
1631 0 : return kFALSE;
1632 : }
1633 :
1634 0 : }
1635 :
|
