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