Main changes:
[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"
35#include "AliMUONVClusterStore.h"
ea94c18b 36#include "AliMUONTrack.h"
37#include "AliMUONTrackParam.h"
38#include "AliMUONTrackExtrap.h"
8cde4af5 39
8cc77df0 40#include "AliLog.h"
41
8cde4af5 42#include <Riostream.h>
ea94c18b 43#include <TMath.h>
44#include <TMatrixD.h>
8cde4af5 45
46/// \cond CLASSIMP
8cc77df0 47ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
8cde4af5 48/// \endcond
8cc77df0 49
019df241 50 //__________________________________________________________________________
ea94c18b 51AliMUONTrackReconstructorK::AliMUONTrackReconstructorK()
52 : AliMUONVTrackReconstructor()
8cc77df0 53{
96ebe67e 54 /// Constructor
8cc77df0 55}
56
019df241 57 //__________________________________________________________________________
ea94c18b 58AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK()
8cc77df0 59{
ea94c18b 60/// Destructor
61}
8cc77df0 62
63 //__________________________________________________________________________
96ebe67e 64void AliMUONTrackReconstructorK::MakeTrackCandidates(const AliMUONVClusterStore& clusterStore)
8cc77df0 65{
019df241 66 /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
ea94c18b 67 /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
96ebe67e 68 /// Good candidates are made of at least three clusters.
ea94c18b 69 /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
b709ac13 70
208f139e 71 TClonesArray *segments;
ea94c18b 72 AliMUONTrack *track;
73 Int_t iCandidate = 0;
96ebe67e 74 Bool_t clusterFound;
8cc77df0 75
ea94c18b 76 AliDebug(1,"Enter MakeTrackCandidates");
8cc77df0 77
ea94c18b 78 // Loop over stations(1..) 5 and 4 and make track candidates
79 for (Int_t istat=4; istat>=3; istat--) {
80
208f139e 81 // Make segments in the station
96ebe67e 82 segments = MakeSegmentsInStation(clusterStore, istat);
ea94c18b 83
84 // Loop over segments
85 for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
86 AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
87
88 // Transform segments to tracks and put them at the end of fRecTracksPtr
89 track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]));
90 fNRecTracks++;
91
ea94c18b 92 // Printout for debuging
93 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
96ebe67e 94 cout<<endl<<"Track parameter covariances at first cluster:"<<endl;
95 ((AliMUONTrackParam*) (track->GetTrackParamAtCluster()->First()))->GetCovariances().Print();
ea94c18b 96 }
97
96ebe67e 98 // Look for compatible cluster(s) in the other station
99 if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast())
100 clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
019df241 101 else {
102 // First recompute track parameters and covariances on station(1..) 5 using Kalman filter
103 // (to make sure all tracks are treated in the same way)
96ebe67e 104 if (istat == 4) RetraceTrack(*track, kFALSE);
105 clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
019df241 106 }
107
96ebe67e 108 // Remove track if no cluster found
109 if (!clusterFound) {
ea94c18b 110 fRecTracksPtr->Remove(track);
111 fNRecTracks--;
112 }
113
114 }
115 // delete the array of segments
7ec3b9cf 116 delete segments;
ea94c18b 117 }
118
ea94c18b 119
019df241 120 // Retrace tracks using Kalman filter and remove bad ones
3304fa09 121 if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast()) {
019df241 122 fRecTracksPtr->Compress(); // this is essential before checking tracks
ea94c18b 123
019df241 124 Int_t currentNRecTracks = fNRecTracks;
125 for (Int_t iRecTrack = 0; iRecTrack < currentNRecTracks; iRecTrack++) {
126 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
127
128 // Recompute track parameters and covariances using Kalman filter
129 RetraceTrack(*track,kTRUE);
130
131 // Remove the track if the normalized chi2 is too high
3304fa09 132 if (track->GetNormalizedChi2() > AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
133 AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking()) {
019df241 134 fRecTracksPtr->Remove(track);
135 fNRecTracks--;
136 }
137
ea94c18b 138 }
139
140 }
141
142 fRecTracksPtr->Compress(); // this is essential before checking tracks
143
144 // Keep all different tracks or only the best ones as required
3304fa09 145 if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
ea94c18b 146 else RemoveDoubleTracks();
147
148 AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
149
8cc77df0 150}
151
ea94c18b 152 //__________________________________________________________________________
153void AliMUONTrackReconstructorK::RetraceTrack(AliMUONTrack &trackCandidate, Bool_t resetSeed)
8cc77df0 154{
96ebe67e 155 /// Re-run the kalman filter from the most downstream cluster to the most uptream one
ea94c18b 156 AliDebug(1,"Enter RetraceTrack");
157
96ebe67e 158 AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Last();
ea94c18b 159
96ebe67e 160 // Reset the "seed" (= track parameters and their covariances at last cluster) if required
ea94c18b 161 if (resetSeed) {
96ebe67e 162 // => Shift track parameters at the position of the last cluster
163 AliMUONVCluster* cluster = startingTrackParam->GetClusterPtr();
164 startingTrackParam->SetNonBendingCoor(cluster->GetX());
165 startingTrackParam->SetBendingCoor(cluster->GetY());
ea94c18b 166
96ebe67e 167 // => Re-compute and reset track parameters covariances at last cluster (as if the other clusters did not exist)
ea94c18b 168 const TMatrixD& kParamCov = startingTrackParam->GetCovariances();
169 TMatrixD newParamCov(5,5);
170 newParamCov.Zero();
171 // Non bending plane
96ebe67e 172 newParamCov(0,0) = cluster->GetErrX2();
ea94c18b 173 newParamCov(1,1) = 100.*kParamCov(1,1);
174 // Bending plane
96ebe67e 175 newParamCov(2,2) = cluster->GetErrY2();
ea94c18b 176 newParamCov(3,3) = 100.*kParamCov(3,3);
177 // Inverse bending momentum
178 newParamCov(4,4) = 0.5*startingTrackParam->GetInverseBendingMomentum() * 0.5*startingTrackParam->GetInverseBendingMomentum();
179 startingTrackParam->SetCovariances(newParamCov);
180
181 // Reset the track chi2
182 startingTrackParam->SetTrackChi2(0.);
183
184 }
185
186 // Redo the tracking
187 RetracePartialTrack(trackCandidate, startingTrackParam);
188
189}
8cc77df0 190
ea94c18b 191 //__________________________________________________________________________
192void AliMUONTrackReconstructorK::RetracePartialTrack(AliMUONTrack &trackCandidate, const AliMUONTrackParam* startingTrackParam)
193{
96ebe67e 194 /// Re-run the kalman filter from the cluster attached to startingTrackParam to the most uptream cluster
ea94c18b 195 AliDebug(1,"Enter RetracePartialTrack");
196
197 // Printout for debuging
198 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 199 cout << "RetracePartialTrack: track chi2 before re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
ea94c18b 200 }
201
202 // Reset the track chi2
96ebe67e 203 trackCandidate.SetGlobalChi2(startingTrackParam->GetTrackChi2());
ea94c18b 204
96ebe67e 205 // loop over attached clusters until the first one and recompute track parameters and covariances using kalman filter
206 Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() - 1;
ea94c18b 207 Int_t currentChamber;
96ebe67e 208 Double_t addChi2TrackAtCluster;
209 AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam);
210 while (trackParamAtCluster) {
ea94c18b 211
212 // reset track parameters and their covariances
96ebe67e 213 trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
214 trackParamAtCluster->SetZ(startingTrackParam->GetZ());
215 trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
ea94c18b 216
217 // add MCS effect
96ebe67e 218 AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
ea94c18b 219
220 // reset propagator for smoother
96ebe67e 221 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) trackParamAtCluster->ResetPropagator();
ea94c18b 222
223 // add MCS in missing chambers if any (at most 2 chambers can be missing according to tracking criteria)
96ebe67e 224 currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
ea94c18b 225 while (currentChamber < expectedChamber) {
226 // extrapolation to the missing chamber (update the propagator)
96ebe67e 227 AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber),
3304fa09 228 AliMUONReconstructor::GetRecoParam()->UseSmoother());
ea94c18b 229 // add MCS effect
96ebe67e 230 AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
ea94c18b 231 expectedChamber--;
232 }
233
96ebe67e 234 // extrapolation to the plane of the cluster attached to the current trackParamAtCluster (update the propagator)
235 AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ(),
3304fa09 236 AliMUONReconstructor::GetRecoParam()->UseSmoother());
ea94c18b 237
3304fa09 238 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
ea94c18b 239 // save extrapolated parameters for smoother
96ebe67e 240 trackParamAtCluster->SetExtrapParameters(trackParamAtCluster->GetParameters());
ea94c18b 241
242 // save extrapolated covariance matrix for smoother
96ebe67e 243 trackParamAtCluster->SetExtrapCovariances(trackParamAtCluster->GetCovariances());
ea94c18b 244 }
245
96ebe67e 246 // Compute new track parameters including "clusterCh2" using kalman filter
247 addChi2TrackAtCluster = RunKalmanFilter(*trackParamAtCluster);
ea94c18b 248
249 // Update the track chi2
96ebe67e 250 trackCandidate.SetGlobalChi2(trackCandidate.GetGlobalChi2() + addChi2TrackAtCluster);
251 trackParamAtCluster->SetTrackChi2(trackCandidate.GetGlobalChi2());
ea94c18b 252
96ebe67e 253 // prepare next step
254 expectedChamber = currentChamber - 1;
255 startingTrackParam = trackParamAtCluster;
256 trackParamAtCluster = (AliMUONTrackParam*) (trackCandidate.GetTrackParamAtCluster()->Before(startingTrackParam));
ea94c18b 257 }
258
259 // Printout for debuging
260 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 261 cout << "RetracePartialTrack: track chi2 after re-tracking: " << trackCandidate.GetGlobalChi2() << endl;
ea94c18b 262 }
263
264}
8cc77df0 265
ea94c18b 266 //__________________________________________________________________________
96ebe67e 267void AliMUONTrackReconstructorK::FollowTracks(const AliMUONVClusterStore& clusterStore)
ea94c18b 268{
269 /// Follow tracks in stations(1..) 3, 2 and 1
270 AliDebug(1,"Enter FollowTracks");
271
272 AliMUONTrack *track;
273 Int_t currentNRecTracks;
96ebe67e 274 Bool_t clusterFound;
ea94c18b 275
276 for (Int_t station = 2; station >= 0; station--) {
277
278 // Save the actual number of reconstructed track in case of
279 // tracks are added or suppressed during the tracking procedure
280 // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
281 currentNRecTracks = fNRecTracks;
282
283 for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
284 AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
285
286 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
287
288 // Printout for debuging
289 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
96ebe67e 290 cout<<endl<<"Track parameter covariances at first cluster:"<<endl;
291 ((AliMUONTrackParam*) (track->GetTrackParamAtCluster()->First()))->GetCovariances().Print();
ea94c18b 292 }
293
96ebe67e 294 // Look for compatible cluster(s) in station(0..) "station"
295 clusterFound = FollowTrackInStation(*track, clusterStore, station);
ea94c18b 296
297 // Try to recover track if required
96ebe67e 298 if (!clusterFound && AliMUONReconstructor::GetRecoParam()->RecoverTracks())
299 clusterFound = RecoverTrack(*track, clusterStore, station);
ea94c18b 300
96ebe67e 301 // remove track if no cluster found
302 if (!clusterFound) {
ea94c18b 303 fRecTracksPtr->Remove(track);
304 fNRecTracks--;
8cc77df0 305 }
ea94c18b 306
8cc77df0 307 }
ea94c18b 308
309 fRecTracksPtr->Compress(); // this is essential before checking tracks
310
ea94c18b 311 // Keep only the best tracks if required
3304fa09 312 if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
ea94c18b 313
314 }
315
316}
b0e4c598 317
ea94c18b 318 //__________________________________________________________________________
96ebe67e 319Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, Int_t nextStation)
ea94c18b 320{
96ebe67e 321 /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
ea94c18b 322 /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
323 /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
324 /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
325 /// Remove the obsolete "trackCandidate" at the end.
96ebe67e 326 /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
327 /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
ea94c18b 328 AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
329
330 // Order the chamber according to the propagation direction (tracking starts with chamber 2):
331 // - nextStation == station(1...) 5 => forward propagation
332 // - nextStation < station(1...) 5 => backward propagation
333 Int_t ch1, ch2;
334 if (nextStation==4) {
335 ch1 = 2*nextStation+1;
336 ch2 = 2*nextStation;
337 } else {
338 ch1 = 2*nextStation;
339 ch2 = 2*nextStation+1;
340 }
341
96ebe67e 342 Double_t chi2OfCluster;
343 Double_t maxChi2OfCluster = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
344 AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
345 Double_t addChi2TrackAtCluster1;
346 Double_t addChi2TrackAtCluster2;
347 Double_t bestAddChi2TrackAtCluster1 = 1.e10;
348 Double_t bestAddChi2TrackAtCluster2 = 1.e10;
349 Bool_t foundOneCluster = kFALSE;
350 Bool_t foundTwoClusters = kFALSE;
ea94c18b 351 AliMUONTrack *newTrack = 0x0;
96ebe67e 352 AliMUONVCluster *clusterCh1, *clusterCh2;
ea94c18b 353 AliMUONTrackParam extrapTrackParam;
96ebe67e 354 AliMUONTrackParam extrapTrackParamAtCluster1;
355 AliMUONTrackParam extrapTrackParamAtCluster2;
356 AliMUONTrackParam bestTrackParamAtCluster1;
357 AliMUONTrackParam bestTrackParamAtCluster2;
358
359 Int_t nClusters = clusterStore.GetSize();
360 Bool_t *clusterCh1Used = new Bool_t[nClusters];
361 for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
362 Int_t iCluster1;
363
019df241 364 // Get track parameters
96ebe67e 365 AliMUONTrackParam extrapTrackParamAtCh(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First());
ea94c18b 366
367 // Add MCS effect
019df241 368 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
ea94c18b 369
370 // reset propagator for smoother
3304fa09 371 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) extrapTrackParamAtCh.ResetPropagator();
ea94c18b 372
373 // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria)
96ebe67e 374 if (ch1 < ch2 && extrapTrackParamAtCh.GetClusterPtr()->GetChamberId() > ch2 + 1) {
ea94c18b 375 // extrapolation to the missing chamber
3304fa09 376 AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2 + 1),
377 AliMUONReconstructor::GetRecoParam()->UseSmoother());
ea94c18b 378 // add MCS effect
019df241 379 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
ea94c18b 380 }
381
382 //Extrapolate trackCandidate to chamber "ch2"
3304fa09 383 AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2),
384 AliMUONReconstructor::GetRecoParam()->UseSmoother());
ea94c18b 385
386 // Printout for debuging
387 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
96ebe67e 388 cout<<endl<<"Track parameter covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
019df241 389 extrapTrackParamAtCh.GetCovariances().Print();
ea94c18b 390 }
391
392 // Printout for debuging
393 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 394 cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
ea94c18b 395 }
396
96ebe67e 397 // Create iterators to loop over clusters in both chambers
398 TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
399 TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
400
401 // look for candidates in chamber 2
402 while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
ea94c18b 403
96ebe67e 404 // try to add the current cluster fast
405 if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
019df241 406
96ebe67e 407 // try to add the current cluster accuratly
408 chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2,
409 AliMUONReconstructor::GetRecoParam()->UseSmoother());
ea94c18b 410
96ebe67e 411 // if good chi2 then try to attach a cluster in the other chamber too
412 if (chi2OfCluster < maxChi2OfCluster) {
413 Bool_t foundSecondCluster = kFALSE;
ea94c18b 414
415 // Printout for debuging
416 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 417 cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
418 << " (Chi2 = " << chi2OfCluster << ")" << endl;
419 cout << " look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
8cc77df0 420 }
ea94c18b 421
3304fa09 422 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
ea94c18b 423 // save extrapolated parameters for smoother
96ebe67e 424 extrapTrackParamAtCluster2.SetExtrapParameters(extrapTrackParamAtCluster2.GetParameters());
ea94c18b 425
426 // save extrapolated covariance matrix for smoother
96ebe67e 427 extrapTrackParamAtCluster2.SetExtrapCovariances(extrapTrackParamAtCluster2.GetCovariances());
ea94c18b 428 }
429
96ebe67e 430 // Compute new track parameters including "clusterCh2" using kalman filter
431 addChi2TrackAtCluster2 = RunKalmanFilter(extrapTrackParamAtCluster2);
ea94c18b 432
433 // copy new track parameters for next step
96ebe67e 434 extrapTrackParam = extrapTrackParamAtCluster2;
ea94c18b 435
436 // add MCS effect
437 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
438
439 // reset propagator for smoother
3304fa09 440 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) extrapTrackParam.ResetPropagator();
ea94c18b 441
019df241 442 //Extrapolate track parameters to chamber "ch1"
3304fa09 443 AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1),
444 AliMUONReconstructor::GetRecoParam()->UseSmoother());
019df241 445
96ebe67e 446 // reset cluster iterator of chamber 1
447 nextInCh1.Reset();
448 iCluster1 = -1;
449
450 // look for second candidates in chamber 1
451 while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
452 iCluster1++;
ea94c18b 453
96ebe67e 454 // try to add the current cluster fast
455 if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
019df241 456
96ebe67e 457 // try to add the current cluster accuratly
458 chi2OfCluster = TryOneCluster(extrapTrackParam, clusterCh1, extrapTrackParamAtCluster1,
3304fa09 459 AliMUONReconstructor::GetRecoParam()->UseSmoother());
019df241 460
96ebe67e 461 // if good chi2 then consider to add the 2 clusters to the "trackCandidate"
462 if (chi2OfCluster < maxChi2OfCluster) {
463 foundSecondCluster = kTRUE;
464 foundTwoClusters = kTRUE;
ea94c18b 465
466 // Printout for debuging
467 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 468 cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
469 << " (Chi2 = " << chi2OfCluster << ")" << endl;
ea94c18b 470 }
019df241 471
3304fa09 472 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
ea94c18b 473 // save extrapolated parameters for smoother
96ebe67e 474 extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
ea94c18b 475
476 // save extrapolated covariance matrix for smoother
96ebe67e 477 extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
ea94c18b 478 }
479
96ebe67e 480 // Compute new track parameters including "clusterCh1" using kalman filter
481 addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
ea94c18b 482
3304fa09 483 if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
96ebe67e 484 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
ea94c18b 485 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
96ebe67e 486 UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2,addChi2TrackAtCluster1,addChi2TrackAtCluster2);
ea94c18b 487 fNRecTracks++;
488
489 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
490 // (going in the right direction)
491 if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
492
96ebe67e 493 // Tag clusterCh1 as used
494 clusterCh1Used[iCluster1] = kTRUE;
ea94c18b 495
496 // Printout for debuging
497 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 498 cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
ea94c18b 499 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
500 }
501
96ebe67e 502 } else if (addChi2TrackAtCluster1+addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1+bestAddChi2TrackAtCluster2) {
503 // keep track of the best couple of clusters
504 bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
505 bestAddChi2TrackAtCluster2 = addChi2TrackAtCluster2;
506 bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
507 bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
ea94c18b 508 }
509
b0e4c598 510 }
ea94c18b 511
8cc77df0 512 }
ea94c18b 513
96ebe67e 514 // if no clusterCh1 found then consider to add clusterCh2 only
515 if (!foundSecondCluster) {
516 foundOneCluster = kTRUE;
ea94c18b 517
3304fa09 518 if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
96ebe67e 519 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
ea94c18b 520 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
96ebe67e 521 UpdateTrack(*newTrack,extrapTrackParamAtCluster2,addChi2TrackAtCluster2);
ea94c18b 522 fNRecTracks++;
523
524 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
525 // (going in the right direction)
526 if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
527
528 // Printout for debuging
529 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 530 cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
ea94c18b 531 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
532 }
533
96ebe67e 534 } else if (!foundTwoClusters && addChi2TrackAtCluster2 < bestAddChi2TrackAtCluster1) {
535 // keep track of the best single cluster except if a couple of clusters has already been found
536 bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster2;
537 bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
ea94c18b 538 }
539
540 }
541
8cc77df0 542 }
ea94c18b 543
544 }
545
546 // look for candidates in chamber 1 not already attached to a track
96ebe67e 547 // if we want to keep all possible tracks or if no good couple of clusters has been found
548 if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
ea94c18b 549
550 // Printout for debuging
551 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 552 cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
8cc77df0 553 }
ea94c18b 554
555 // add MCS effect for next step
019df241 556 AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
557
558 //Extrapolate trackCandidate to chamber "ch1"
3304fa09 559 AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1),
560 AliMUONReconstructor::GetRecoParam()->UseSmoother());
019df241 561
96ebe67e 562 // reset cluster iterator of chamber 1
563 nextInCh1.Reset();
564 iCluster1 = -1;
565
566 // look for second candidates in chamber 1
567 while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
568 iCluster1++;
ea94c18b 569
96ebe67e 570 if (clusterCh1Used[iCluster1]) continue; // Skip clusters already used
ea94c18b 571
96ebe67e 572 // try to add the current cluster fast
573 if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
019df241 574
96ebe67e 575 // try to add the current cluster accuratly
576 chi2OfCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1,
577 AliMUONReconstructor::GetRecoParam()->UseSmoother());
019df241 578
96ebe67e 579 // if good chi2 then consider to add clusterCh1
580 // We do not try to attach a cluster in the other chamber too since it has already been done above
581 if (chi2OfCluster < maxChi2OfCluster) {
582 foundOneCluster = kTRUE;
ea94c18b 583
584 // Printout for debuging
585 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 586 cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
587 << " (Chi2 = " << chi2OfCluster << ")" << endl;
ea94c18b 588 }
019df241 589
3304fa09 590 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) {
ea94c18b 591 // save extrapolated parameters for smoother
96ebe67e 592 extrapTrackParamAtCluster1.SetExtrapParameters(extrapTrackParamAtCluster1.GetParameters());
ea94c18b 593
594 // save extrapolated covariance matrix for smoother
96ebe67e 595 extrapTrackParamAtCluster1.SetExtrapCovariances(extrapTrackParamAtCluster1.GetCovariances());
ea94c18b 596 }
597
96ebe67e 598 // Compute new track parameters including "clusterCh1" using kalman filter
599 addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
ea94c18b 600
3304fa09 601 if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
96ebe67e 602 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
ea94c18b 603 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
96ebe67e 604 UpdateTrack(*newTrack,extrapTrackParamAtCluster1,addChi2TrackAtCluster1);
ea94c18b 605 fNRecTracks++;
606
607 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
608 // (going in the right direction)
609 if (nextStation == 4) RetraceTrack(*newTrack,kTRUE);
610
611 // Printout for debuging
612 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 613 cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
ea94c18b 614 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
615 }
616
96ebe67e 617 } else if (addChi2TrackAtCluster1 < bestAddChi2TrackAtCluster1) {
618 // keep track of the best single cluster except if a couple of clusters has already been found
619 bestAddChi2TrackAtCluster1 = addChi2TrackAtCluster1;
620 bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
ea94c18b 621 }
622
623 }
624
8cc77df0 625 }
ea94c18b 626
627 }
628
629 // fill out the best track if required else clean up the fRecTracksPtr array
3304fa09 630 if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) {
96ebe67e 631 if (foundTwoClusters) {
632 UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2,bestAddChi2TrackAtCluster1,bestAddChi2TrackAtCluster2);
ea94c18b 633
634 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
635 // (going in the right direction)
636 if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE);
8cc77df0 637
ea94c18b 638 // Printout for debuging
639 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 640 cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
ea94c18b 641 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
642 }
643
96ebe67e 644 } else if (foundOneCluster) {
645 UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestAddChi2TrackAtCluster1);
ea94c18b 646
647 // if we are arrived on station(1..) 5, recompute track parameters and covariances starting from this station
648 // (going in the right direction)
649 if (nextStation == 4) RetraceTrack(trackCandidate,kTRUE);
650
651 // Printout for debuging
652 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
96ebe67e 653 cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
ea94c18b 654 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
655 }
656
17773e2e 657 } else {
96ebe67e 658 delete [] clusterCh1Used;
17773e2e 659 return kFALSE;
660 }
96ebe67e 661 } else if (foundOneCluster || foundTwoClusters) {
ea94c18b 662
663 // remove obsolete track
664 fRecTracksPtr->Remove(&trackCandidate);
665 fNRecTracks--;
666
17773e2e 667 } else {
96ebe67e 668 delete [] clusterCh1Used;
17773e2e 669 return kFALSE;
670 }
96ebe67e 671 delete [] clusterCh1Used;
ea94c18b 672 return kTRUE;
673
674}
8cc77df0 675
ea94c18b 676 //__________________________________________________________________________
96ebe67e 677Double_t AliMUONTrackReconstructorK::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster)
ea94c18b 678{
96ebe67e 679 /// Compute new track parameters and their covariances including new cluster using kalman filter
ea94c18b 680 /// return the additional track chi2
681 AliDebug(1,"Enter RunKalmanFilter");
682
683 // Get actual track parameters (p)
96ebe67e 684 TMatrixD param(trackParamAtCluster.GetParameters());
ea94c18b 685
96ebe67e 686 // Get new cluster parameters (m)
687 AliMUONVCluster *cluster = trackParamAtCluster.GetClusterPtr();
688 TMatrixD clusterParam(5,1);
689 clusterParam.Zero();
690 clusterParam(0,0) = cluster->GetX();
691 clusterParam(2,0) = cluster->GetY();
ea94c18b 692
693 // Compute the actual parameter weight (W)
96ebe67e 694 TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
ea94c18b 695 if (paramWeight.Determinant() != 0) {
696 paramWeight.Invert();
697 } else {
698 AliWarning(" Determinant = 0");
699 return 1.e10;
8cc77df0 700 }
ea94c18b 701
96ebe67e 702 // Compute the new cluster weight (U)
703 TMatrixD clusterWeight(5,5);
704 clusterWeight.Zero();
705 clusterWeight(0,0) = 1. / cluster->GetErrX2();
706 clusterWeight(2,2) = 1. / cluster->GetErrY2();
8cc77df0 707
ea94c18b 708 // Compute the new parameters covariance matrix ( (W+U)^-1 )
96ebe67e 709 TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
ea94c18b 710 if (newParamCov.Determinant() != 0) {
711 newParamCov.Invert();
712 } else {
713 AliWarning(" Determinant = 0");
714 return 1.e10;
715 }
716
717 // Save the new parameters covariance matrix
96ebe67e 718 trackParamAtCluster.SetCovariances(newParamCov);
ea94c18b 719
720 // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
96ebe67e 721 TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
722 TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
ea94c18b 723 TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
724 newParam += param; // ((W+U)^-1)U(m-p) + p
725
726 // Save the new parameters
96ebe67e 727 trackParamAtCluster.SetParameters(newParam);
ea94c18b 728
729 // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
730 tmp = newParam; // p'
731 tmp -= param; // (p'-p)
732 TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
733 TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
734 tmp = newParam; // p'
96ebe67e 735 tmp -= clusterParam; // (p'-m)
736 TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
ea94c18b 737 addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
738
739 return addChi2Track(0,0);
740
8cc77df0 741}
742
ea94c18b 743 //__________________________________________________________________________
96ebe67e 744void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster, Double_t addChi2)
8cc77df0 745{
96ebe67e 746 /// Add 1 cluster to the track candidate
ea94c18b 747 /// Update chi2 of the track
748
96ebe67e 749 // Flag cluster as being not removable
750 trackParamAtCluster.SetRemovable(kFALSE);
751 trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
5f4ceff2 752
96ebe67e 753 // Update the track chi2 into trackParamAtCluster
754 trackParamAtCluster.SetTrackChi2(track.GetGlobalChi2() + addChi2);
ea94c18b 755
756 // Update the chi2 of the new track
96ebe67e 757 track.SetGlobalChi2(trackParamAtCluster.GetTrackChi2());
ea94c18b 758
96ebe67e 759 // Update array of TrackParamAtCluster
760 track.AddTrackParamAtCluster(trackParamAtCluster,*(trackParamAtCluster.GetClusterPtr()));
761 track.GetTrackParamAtCluster()->Sort();
ea94c18b 762
763}
8cc77df0 764
ea94c18b 765 //__________________________________________________________________________
96ebe67e 766void AliMUONTrackReconstructorK::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2,
767 Double_t addChi2AtCluster1, Double_t addChi2AtCluster2)
ea94c18b 768{
96ebe67e 769 /// Add 2 clusters to the track candidate (order is important)
ea94c18b 770 /// Update track and local chi2
771
96ebe67e 772 // Update local chi2 at first cluster
773 AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
774 Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
775 Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
776 Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
777 deltaY*deltaY / cluster1->GetErrY2();
778 trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
ea94c18b 779
96ebe67e 780 // Flag first cluster as being removable
781 trackParamAtCluster1.SetRemovable(kTRUE);
ea94c18b 782
96ebe67e 783 // Update local chi2 at second cluster
784 AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
785 AliMUONTrackParam extrapTrackParamAtCluster2(trackParamAtCluster1);
786 AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParamAtCluster2, trackParamAtCluster2.GetZ());
787 deltaX = extrapTrackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
788 deltaY = extrapTrackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
789 Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
790 deltaY*deltaY / cluster2->GetErrY2();
791 trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
ea94c18b 792
96ebe67e 793 // Flag second cluster as being removable
794 trackParamAtCluster2.SetRemovable(kTRUE);
ea94c18b 795
96ebe67e 796 // Update the track chi2 into trackParamAtCluster2
797 trackParamAtCluster2.SetTrackChi2(track.GetGlobalChi2() + addChi2AtCluster2);
ea94c18b 798
96ebe67e 799 // Update the track chi2 into trackParamAtCluster1
800 trackParamAtCluster1.SetTrackChi2(trackParamAtCluster2.GetTrackChi2() + addChi2AtCluster1);
ea94c18b 801
802 // Update the chi2 of the new track
96ebe67e 803 track.SetGlobalChi2(trackParamAtCluster1.GetTrackChi2());
ea94c18b 804
96ebe67e 805 // Update array of trackParamAtCluster
806 track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
807 track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
808 track.GetTrackParamAtCluster()->Sort();
ea94c18b 809
810}
8cc77df0 811
ea94c18b 812 //__________________________________________________________________________
96ebe67e 813Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, Int_t nextStation)
ea94c18b 814{
815 /// Try to recover the track candidate in the next station
96ebe67e 816 /// by removing the worst of the two clusters attached in the current station
ea94c18b 817 /// Return kTRUE if recovering succeeds
818 AliDebug(1,"Enter RecoverTrack");
819
96ebe67e 820 // Do not try to recover track until we have attached cluster(s) on station(1..) 3
ea94c18b 821 if (nextStation > 1) return kFALSE;
822
96ebe67e 823 Int_t worstClusterNumber = -1;
824 Double_t localChi2, worstLocalChi2 = 0.;
ea94c18b 825
96ebe67e 826 // Look for the cluster to remove
827 for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
828 AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
ea94c18b 829
96ebe67e 830 // check if current cluster is removable
831 if (!trackParamAtCluster->IsRemovable()) return kFALSE;
ea94c18b 832
96ebe67e 833 // Pick up cluster with the worst chi2
834 localChi2 = trackParamAtCluster->GetLocalChi2();
835 if (localChi2 > worstLocalChi2) {
836 worstLocalChi2 = localChi2;
837 worstClusterNumber = clusterNumber;
8cc77df0 838 }
ea94c18b 839 }
840
96ebe67e 841 // Reset best cluster as being NOT removable
842 ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt((worstClusterNumber+1)%2))->SetRemovable(kFALSE);
ea94c18b 843
96ebe67e 844 // Remove the worst cluster
845 trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
ea94c18b 846
96ebe67e 847 // Re-calculate track parameters at the (new) first cluster
848 RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(1));
ea94c18b 849
96ebe67e 850 // Look for new cluster(s) in next station
851 return FollowTrackInStation(trackCandidate, clusterStore, nextStation);
ea94c18b 852
8cc77df0 853}
854
855 //__________________________________________________________________________
ea94c18b 856Bool_t AliMUONTrackReconstructorK::RunSmoother(AliMUONTrack &track)
8cc77df0 857{
ea94c18b 858 /// Compute new track parameters and their covariances using smoother
3304fa09 859 AliDebug(1,"Enter UseSmoother");
ea94c18b 860
96ebe67e 861 AliMUONTrackParam *previousTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
ea94c18b 862
96ebe67e 863 // Smoothed parameters and covariances at first cluster = filtered parameters and covariances
ea94c18b 864 previousTrackParam->SetSmoothParameters(previousTrackParam->GetParameters());
865 previousTrackParam->SetSmoothCovariances(previousTrackParam->GetCovariances());
866
96ebe67e 867 // Compute local chi2 at first cluster
868 AliMUONVCluster *cluster = previousTrackParam->GetClusterPtr();
869 Double_t dX = cluster->GetX() - previousTrackParam->GetNonBendingCoor();
870 Double_t dY = cluster->GetY() - previousTrackParam->GetBendingCoor();
871 Double_t localChi2 = dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
ea94c18b 872
96ebe67e 873 // Save local chi2 at first cluster
ea94c18b 874 previousTrackParam->SetLocalChi2(localChi2);
875
96ebe67e 876 AliMUONTrackParam *currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
ea94c18b 877 while (currentTrackParam) {
878
879 // Get variables
880 const TMatrixD &extrapParameters = previousTrackParam->GetExtrapParameters(); // X(k+1 k)
881 const TMatrixD &filteredParameters = currentTrackParam->GetParameters(); // X(k k)
882 const TMatrixD &previousSmoothParameters = previousTrackParam->GetSmoothParameters(); // X(k+1 n)
883 const TMatrixD &propagator = previousTrackParam->GetPropagator(); // F(k)
884 const TMatrixD &extrapCovariances = previousTrackParam->GetExtrapCovariances(); // C(k+1 k)
885 const TMatrixD &filteredCovariances = currentTrackParam->GetCovariances(); // C(k k)
886 const TMatrixD &previousSmoothCovariances = previousTrackParam->GetSmoothCovariances(); // C(k+1 n)
887
888 // Compute smoother gain: A(k) = C(kk) * F(f)^t * (C(k+1 k))^-1
889 TMatrixD extrapWeight(extrapCovariances);
890 if (extrapWeight.Determinant() != 0) {
891 extrapWeight.Invert(); // (C(k+1 k))^-1
892 } else {
893 AliWarning(" Determinant = 0");
894 return kFALSE;
895 }
896 TMatrixD smootherGain(filteredCovariances,TMatrixD::kMultTranspose,propagator); // C(kk) * F(f)^t
897 smootherGain *= extrapWeight; // C(kk) * F(f)^t * (C(k+1 k))^-1
898
899 // Compute smoothed parameters: X(k n) = X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
900 TMatrixD tmpParam(previousSmoothParameters,TMatrixD::kMinus,extrapParameters); // X(k+1 n) - X(k+1 k)
901 TMatrixD smoothParameters(smootherGain,TMatrixD::kMult,tmpParam); // A(k) * (X(k+1 n) - X(k+1 k))
902 smoothParameters += filteredParameters; // X(k k) + A(k) * (X(k+1 n) - X(k+1 k))
903
904 // Save smoothed parameters
905 currentTrackParam->SetSmoothParameters(smoothParameters);
906
907 // Compute smoothed covariances: C(k n) = C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
908 TMatrixD tmpCov(previousSmoothCovariances,TMatrixD::kMinus,extrapCovariances); // C(k+1 n) - C(k+1 k)
909 TMatrixD tmpCov2(tmpCov,TMatrixD::kMultTranspose,smootherGain); // (C(k+1 n) - C(k+1 k)) * (A(k))^t
910 TMatrixD smoothCovariances(smootherGain,TMatrixD::kMult,tmpCov2); // A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
911 smoothCovariances += filteredCovariances; // C(k k) + A(k) * (C(k+1 n) - C(k+1 k)) * (A(k))^t
912
913 // Save smoothed covariances
914 currentTrackParam->SetSmoothCovariances(smoothCovariances);
915
96ebe67e 916 // Compute smoothed residual: r(k n) = cluster - X(k n)
917 cluster = currentTrackParam->GetClusterPtr();
ea94c18b 918 TMatrixD smoothResidual(2,1);
919 smoothResidual.Zero();
96ebe67e 920 smoothResidual(0,0) = cluster->GetX() - smoothParameters(0,0);
921 smoothResidual(1,0) = cluster->GetY() - smoothParameters(2,0);
ea94c18b 922
96ebe67e 923 // Compute weight of smoothed residual: W(k n) = (clusterCov - C(k n))^-1
ea94c18b 924 TMatrixD smoothResidualWeight(2,2);
96ebe67e 925 smoothResidualWeight(0,0) = cluster->GetErrX2() - smoothCovariances(0,0);
ea94c18b 926 smoothResidualWeight(0,1) = - smoothCovariances(0,2);
927 smoothResidualWeight(1,0) = - smoothCovariances(2,0);
96ebe67e 928 smoothResidualWeight(1,1) = cluster->GetErrY2() - smoothCovariances(2,2);
ea94c18b 929 if (smoothResidualWeight.Determinant() != 0) {
930 smoothResidualWeight.Invert();
931 } else {
932 AliWarning(" Determinant = 0");
933 return kFALSE;
934 }
935
936 // Compute local chi2 = (r(k n))^t * W(k n) * r(k n)
937 TMatrixD tmpChi2(smoothResidual,TMatrixD::kTransposeMult,smoothResidualWeight); // (r(k n))^t * W(k n)
938 TMatrixD localChi2(tmpChi2,TMatrixD::kMult,smoothResidual); // (r(k n))^t * W(k n) * r(k n)
939
940 // Save local chi2
941 currentTrackParam->SetLocalChi2(localChi2(0,0));
942
943 previousTrackParam = currentTrackParam;
96ebe67e 944 currentTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(previousTrackParam);
ea94c18b 945 }
946
947 return kTRUE;
948
8cc77df0 949}
950
951 //__________________________________________________________________________
96ebe67e 952void AliMUONTrackReconstructorK::ComplementTracks(const AliMUONVClusterStore& clusterStore)
b709ac13 953{
954 /// Complete tracks by adding missing clusters (if there is an overlap between
955 /// two detection elements, the track may have two clusters in the same chamber)
956 /// Recompute track parameters and covariances at each clusters
957 AliDebug(1,"Enter ComplementTracks");
958
959 Int_t chamberId, detElemId;
96ebe67e 960 Double_t chi2OfCluster, addChi2TrackAtCluster, bestAddChi2TrackAtCluster;
961 Double_t maxChi2OfCluster = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() *
962 AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
963 Bool_t foundOneCluster, trackModified;
964 AliMUONVCluster *cluster;
965 AliMUONTrackParam *trackParam, *previousTrackParam, *nextTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
b709ac13 966
967 // Remove double track to complete only "good" tracks
968 RemoveDoubleTracks();
969
970 AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
971 while (track) {
972 trackModified = kFALSE;
973
96ebe67e 974 trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
b709ac13 975 previousTrackParam = trackParam;
976 while (trackParam) {
96ebe67e 977 foundOneCluster = kFALSE;
978 bestAddChi2TrackAtCluster = 1.e10;
979 chamberId = trackParam->GetClusterPtr()->GetChamberId();
980 detElemId = trackParam->GetClusterPtr()->GetDetElemId();
b709ac13 981
982 // prepare nextTrackParam before adding new cluster because of the sorting
96ebe67e 983 nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
b709ac13 984
96ebe67e 985 // Create iterators to loop over clusters in current chamber
986 TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
b709ac13 987
988 // look for one second candidate in the same chamber
96ebe67e 989 while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
b709ac13 990
991 // look for a cluster in another detection element
96ebe67e 992 if (cluster->GetDetElemId() == detElemId) continue;
b709ac13 993
96ebe67e 994 // try to add the current cluster fast
995 if (!TryOneClusterFast(*trackParam, cluster)) continue;
b709ac13 996
96ebe67e 997 // try to add the current cluster accurately
b709ac13 998 // never use track parameters at last cluster because the covariance matrix is meaningless
96ebe67e 999 if (nextTrackParam) chi2OfCluster = TryOneCluster(*trackParam, cluster, trackParamAtCluster);
1000 else chi2OfCluster = TryOneCluster(*previousTrackParam, cluster, trackParamAtCluster);
b709ac13 1001
1002 // if good chi2 then consider to add this cluster to the track
96ebe67e 1003 if (chi2OfCluster < maxChi2OfCluster) {
b709ac13 1004
96ebe67e 1005 // Compute local track parameters including current cluster using kalman filter
1006 addChi2TrackAtCluster = RunKalmanFilter(trackParamAtCluster);
b709ac13 1007
1008 // keep track of the best cluster
96ebe67e 1009 if (addChi2TrackAtCluster < bestAddChi2TrackAtCluster) {
1010 bestAddChi2TrackAtCluster = addChi2TrackAtCluster;
1011 bestTrackParamAtCluster = trackParamAtCluster;
1012 foundOneCluster = kTRUE;
b709ac13 1013 }
1014
1015 }
1016
1017 }
1018
1019 // add new cluster if any
96ebe67e 1020 if (foundOneCluster) {
1021 UpdateTrack(*track,bestTrackParamAtCluster,bestAddChi2TrackAtCluster);
1022 bestTrackParamAtCluster.SetAloneInChamber(kFALSE);
b709ac13 1023 trackParam->SetAloneInChamber(kFALSE);
1024 trackModified = kTRUE;
1025 }
1026
1027 previousTrackParam = trackParam;
1028 trackParam = nextTrackParam;
1029 }
1030
1031 // re-compute track parameters using kalman filter if needed
1032 if (trackModified) RetraceTrack(*track,kTRUE);
1033
1034 track = (AliMUONTrack*) fRecTracksPtr->After(track);
1035 }
1036
1037}
1038
1039//__________________________________________________________________________
ea94c18b 1040void AliMUONTrackReconstructorK::ImproveTracks()
8cc77df0 1041{
ea94c18b 1042 /// Improve tracks by removing clusters with local chi2 highter than the defined cut
1043 /// Recompute track parameters and covariances at the remaining clusters
1044 AliDebug(1,"Enter ImproveTracks");
1045
1046 Double_t localChi2, worstLocalChi2;
b709ac13 1047 Int_t worstChamber, previousChamber;
1048 AliMUONTrack *track, *nextTrack;
96ebe67e 1049 AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster, *previousTrackParam, *nextTrackParam;
ea94c18b 1050 Bool_t smoothed;
3304fa09 1051 Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement() *
1052 AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement();
ea94c18b 1053
1054 // Remove double track to improve only "good" tracks
1055 RemoveDoubleTracks();
1056
b709ac13 1057 track = (AliMUONTrack*) fRecTracksPtr->First();
8cc77df0 1058 while (track) {
ea94c18b 1059
b709ac13 1060 // prepare next track in case the actual track is suppressed
1061 nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track);
1062
ea94c18b 1063 while (!track->IsImproved()) {
1064
1065 // Run smoother if required
1066 smoothed = kFALSE;
3304fa09 1067 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) smoothed = RunSmoother(*track);
ea94c18b 1068
1069 // Use standard procedure to compute local chi2 if smoother not required or not working
1070 if (!smoothed) {
1071
1072 // Update track parameters and covariances
96ebe67e 1073 track->UpdateCovTrackParamAtCluster();
ea94c18b 1074
96ebe67e 1075 // Compute local chi2 of each clusters
ea94c18b 1076 track->ComputeLocalChi2(kTRUE);
1077 }
1078
96ebe67e 1079 // Look for the cluster to remove
1080 worstTrackParamAtCluster = NULL;
ea94c18b 1081 worstLocalChi2 = 0.;
96ebe67e 1082 trackParamAtCluster = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1083 while (trackParamAtCluster) {
ea94c18b 1084
96ebe67e 1085 // Pick up cluster with the worst chi2
1086 localChi2 = trackParamAtCluster->GetLocalChi2();
ea94c18b 1087 if (localChi2 > worstLocalChi2) {
1088 worstLocalChi2 = localChi2;
96ebe67e 1089 worstTrackParamAtCluster = trackParamAtCluster;
ea94c18b 1090 }
1091
96ebe67e 1092 trackParamAtCluster = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParamAtCluster);
ea94c18b 1093 }
1094
96ebe67e 1095 // Check if bad removable cluster found
1096 if (!worstTrackParamAtCluster) {
ea94c18b 1097 track->SetImproved(kTRUE);
1098 break;
1099 }
1100
b709ac13 1101 // Check whether the worst chi2 is under requirement or not
3304fa09 1102 if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
ea94c18b 1103 track->SetImproved(kTRUE);
1104 break;
1105 }
1106
96ebe67e 1107 // if the worst cluster is not removable then remove the entire track
1108 if (!worstTrackParamAtCluster->IsRemovable() && worstTrackParamAtCluster->IsAloneInChamber()) {
b709ac13 1109 fRecTracksPtr->Remove(track);
1110 fNRecTracks--;
ea94c18b 1111 break;
1112 }
1113
96ebe67e 1114 // Reset the second cluster in the same station as being not removable
1115 // or reset the second cluster in the same chamber as being alone
1116 worstChamber = worstTrackParamAtCluster->GetClusterPtr()->GetChamberId();
1117 previousTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->Before(worstTrackParamAtCluster);
1118 nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(worstTrackParamAtCluster);
1119 if (worstTrackParamAtCluster->IsAloneInChamber()) { // Worst cluster removable and alone in chamber
b709ac13 1120
1121 if (worstChamber%2 == 0) { // Modify flags in next chamber
1122
1123 nextTrackParam->SetRemovable(kFALSE);
96ebe67e 1124 if (!nextTrackParam->IsAloneInChamber()) // Make sure both clusters in second chamber are not removable anymore
1125 ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(nextTrackParam))->SetRemovable(kFALSE);
b709ac13 1126
1127 } else { // Modify flags in previous chamber
1128
1129 previousTrackParam->SetRemovable(kFALSE);
96ebe67e 1130 if (!previousTrackParam->IsAloneInChamber()) // Make sure both clusters in second chamber are not removable anymore
1131 ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->Before(previousTrackParam))->SetRemovable(kFALSE);
b709ac13 1132
1133 }
1134
96ebe67e 1135 } else { // Worst cluster not alone in its chamber
b709ac13 1136
96ebe67e 1137 if (previousTrackParam) previousChamber = previousTrackParam->GetClusterPtr()->GetChamberId();
b709ac13 1138 else previousChamber = -1;
1139
96ebe67e 1140 if (previousChamber == worstChamber) { // the second cluster on the same chamber is the previous one
b709ac13 1141
1142 previousTrackParam->SetAloneInChamber(kTRUE);
96ebe67e 1143 // transfert the removability to the second cluster
1144 if (worstTrackParamAtCluster->IsRemovable()) previousTrackParam->SetRemovable(kTRUE);
b709ac13 1145
96ebe67e 1146 } else { // the second cluster on the same chamber is the next one
b709ac13 1147
1148 nextTrackParam->SetAloneInChamber(kTRUE);
96ebe67e 1149 // transfert the removability to the second cluster
1150 if (worstTrackParamAtCluster->IsRemovable()) nextTrackParam->SetRemovable(kTRUE);
b709ac13 1151
1152 }
1153
1154 }
ea94c18b 1155
96ebe67e 1156 // Remove the worst cluster
1157 track->RemoveTrackParamAtCluster(worstTrackParamAtCluster);
ea94c18b 1158
1159 // Re-calculate track parameters
96ebe67e 1160 // - from the cluster immediately downstream the one suppressed
ea94c18b 1161 // - or from the begining - if parameters have been re-computed using the standard method (kalman parameters have been lost)
96ebe67e 1162 // - or if the removed cluster was the last one
b709ac13 1163 if (smoothed && nextTrackParam) RetracePartialTrack(*track,nextTrackParam);
ea94c18b 1164 else RetraceTrack(*track,kTRUE);
1165
1166 // Printout for debuging
1167 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1168 cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(track)+1 << " improved " << endl;
1169 }
1170
1171 }
1172
b709ac13 1173 track = nextTrack;
ea94c18b 1174 }
1175
b709ac13 1176 // compress the array in case of some tracks have been removed
1177 fRecTracksPtr->Compress();
1178
8cc77df0 1179}
1180
b0e4c598 1181 //__________________________________________________________________________
ea94c18b 1182void AliMUONTrackReconstructorK::Finalize()
b0e4c598 1183{
96ebe67e 1184 /// Update track parameters and covariances at each attached cluster
b709ac13 1185
ea94c18b 1186 AliMUONTrack *track;
96ebe67e 1187 AliMUONTrackParam *trackParamAtCluster;
ea94c18b 1188 Bool_t smoothed = kFALSE;
1189
1190 track = (AliMUONTrack*) fRecTracksPtr->First();
1191 while (track) {
1192
1193 // update track parameters (using smoother if required) if not already done
1194 if (!track->IsImproved()) {
1195 smoothed = kFALSE;
3304fa09 1196 if (AliMUONReconstructor::GetRecoParam()->UseSmoother()) smoothed = RunSmoother(*track);
96ebe67e 1197 if (!smoothed) track->UpdateCovTrackParamAtCluster();
ea94c18b 1198 }
1199
96ebe67e 1200 // copy smoothed parameters and covariances if any
1201 if (smoothed) {
ea94c18b 1202
96ebe67e 1203 trackParamAtCluster = (AliMUONTrackParam*) (track->GetTrackParamAtCluster()->First());
1204 while (trackParamAtCluster) {
1205
1206 trackParamAtCluster->SetParameters(trackParamAtCluster->GetSmoothParameters());
1207 trackParamAtCluster->SetCovariances(trackParamAtCluster->GetSmoothCovariances());
1208
1209 trackParamAtCluster = (AliMUONTrackParam*) (track->GetTrackParamAtCluster()->After(trackParamAtCluster));
ea94c18b 1210 }
1211
ea94c18b 1212 }
1213
1214 track = (AliMUONTrack*) fRecTracksPtr->After(track);
b0e4c598 1215 }
ea94c18b 1216
b0e4c598 1217}
ea94c18b 1218