]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/AliMUONTrackReconstructorK.cxx
- Removing use of volpath.dat file, now we can use volume symnames
[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
06ca6d7b 18/// \class AliMUONTrackReconstructorK
19/////////////////////////////////////
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///
8cc77df0 29////////////////////////////////////
30
8cc77df0 31#include "AliMUONTrackReconstructorK.h"
8cc77df0 32#include "AliMUONConstants.h"
33#include "AliMUONHitForRec.h"
ea94c18b 34#include "AliMUONTrack.h"
35#include "AliMUONTrackParam.h"
36#include "AliMUONTrackExtrap.h"
8cde4af5 37
8cc77df0 38#include "AliLog.h"
39
8cde4af5 40#include <Riostream.h>
ea94c18b 41#include <TMath.h>
42#include <TMatrixD.h>
8cde4af5 43
44/// \cond CLASSIMP
8cc77df0 45ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
8cde4af5 46/// \endcond
8cc77df0 47
ea94c18b 48//************* Defaults parameters for reconstruction
49const Bool_t AliMUONTrackReconstructorK::fgkRunSmoother = kTRUE;
8cc77df0 50
8cc77df0 51
7ec3b9cf 52//__________________________________________________________________________
ea94c18b 53AliMUONTrackReconstructorK::AliMUONTrackReconstructorK()
54 : AliMUONVTrackReconstructor()
8cc77df0 55{
ea94c18b 56 /// Constructor for class AliMUONTrackReconstructorK
57 AliInfo("*** Tracking with Kalman Filter ***");
8cc77df0 58}
59
7ec3b9cf 60//__________________________________________________________________________
ea94c18b 61AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
8cc77df0 62{
ea94c18b 63/// Destructor
64}
8cc77df0 65
66 //__________________________________________________________________________
ea94c18b 67void AliMUONTrackReconstructorK::MakeTrackCandidates()
8cc77df0 68{
ea94c18b 69 /// To make track candidates:
70 /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
71 /// Good candidates are made of at least three hitForRec's.
72 /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
208f139e 73 TClonesArray *segments;
ea94c18b 74 AliMUONTrack *track;
75 Int_t iCandidate = 0;
8cc77df0 76
ea94c18b 77 AliDebug(1,"Enter MakeTrackCandidates");
8cc77df0 78
ea94c18b 79 // Loop over stations(1..) 5 and 4 and make track candidates
80 for (Int_t istat=4; istat>=3; istat--) {
81
208f139e 82 // Make segments in the station
83 segments = MakeSegmentsInStation(istat);
ea94c18b 84
85 // Loop over segments
86 for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
87 AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
88
89 // Transform segments to tracks and put them at the end of fRecTracksPtr
90 track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]));
91 fNRecTracks++;
92
93 // Recompute track parameters and covariances on station(1..) 5 using Kalman filter
94 // (to make sure all tracks are treated in the same way)
95 if (istat == 4) RetraceTrack(*track,kFALSE);
96
97 // Printout for debuging
98 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
99 cout<<endl<<"Track parameter covariances at first hit:"<<endl;
100 ((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()))->GetCovariances().Print();
101 }
102
103 // Look for compatible hitForRec(s) in the other station
104 // Remove track if no hit found
105 if (!FollowTrackInStation(*track,7-istat)) {
106 fRecTracksPtr->Remove(track);
107 fNRecTracks--;
108 }
109
110 }
111 // delete the array of segments
7ec3b9cf 112 delete segments;
ea94c18b 113 }
114
115 fRecTracksPtr->Compress(); // this is essential before checking tracks
116
117 // Remove bad tracks
118 Int_t currentNRecTracks = fNRecTracks;
119 for (Int_t iRecTrack = 0; iRecTrack < currentNRecTracks; iRecTrack++) {
120
121 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
122
123 // Remove the track if the normalized chi2 is too high
124 if (track->GetNormalizedChi2() > fgkSigmaToCutForTracking * fgkSigmaToCutForTracking) {
125 fRecTracksPtr->Remove(track);
126 fNRecTracks--;
127 }
128
129 }
130
131 fRecTracksPtr->Compress(); // this is essential before checking tracks
132
133 // Keep all different tracks or only the best ones as required
134 if (fgkTrackAllTracks) RemoveIdenticalTracks();
135 else RemoveDoubleTracks();
136
137 AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
138
8cc77df0 139}
140
ea94c18b 141 //__________________________________________________________________________
142void AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed)
8cc77df0 143{
ea94c18b 144 /// Re-run the kalman filter from the most downstream hit to the most uptream one
145 AliDebug(1,"Enter RetraceTrack");
146
147 AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtHit()->Last();
148
149 // Reset the "seed" (= track parameters and their covariances at last hit) if required
150 if (resetSeed) {
151 // => Shift track parameters at the position of the last hit
152 AliMUONHitForRec* hitForRecAtHit = startingTrackParam->GetHitForRecPtr();
153 startingTrackParam->SetNonBendingCoor(hitForRecAtHit->GetNonBendingCoor());
154 startingTrackParam->SetBendingCoor(hitForRecAtHit->GetBendingCoor());
155
156 // => Re-compute and reset track parameters covariances at last hit (as if the other hits did not exist)
157 const TMatrixD& kParamCov = startingTrackParam->GetCovariances();
158 TMatrixD newParamCov(5,5);
159 newParamCov.Zero();
160 // Non bending plane
161 newParamCov(0,0) = hitForRecAtHit->GetNonBendingReso2();
162 newParamCov(1,1) = 100.*kParamCov(1,1);
163 // Bending plane
164 newParamCov(2,2) = hitForRecAtHit->GetBendingReso2();
165 newParamCov(3,3) = 100.*kParamCov(3,3);
166 // Inverse bending momentum
167 newParamCov(4,4) = 0.5*startingTrackParam->GetInverseBendingMomentum() * 0.5*startingTrackParam->GetInverseBendingMomentum();
168 startingTrackParam->SetCovariances(newParamCov);
169
170 // Reset the track chi2
171 startingTrackParam->SetTrackChi2(0.);
172
173 }
174
175 // Redo the tracking
176 RetracePartialTrack(trackCandidate, startingTrackParam);
177
178}
8cc77df0 179
ea94c18b 180 //__________________________________________________________________________
181void AliMUONTrackReconstructorK::RetracePartialTrack(AliMUONTrack &trackCandidate, const AliMUONTrackParam* startingTrackParam)
182{
183 /// Re-run the kalman filter from the hit attached to startingTrackParam to the most uptream hit
184 AliDebug(1,"Enter RetracePartialTrack");
185
186 // Printout for debuging
187 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
188 cout << "RetracePartialTrack: track chi2 before re-tracking: " << trackCandidate.GetFitFMin() << endl;
189 }
190
191 // Reset the track chi2
192 trackCandidate.SetFitFMin(startingTrackParam->GetTrackChi2());
193
194 // loop over attached hits until the first one and recompute track parameters and covariances using kalman filter
195 Int_t expectedChamber = startingTrackParam->GetHitForRecPtr()->GetChamberNumber() - 1;
196 Int_t currentChamber;
197 Double_t addChi2TrackAtHit;
198 AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtHit()->Before(startingTrackParam);
199 while (trackParamAtHit) {
200
201 // reset track parameters and their covariances
202 trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
203 trackParamAtHit->SetZ(startingTrackParam->GetZ());
204 trackParamAtHit->SetCovariances(startingTrackParam->GetCovariances());
205
206 // add MCS effect
207 AliMUONTrackExtrap::AddMCSEffect(trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
208
209 // reset propagator for smoother
210 if (fgkRunSmoother) trackParamAtHit->ResetPropagator();
211
212 // add MCS in missing chambers if any (at most 2 chambers can be missing according to tracking criteria)
213 currentChamber = trackParamAtHit->GetHitForRecPtr()->GetChamberNumber();
214 while (currentChamber < expectedChamber) {
215 // extrapolation to the missing chamber (update the propagator)
216 AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, AliMUONConstants::DefaultChamberZ(expectedChamber), fgkRunSmoother);
217 // add MCS effect
218 AliMUONTrackExtrap::AddMCSEffect(trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
219 expectedChamber--;
220 }
221
222 // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit (update the propagator)
223 AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, trackParamAtHit->GetHitForRecPtr()->GetZ(), fgkRunSmoother);
224
225 if (fgkRunSmoother) {
226 // save extrapolated parameters for smoother
227 trackParamAtHit->SetExtrapParameters(trackParamAtHit->GetParameters());
228
229 // save extrapolated covariance matrix for smoother
230 trackParamAtHit->SetExtrapCovariances(trackParamAtHit->GetCovariances());
231 }
232
233 // Compute new track parameters including "hitForRecCh2" using kalman filter
234 addChi2TrackAtHit = RunKalmanFilter(*trackParamAtHit);
235
236 // Update the track chi2
237 trackCandidate.SetFitFMin(trackCandidate.GetFitFMin() + addChi2TrackAtHit);
238 trackParamAtHit->SetTrackChi2(trackCandidate.GetFitFMin());
239
240 // prepare next step, add MCS effects in parameter covariances
241 expectedChamber--;
242 startingTrackParam = trackParamAtHit;
243 trackParamAtHit = (AliMUONTrackParam*) (trackCandidate.GetTrackParamAtHit()->Before(startingTrackParam));
244 }
245
246 // Printout for debuging
247 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
248 cout << "RetracePartialTrack: track chi2 after re-tracking: " << trackCandidate.GetFitFMin() << endl;
249 }
250
251}
8cc77df0 252
ea94c18b 253 //__________________________________________________________________________
254void AliMUONTrackReconstructorK::FollowTracks()
255{
256 /// Follow tracks in stations(1..) 3, 2 and 1
257 AliDebug(1,"Enter FollowTracks");
258
259 AliMUONTrack *track;
260 Int_t currentNRecTracks;
261 Bool_t hitFound;
262
263 for (Int_t station = 2; station >= 0; station--) {
264
265 // Save the actual number of reconstructed track in case of
266 // tracks are added or suppressed during the tracking procedure
267 // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
268 currentNRecTracks = fNRecTracks;
269
270 for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
271 AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
272
273 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
274
275 // Printout for debuging
276 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
277 cout<<endl<<"Track parameter covariances at first hit:"<<endl;
278 ((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()))->GetCovariances().Print();
279 }
280
281 // Look for compatible hitForRec in station(0..) "station"
282 hitFound = FollowTrackInStation(*track,station);
283
284 // Try to recover track if required
285 if (!hitFound && fgkRecoverTracks) hitFound = RecoverTrack(*track,station);
286
287 // remove track if no hit found
288 if (!hitFound) {
289 fRecTracksPtr->Remove(track);
290 fNRecTracks--;
8cc77df0 291 }
ea94c18b 292
8cc77df0 293 }
ea94c18b 294
295 fRecTracksPtr->Compress(); // this is essential before checking tracks
296
297 // Remove bad tracks
298 Int_t currentNRecTracks = fNRecTracks;
299 for (Int_t iRecTrack = 0; iRecTrack < currentNRecTracks; iRecTrack++) {
300 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
301
302 // Remove the track if the normalized chi2 is too high
303 if (track->GetNormalizedChi2() > fgkSigmaToCutForTracking * fgkSigmaToCutForTracking) {
304 fRecTracksPtr->Remove(track);
305 fNRecTracks--;
306 }
307
308 }
309
310 fRecTracksPtr->Compress(); // this is essential before checking tracks
311
312 // Keep only the best tracks if required
313 if (!fgkTrackAllTracks) RemoveDoubleTracks();
314
315 }
316
317}
b0e4c598 318
ea94c18b 319 //__________________________________________________________________________
320Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandidate, Int_t nextStation)
321{
322 /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s)
323 /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
324 /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
325 /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
326 /// Remove the obsolete "trackCandidate" at the end.
327 /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority.
328 /// return kTRUE if new hits have been found (otherwise return kFALSE)
329 AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
330
331 // Order the chamber according to the propagation direction (tracking starts with chamber 2):
332 // - nextStation == station(1...) 5 => forward propagation
333 // - nextStation < station(1...) 5 => backward propagation
334 Int_t ch1, ch2;
335 if (nextStation==4) {
336 ch1 = 2*nextStation+1;
337 ch2 = 2*nextStation;
338 } else {
339 ch1 = 2*nextStation;
340 ch2 = 2*nextStation+1;
341 }
342
343 Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2);
344 Double_t chi2OfHitForRec;
345 Double_t maxChi2OfHitForRec = 2. * fgkSigmaToCutForTracking * fgkSigmaToCutForTracking; // 2 because 2 quantities in chi2
346 Double_t addChi2TrackAtHit1;
347 Double_t addChi2TrackAtHit2;
348 Double_t bestAddChi2TrackAtHit1 = 1.e10;
349 Double_t bestAddChi2TrackAtHit2 = 1.e10;
350 Bool_t foundOneHit = kFALSE;
351 Bool_t foundTwoHits = kFALSE;
352 AliMUONTrack *newTrack = 0x0;
353 AliMUONHitForRec *hitForRecCh1, *hitForRecCh2;
354 AliMUONTrackParam extrapTrackParam;
355 AliMUONTrackParam extrapTrackParamAtHit1;
356 AliMUONTrackParam extrapTrackParamAtHit2;
357 AliMUONTrackParam bestTrackParamAtHit1;
358 AliMUONTrackParam bestTrackParamAtHit2;
359 Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]];
360 for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE;
361
362 // Get track parameters and extrapolate them to chamber "ch2" to save computing time in the next steps
363 AliMUONTrackParam extrapTrackParamAtCh2(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->First());
364
365 // Add MCS effect
366 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.);
367
368 // reset propagator for smoother
369 if (fgkRunSmoother) extrapTrackParamAtCh2.ResetPropagator();
370
371 // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria)
372 if (ch1 < ch2 && extrapTrackParamAtCh2.GetHitForRecPtr()->GetChamberNumber() > ch2 + 1) {
373 // extrapolation to the missing chamber
374 AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh2, AliMUONConstants::DefaultChamberZ(ch2 + 1), fgkRunSmoother);
375 // add MCS effect
376 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.);
377 }
378
379 //Extrapolate trackCandidate to chamber "ch2"
380 AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh2, zCh2, fgkRunSmoother);
381
382 // Printout for debuging
383 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
384 cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl;
385 extrapTrackParamAtCh2.GetCovariances().Print();
386 }
387
388 // Printout for debuging
389 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
390 cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl;
391 }
392
393 // look for candidates in chamber 2
394 for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) {
395
396 hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2);
397
398 // try to add the current hit
399 chi2OfHitForRec = TryOneHitForRec(extrapTrackParamAtCh2, hitForRecCh2, extrapTrackParamAtHit2, fgkRunSmoother);
400
401 // if good chi2 then try to attach a hitForRec in the other chamber too
402 if (chi2OfHitForRec < maxChi2OfHitForRec) {
403 Bool_t foundSecondHit = kFALSE;
404
405 // Printout for debuging
406 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
407 cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch2+1
408 << " (Chi2 = " << chi2OfHitForRec << ")" << endl;
409 cout << " look for second hits in chamber(1..): " << ch1+1 << " ..." << endl;
8cc77df0 410 }
ea94c18b 411
412 if (fgkRunSmoother) {
413 // save extrapolated parameters for smoother
414 extrapTrackParamAtHit2.SetExtrapParameters(extrapTrackParamAtHit2.GetParameters());
415
416 // save extrapolated covariance matrix for smoother
417 extrapTrackParamAtHit2.SetExtrapCovariances(extrapTrackParamAtHit2.GetCovariances());
418 }
419
420 // Compute new track parameters including "hitForRecCh2" using kalman filter
421 addChi2TrackAtHit2 = RunKalmanFilter(extrapTrackParamAtHit2);
422
423 // copy new track parameters for next step
424 extrapTrackParam = extrapTrackParamAtHit2;
425
426 // add MCS effect
427 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
428
429 // reset propagator for smoother
430 if (fgkRunSmoother) extrapTrackParam.ResetPropagator();
431
432 for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
433
434 hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
435
436 // try to add the current hit
437 chi2OfHitForRec = TryOneHitForRec(extrapTrackParam, hitForRecCh1, extrapTrackParamAtHit1, fgkRunSmoother);
438
439 // if good chi2 then consider to add the 2 hitForRec to the "trackCandidate"
440 if (chi2OfHitForRec < maxChi2OfHitForRec) {
441 foundSecondHit = kTRUE;
442 foundTwoHits = kTRUE;
443
444 // Printout for debuging
445 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
446 cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch1+1
447 << " (Chi2 = " << chi2OfHitForRec << ")" << endl;
448 }
449
450 if (fgkRunSmoother) {
451 // save extrapolated parameters for smoother
452 extrapTrackParamAtHit1.SetExtrapParameters(extrapTrackParamAtHit1.GetParameters());
453
454 // save extrapolated covariance matrix for smoother
455 extrapTrackParamAtHit1.SetExtrapCovariances(extrapTrackParamAtHit1.GetCovariances());
456 }
457
458 // Compute new track parameters including "hitForRecCh1" using kalman filter
459 addChi2TrackAtHit1 = RunKalmanFilter(extrapTrackParamAtHit1);
460
461 if (fgkTrackAllTracks) {
462 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
463 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
464 UpdateTrack(*newTrack,extrapTrackParamAtHit1,extrapTrackParamAtHit2,addChi2TrackAtHit1,addChi2TrackAtHit2);
465 fNRecTracks++;
466
467 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
468 // (going in the right direction)
469 if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
470
471 // Tag hitForRecCh1 as used
472 hitForRecCh1Used[hit1] = kTRUE;
473
474 // Printout for debuging
475 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
476 cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1 << endl;
477 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
478 }
479
480 } else if (addChi2TrackAtHit1+addChi2TrackAtHit2 < bestAddChi2TrackAtHit1+bestAddChi2TrackAtHit2) {
481 // keep track of the best couple of hits
482 bestAddChi2TrackAtHit1 = addChi2TrackAtHit1;
483 bestAddChi2TrackAtHit2 = addChi2TrackAtHit2;
484 bestTrackParamAtHit1 = extrapTrackParamAtHit1;
485 bestTrackParamAtHit2 = extrapTrackParamAtHit2;
486 }
487
b0e4c598 488 }
ea94c18b 489
8cc77df0 490 }
ea94c18b 491
492 // if no hitForRecCh1 found then consider to add hitForRecCh2 only
493 if (!foundSecondHit) {
494 foundOneHit = kTRUE;
495
496 if (fgkTrackAllTracks) {
497 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
498 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
499 UpdateTrack(*newTrack,extrapTrackParamAtHit2,addChi2TrackAtHit2);
500 fNRecTracks++;
501
502 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
503 // (going in the right direction)
504 if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
505
506 // Printout for debuging
507 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
508 cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1 << endl;
509 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
510 }
511
512 } else if (!foundTwoHits && addChi2TrackAtHit2 < bestAddChi2TrackAtHit1) {
513 // keep track of the best single hitForRec except if a couple
514 // of hits has already been found (i.e. bestHitForRec2!=0x0)
515 bestAddChi2TrackAtHit1 = addChi2TrackAtHit2;
516 bestTrackParamAtHit1 = extrapTrackParamAtHit2;
517 }
518
519 }
520
8cc77df0 521 }
ea94c18b 522
523 }
524
525 // look for candidates in chamber 1 not already attached to a track
526 // if we want to keep all possible tracks or if no good couple of hitForRec has been found
527 if (fgkTrackAllTracks || !foundTwoHits) {
528
529 // Printout for debuging
530 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
531 cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl;
8cc77df0 532 }
ea94c18b 533
534 // add MCS effect for next step
535 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.);
536
537 for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
538
539 hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
540
541 if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used
542
543 // try to add the current hit
544 chi2OfHitForRec = TryOneHitForRec(extrapTrackParamAtCh2, hitForRecCh1, extrapTrackParamAtHit1, fgkRunSmoother);
545
546 // if good chi2 then consider to add hitForRecCh1
547 // We do not try to attach a hitForRec in the other chamber too since it has already been done above
548 if (chi2OfHitForRec < maxChi2OfHitForRec) {
549 foundOneHit = kTRUE;
550
551 // Printout for debuging
552 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
553 cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch1+1
554 << " (Chi2 = " << chi2OfHitForRec << ")" << endl;
555 }
556
557 if (fgkRunSmoother) {
558 // save extrapolated parameters for smoother
559 extrapTrackParamAtHit1.SetExtrapParameters(extrapTrackParamAtHit1.GetParameters());
560
561 // save extrapolated covariance matrix for smoother
562 extrapTrackParamAtHit1.SetExtrapCovariances(extrapTrackParamAtHit1.GetCovariances());
563 }
564
565 // Compute new track parameters including "hitForRecCh1" using kalman filter
566 addChi2TrackAtHit1 = RunKalmanFilter(extrapTrackParamAtHit1);
567
568 if (fgkTrackAllTracks) {
569 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
570 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
571 UpdateTrack(*newTrack,extrapTrackParamAtHit1,addChi2TrackAtHit1);
572 fNRecTracks++;
573
574 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
575 // (going in the right direction)
576 if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
577
578 // Printout for debuging
579 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
580 cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1 << endl;
581 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
582 }
583
584 } else if (addChi2TrackAtHit1 < bestAddChi2TrackAtHit1) {
585 // keep track of the best single hitForRec except if a couple
586 // of hits has already been found (i.e. bestHitForRec1!=0x0)
587 bestAddChi2TrackAtHit1 = addChi2TrackAtHit1;
588 bestTrackParamAtHit1 = extrapTrackParamAtHit1;
589 }
590
591 }
592
8cc77df0 593 }
ea94c18b 594
595 }
596
597 // fill out the best track if required else clean up the fRecTracksPtr array
598 if (!fgkTrackAllTracks) {
599 if (foundTwoHits) {
600 UpdateTrack(trackCandidate,bestTrackParamAtHit1,bestTrackParamAtHit2,bestAddChi2TrackAtHit1,bestAddChi2TrackAtHit2);
601
602 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
603 // (going in the right direction)
604 if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE);
8cc77df0 605
ea94c18b 606 // Printout for debuging
607 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
608 cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1 << endl;
609 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
610 }
611
612 } else if (foundOneHit) {
613 UpdateTrack(trackCandidate,bestTrackParamAtHit1,bestAddChi2TrackAtHit1);
614
615 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
616 // (going in the right direction)
617 if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE);
618
619 // Printout for debuging
620 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
621 cout << "FollowTrackInStation: added the best hit in chamber(1..): " << bestTrackParamAtHit1.GetHitForRecPtr()->GetChamberNumber()+1 << endl;
622 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
623 }
624
625 } else return kFALSE;
626
627 } else if (foundOneHit || foundTwoHits) {
628
629 // remove obsolete track
630 fRecTracksPtr->Remove(&trackCandidate);
631 fNRecTracks--;
632
633 } else return kFALSE;
634
635 return kTRUE;
636
637}
8cc77df0 638
ea94c18b 639 //__________________________________________________________________________
640Double_t AliMUONTrackReconstructorK::RunKalmanFilter(AliMUONTrackParam &trackParamAtHit)
641{
642 /// Compute new track parameters and their covariances including new hit using kalman filter
643 /// return the additional track chi2
644 AliDebug(1,"Enter RunKalmanFilter");
645
646 // Get actual track parameters (p)
647 TMatrixD param(trackParamAtHit.GetParameters());
648
649 // Get new hit parameters (m)
650 AliMUONHitForRec *hitForRecAtHit = trackParamAtHit.GetHitForRecPtr();
651 TMatrixD hit(5,1);
652 hit.Zero();
653 hit(0,0) = hitForRecAtHit->GetNonBendingCoor();
654 hit(2,0) = hitForRecAtHit->GetBendingCoor();
655
656 // Compute the actual parameter weight (W)
657 TMatrixD paramWeight(trackParamAtHit.GetCovariances());
658 if (paramWeight.Determinant() != 0) {
659 paramWeight.Invert();
660 } else {
661 AliWarning(" Determinant = 0");
662 return 1.e10;
8cc77df0 663 }
ea94c18b 664
665 // Compute the new hit weight (U)
666 TMatrixD hitWeight(5,5);
667 hitWeight.Zero();
668 hitWeight(0,0) = 1. / hitForRecAtHit->GetNonBendingReso2();
669 hitWeight(2,2) = 1. / hitForRecAtHit->GetBendingReso2();
8cc77df0 670
ea94c18b 671 // Compute the new parameters covariance matrix ( (W+U)^-1 )
672 TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,hitWeight);
673 if (newParamCov.Determinant() != 0) {
674 newParamCov.Invert();
675 } else {
676 AliWarning(" Determinant = 0");
677 return 1.e10;
678 }
679
680 // Save the new parameters covariance matrix
681 trackParamAtHit.SetCovariances(newParamCov);
682
683 // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
684 TMatrixD tmp(hit,TMatrixD::kMinus,param);
685 TMatrixD tmp2(hitWeight,TMatrixD::kMult,tmp); // U(m-p)
686 TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
687 newParam += param; // ((W+U)^-1)U(m-p) + p
688
689 // Save the new parameters
690 trackParamAtHit.SetParameters(newParam);
691
692 // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
693 tmp = newParam; // p'
694 tmp -= param; // (p'-p)
695 TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
696 TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
697 tmp = newParam; // p'
698 tmp -= hit; // (p'-m)
699 TMatrixD tmp4(hitWeight,TMatrixD::kMult,tmp); // U(p'-m)
700 addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
701
702 return addChi2Track(0,0);
703
8cc77df0 704}
705
ea94c18b 706 //__________________________________________________________________________
707void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit, Double_t addChi2)
8cc77df0 708{
ea94c18b 709 /// Add 1 hit to the track candidate
710 /// Update chi2 of the track
711
712 // Update the track chi2 into TrackParamAtHit
713 trackParamAtHit.SetTrackChi2(track.GetFitFMin() + addChi2);
714
715 // Update the chi2 of the new track
716 track.SetFitFMin(trackParamAtHit.GetTrackChi2());
717
718 // Update array of TrackParamAtHit
719 track.AddTrackParamAtHit(&trackParamAtHit,trackParamAtHit.GetHitForRecPtr());
720 track.GetTrackParamAtHit()->Sort();
721
722}
8cc77df0 723
ea94c18b 724 //__________________________________________________________________________
725void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit1, AliMUONTrackParam &trackParamAtHit2,
726 Double_t addChi2AtHit1, Double_t addChi2AtHit2)
727{
728 /// Add 2 hits to the track candidate
729 /// Update track and local chi2
730
731 // Update local chi2 at first hit
732 AliMUONHitForRec* hit1 = trackParamAtHit1.GetHitForRecPtr();
733 Double_t deltaX = trackParamAtHit1.GetNonBendingCoor() - hit1->GetNonBendingCoor();
734 Double_t deltaY = trackParamAtHit1.GetBendingCoor() - hit1->GetBendingCoor();
735 Double_t localChi2AtHit1 = deltaX*deltaX / hit1->GetNonBendingReso2() +
736 deltaY*deltaY / hit1->GetBendingReso2();
737 trackParamAtHit1.SetLocalChi2(localChi2AtHit1);
738
739 // Flag first hit as being removable
740 trackParamAtHit1.SetRemovable(kTRUE);
741
742 // Update local chi2 at second hit
743 AliMUONHitForRec* hit2 = trackParamAtHit2.GetHitForRecPtr();
744 AliMUONTrackParam extrapTrackParamAtHit2(trackParamAtHit1);
745 AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParamAtHit2, trackParamAtHit2.GetZ());
746 deltaX = extrapTrackParamAtHit2.GetNonBendingCoor() - hit2->GetNonBendingCoor();
747 deltaY = extrapTrackParamAtHit2.GetBendingCoor() - hit2->GetBendingCoor();
748 Double_t localChi2AtHit2 = deltaX*deltaX / hit2->GetNonBendingReso2() +
749 deltaY*deltaY / hit2->GetBendingReso2();
750 trackParamAtHit2.SetLocalChi2(localChi2AtHit2);
751
752 // Flag second hit as being removable
753 trackParamAtHit2.SetRemovable(kTRUE);
754
755 // Update the track chi2 into TrackParamAtHit1
756 trackParamAtHit1.SetTrackChi2(track.GetFitFMin() + addChi2AtHit1);
757
758 // Update the track chi2 into TrackParamAtHit2
759 trackParamAtHit2.SetTrackChi2(trackParamAtHit1.GetTrackChi2() + addChi2AtHit2);
760
761 // Update the chi2 of the new track
762 track.SetFitFMin(trackParamAtHit2.GetTrackChi2());
763
764 // Update array of TrackParamAtHit
765 track.AddTrackParamAtHit(&trackParamAtHit1,trackParamAtHit1.GetHitForRecPtr());
766 track.AddTrackParamAtHit(&trackParamAtHit2,trackParamAtHit2.GetHitForRecPtr());
767 track.GetTrackParamAtHit()->Sort();
768
769}
8cc77df0 770
ea94c18b 771 //__________________________________________________________________________
772Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, Int_t nextStation)
773{
774 /// Try to recover the track candidate in the next station
775 /// by removing the worst of the two hits attached in the current station
776 /// Return kTRUE if recovering succeeds
777 AliDebug(1,"Enter RecoverTrack");
778
779 // Do not try to recover track until we have attached hit(s) on station(1..) 3
780 if (nextStation > 1) return kFALSE;
781
782 Int_t worstHitNumber = -1;
783 Double_t localChi2, worstChi2 = 0.;
784
785 // Look for the hit to remove
786 for (Int_t hitNumber = 0; hitNumber < 2; hitNumber++) {
787 AliMUONTrackParam *trackParamAtHit = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(hitNumber);
788
789 // check if current hit is removable
790 if (!trackParamAtHit->IsRemovable()) return kFALSE;
791
792 // Pick up hit with the worst chi2
793 localChi2 = trackParamAtHit->GetLocalChi2();
794 if (localChi2 > worstChi2) {
795 worstChi2 = localChi2;
796 worstHitNumber = hitNumber;
8cc77df0 797 }
ea94c18b 798 }
799
800 // Reset best hit as being NOT removable
801 ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt((worstHitNumber+1)%2))->SetRemovable(kFALSE);
802
803 // Remove the worst hit
804 trackCandidate.RemoveTrackParamAtHit((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(worstHitNumber));
805
806 // Re-calculate track parameters at the (new) first hit
807 RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(1));
808
809 // Look for new hit(s) in next station
810 return FollowTrackInStation(trackCandidate,nextStation);
811
8cc77df0 812}
813
814 //__________________________________________________________________________
ea94c18b 815Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track)
8cc77df0 816{
ea94c18b 817 /// Compute new track parameters and their covariances using smoother
818 AliDebug(1,"Enter RunSmoother");
819
820 AliMUONTrackParam *previousTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtHit()->First();
821
822 // Smoothed parameters and covariances at first hit = filtered parameters and covariances
823 previousTrackParam->SetSmoothParameters(previousTrackParam->GetParameters());
824 previousTrackParam->SetSmoothCovariances(previousTrackParam->GetCovariances());
825
826 // Compute local chi2 at first hit
827 AliMUONHitForRec *hitForRecAtHit = previousTrackParam->GetHitForRecPtr();
828 Double_t dX = hitForRecAtHit->GetNonBendingCoor() - previousTrackParam->GetNonBendingCoor();
829 Double_t dY = hitForRecAtHit->GetBendingCoor() - previousTrackParam->GetBendingCoor();
830 Double_t localChi2 = dX * dX / hitForRecAtHit->GetNonBendingReso2() + dY * dY / hitForRecAtHit->GetBendingReso2();
831
832 // Save local chi2 at first hit
833 previousTrackParam->SetLocalChi2(localChi2);
834
835 AliMUONTrackParam *currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtHit()->After(previousTrackParam);
836 while (currentTrackParam) {
837
838 // Get variables
839 const TMatrixD &extrapParameters = previousTrackParam->GetExtrapParameters(); // X(k+1 k)
840 const TMatrixD &filteredParameters = currentTrackParam->GetParameters(); // X(k k)
841 const TMatrixD &previousSmoothParameters = previousTrackParam->GetSmoothParameters(); // X(k+1 n)
842 const TMatrixD &propagator = previousTrackParam->GetPropagator(); // F(k)
843 const TMatrixD &extrapCovariances = previousTrackParam->GetExtrapCovariances(); // C(k+1 k)
844 const TMatrixD &filteredCovariances = currentTrackParam->GetCovariances(); // C(k k)
845 const TMatrixD &previousSmoothCovariances = previousTrackParam->GetSmoothCovariances(); // C(k+1 n)
846
847 // Compute smoother gain: A(k) = C(kk) * F(f)^t * (C(k+1 k))^-1
848 TMatrixD extrapWeight(extrapCovariances);
849 if (extrapWeight.Determinant() != 0) {
850 extrapWeight.Invert(); // (C(k+1 k))^-1
851 } else {
852 AliWarning(" Determinant = 0");
853 return kFALSE;
854 }
855 TMatrixD smootherGain(filteredCovariances,TMatrixD::kMultTranspose,propagator); // C(kk) * F(f)^t
856 smootherGain *= extrapWeight; // C(kk) * F(f)^t * (C(k+1 k))^-1
857
858 // Compute smoothed parameters: X(k n) = X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
859 TMatrixD tmpParam(previousSmoothParameters,TMatrixD::kMinus,extrapParameters); // X(k+1 n) - X(k+1 k)
860 TMatrixD smoothParameters(smootherGain,TMatrixD::kMult,tmpParam); // A(k) * (X(k+1 n) - X(k+1 k))
861 smoothParameters += filteredParameters; // X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
862
863 // Save smoothed parameters
864 currentTrackParam->SetSmoothParameters(smoothParameters);
865
866 // Compute smoothed covariances: C(k n) = C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
867 TMatrixD tmpCov(previousSmoothCovariances,TMatrixD::kMinus,extrapCovariances); // C(k+1 n) - C(k+1 k)
868 TMatrixD tmpCov2(tmpCov,TMatrixD::kMultTranspose,smootherGain); // (C(k+1 n) - C(k+1 k)) * (A(k))^t
869 TMatrixD smoothCovariances(smootherGain,TMatrixD::kMult,tmpCov2); // A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
870 smoothCovariances += filteredCovariances; // C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
871
872 // Save smoothed covariances
873 currentTrackParam->SetSmoothCovariances(smoothCovariances);
874
875 // Compute smoothed residual: r(k n) = hit - X(k n)
876 hitForRecAtHit = currentTrackParam->GetHitForRecPtr();
877 TMatrixD smoothResidual(2,1);
878 smoothResidual.Zero();
879 smoothResidual(0,0) = hitForRecAtHit->GetNonBendingCoor() - smoothParameters(0,0);
880 smoothResidual(1,0) = hitForRecAtHit->GetBendingCoor() - smoothParameters(2,0);
881
882 // Compute weight of smoothed residual: W(k n) = (hitCov - C(k n))^-1
883 TMatrixD smoothResidualWeight(2,2);
884 smoothResidualWeight(0,0) = hitForRecAtHit->GetNonBendingReso2() - smoothCovariances(0,0);
885 smoothResidualWeight(0,1) = - smoothCovariances(0,2);
886 smoothResidualWeight(1,0) = - smoothCovariances(2,0);
887 smoothResidualWeight(1,1) = hitForRecAtHit->GetBendingReso2() - smoothCovariances(2,2);
888 if (smoothResidualWeight.Determinant() != 0) {
889 smoothResidualWeight.Invert();
890 } else {
891 AliWarning(" Determinant = 0");
892 return kFALSE;
893 }
894
895 // Compute local chi2 = (r(k n))^t * W(k n) * r(k n)
896 TMatrixD tmpChi2(smoothResidual,TMatrixD::kTransposeMult,smoothResidualWeight); // (r(k n))^t * W(k n)
897 TMatrixD localChi2(tmpChi2,TMatrixD::kMult,smoothResidual); // (r(k n))^t * W(k n) * r(k n)
898
899 // Save local chi2
900 currentTrackParam->SetLocalChi2(localChi2(0,0));
901
902 previousTrackParam = currentTrackParam;
903 currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtHit()->After(previousTrackParam);
904 }
905
906 return kTRUE;
907
8cc77df0 908}
909
8cc77df0 910 //__________________________________________________________________________
ea94c18b 911void AliMUONTrackReconstructorK::ImproveTracks()
8cc77df0 912{
ea94c18b 913 /// Improve tracks by removing clusters with local chi2 highter than the defined cut
914 /// Recompute track parameters and covariances at the remaining clusters
915 AliDebug(1,"Enter ImproveTracks");
916
917 Double_t localChi2, worstLocalChi2;
918 Int_t worstChamber;
919 AliMUONTrackParam *trackParamAtHit, *worstTrackParamAtHit;
920 Bool_t smoothed;
921
922 // Remove double track to improve only "good" tracks
923 RemoveDoubleTracks();
924
925 AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
8cc77df0 926 while (track) {
ea94c18b 927
928 while (!track->IsImproved()) {
929
930 // Run smoother if required
931 smoothed = kFALSE;
932 if (fgkRunSmoother) smoothed = RunSmoother(*track);
933
934 // Use standard procedure to compute local chi2 if smoother not required or not working
935 if (!smoothed) {
936
937 // Update track parameters and covariances
938 track->UpdateCovTrackParamAtHit();
939
940 // Compute local chi2 of each hits
941 track->ComputeLocalChi2(kTRUE);
942 }
943
944 // Look for the hit to remove
945 worstTrackParamAtHit = NULL;
946 worstLocalChi2 = 0.;
947 trackParamAtHit = (AliMUONTrackParam*)track->GetTrackParamAtHit()->First();
948 while (trackParamAtHit) {
949
950 // Pick up hit with the worst chi2
951 localChi2 = trackParamAtHit->GetLocalChi2();
952 if (localChi2 > worstLocalChi2) {
953 worstLocalChi2 = localChi2;
954 worstTrackParamAtHit = trackParamAtHit;
955 }
956
957 trackParamAtHit = (AliMUONTrackParam*)track->GetTrackParamAtHit()->After(trackParamAtHit);
958 }
959
960 // Check if bad removable hit found
961 if (!worstTrackParamAtHit) {
962 track->SetImproved(kTRUE);
963 break;
964 }
965
966 // check whether the worst hit is removable or not
967 if (!worstTrackParamAtHit->IsRemovable()) {
968 track->SetImproved(kTRUE);
969 break;
970 }
971
972 // Check whether the worst chi2 is under requirement or not
973 if (worstLocalChi2 < 2. * fgkSigmaToCutForImprovement * fgkSigmaToCutForImprovement) { // 2 because 2 quantities in chi2
974 track->SetImproved(kTRUE);
975 break;
976 }
977
978 // Reset the second hit in the same station as the bad one as being NOT removable
979 worstChamber = worstTrackParamAtHit->GetHitForRecPtr()->GetChamberNumber();
980 if (worstChamber%2 == 0) ((AliMUONTrackParam*)track->GetTrackParamAtHit()->After(worstTrackParamAtHit))->SetRemovable(kFALSE);
981 else ((AliMUONTrackParam*)track->GetTrackParamAtHit()->Before(worstTrackParamAtHit))->SetRemovable(kFALSE);
982
983 // Save pointer to the trackParamAtHit next to the one to be removed
984 trackParamAtHit = (AliMUONTrackParam*)track->GetTrackParamAtHit()->After(worstTrackParamAtHit);
985
986 // Remove the worst hit
987 track->RemoveTrackParamAtHit(worstTrackParamAtHit);
988
989 // Re-calculate track parameters
990 // - from the hit immediately downstream the one suppressed
991 // - or from the begining - if parameters have been re-computed using the standard method (kalman parameters have been lost)
992 // - or if the removed hit was the last one
993 if (smoothed && trackParamAtHit) RetracePartialTrack(*track,trackParamAtHit);
994 else RetraceTrack(*track,kTRUE);
995
996 // Printout for debuging
997 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
998 cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(track)+1 << " improved " << endl;
999 }
1000
1001 }
1002
1003 track = (AliMUONTrack*) fRecTracksPtr->After(track);
1004 }
1005
8cc77df0 1006}
1007
b0e4c598 1008 //__________________________________________________________________________
ea94c18b 1009void AliMUONTrackReconstructorK::Finalize()
b0e4c598 1010{
ea94c18b 1011 /// Fill AliMUONTrack's fHitForRecAtHit array
1012 AliMUONTrack *track;
1013 AliMUONTrackParam *trackParamAtHit;
1014 Bool_t smoothed = kFALSE;
1015
1016 track = (AliMUONTrack*) fRecTracksPtr->First();
1017 while (track) {
1018
1019 // update track parameters (using smoother if required) if not already done
1020 if (!track->IsImproved()) {
1021 smoothed = kFALSE;
1022 if (fgkRunSmoother) smoothed = RunSmoother(*track);
1023 if (!smoothed) track->UpdateCovTrackParamAtHit();
1024 }
1025
1026 trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
1027 while (trackParamAtHit) {
1028
1029 // copy smoothed parameters and covariances if any
1030 if (smoothed) {
1031 trackParamAtHit->SetParameters(trackParamAtHit->GetSmoothParameters());
1032 trackParamAtHit->SetCovariances(trackParamAtHit->GetSmoothCovariances());
1033 }
1034
1035 // update array of track hits
1036 track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
1037
1038 trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
1039 }
1040
1041 track = (AliMUONTrack*) fRecTracksPtr->After(track);
1042
b0e4c598 1043 }
ea94c18b 1044
b0e4c598 1045}
ea94c18b 1046