MUON/AliMUONTrackReconstructor.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 /// \class AliMUONTrackReconstructor
  20 /// MUON track reconstructor using the original method
  21 ///
  22 /// This class contains as data:
  23 /// - the parameters for the track reconstruction
  24 ///
  25 /// It contains as methods, among others:
  26 /// - MakeTracks to build the tracks
  27 //-----------------------------------------------------------------------------
  28
  29 #include "AliMUONTrackReconstructor.h"
  30
  31 #include "AliMUONConstants.h"
  32 #include "AliMUONVCluster.h"
  33 #include "AliMUONVClusterServer.h"
  34 #include "AliMUONVClusterStore.h"
  35 #include "AliMUONTrack.h"
  36 #include "AliMUONTrackParam.h"
  37 #include "AliMUONTrackExtrap.h"
  38 #include "AliMUONRecoParam.h"
  39
  40 #include "AliMpArea.h"
  41
  42 #include "AliLog.h"
  43
  44 #include <TMinuit.h>
  45 #include <Riostream.h>
  46 #include <TMath.h>
  47 #include <TMatrixD.h>
  48
  49 // Functions to be minimized with Minuit
  50 void TrackChi2(Int_t &nParam, Double_t *gradient, Double_t &chi2, Double_t *param, Int_t flag);
  51
  52 /// \cond CLASSIMP
  53 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
  54 /// \endcond
  55
  56   //__________________________________________________________________________
  57 AliMUONTrackReconstructor::AliMUONTrackReconstructor(const AliMUONRecoParam* recoParam, AliMUONVClusterServer* clusterServer)
  58   : AliMUONVTrackReconstructor(recoParam,clusterServer)
  59 {
  60   /// Constructor
  61 }
  62
  63   //__________________________________________________________________________
  64 AliMUONTrackReconstructor::~AliMUONTrackReconstructor()
  65 {
  66 /// Destructor
  67 }
  68
  69   //__________________________________________________________________________
  70 void AliMUONTrackReconstructor::MakeTrackCandidates(AliMUONVClusterStore& clusterStore)
  71 {
  72   /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
  73   /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
  74   /// Good candidates are made of at least three clusters.
  75   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
  76
  77   TClonesArray *segments;
  78   AliMUONTrack *track;
  79   Int_t iCandidate = 0;
  80   Bool_t clusterFound;
  81
  82   AliDebug(1,"Enter MakeTrackCandidates");
  83
  84   // Unless we're doing combined tracking, we'll clusterize all stations at once
  85   Int_t firstChamber(0);
  86   Int_t lastChamber(9);
  87
  88   if (GetRecoParam()->CombineClusterTrackReco()) {
  89     // ... Here's the exception : ask the clustering to reconstruct
  90     // clusters *only* in station 4 and 5 for combined tracking
  91     firstChamber = 6;
  92   }
  93
  94   for (Int_t i = firstChamber; i <= lastChamber; ++i )
  95   {
  96     if (fClusterServer && GetRecoParam()->UseChamber(i)) fClusterServer->Clusterize(i, clusterStore, AliMpArea(), GetRecoParam());
  97   }
  98
  99   // Loop over stations(1..) 5 and 4 and make track candidates
 100   for (Int_t istat=4; istat>=3; istat--) {
 101
 102     // Make segments in the station
 103     segments = MakeSegmentsBetweenChambers(clusterStore, 2*istat, 2*istat+1);
 104
 105     // Loop over segments
 106     for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++)
 107     {
 108       AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 109
 110       // Transform segments to tracks and put them at the end of fRecTracksPtr
 111       track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]),GetRecoParam()->GetBendingVertexDispersion());
 112       fNRecTracks++;
 113
 114       // Look for compatible cluster(s) in the other station
 115       if (GetRecoParam()->MakeTrackCandidatesFast())
 116         clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat);
 117       else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat);
 118
 119       // Remove track if no cluster found
 120       if (!clusterFound && GetRecoParam()->RequestStation(7-istat)) {
 121         fRecTracksPtr->Remove(track);
 122         fNRecTracks--;
 123       }
 124
 125     }
 126
 127     // delete the array of segments
 128     delete segments;
 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 (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
 135   else RemoveDoubleTracks();
 136
 137   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 138
 139 }
 140
 141   //__________________________________________________________________________
 142 void AliMUONTrackReconstructor::MakeMoreTrackCandidates(AliMUONVClusterStore& clusterStore)
 143 {
 144   /// To make extra track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE):
 145   /// clustering is supposed to be already done
 146   /// Start with segments made of 1 cluster in each of the stations 4 and 5 then follow track in remaining chambers.
 147   /// Good candidates are made of at least three clusters if both station are requested (two otherwise).
 148   /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
 149
 150   TClonesArray *segments;
 151   AliMUONTrack *track;
 152   Int_t iCandidate = 0, iCurrentTrack, nCurrentTracks;
 153   Bool_t clusterFound;
 154
 155   AliDebug(1,"Enter MakeMoreTrackCandidates");
 156
 157   // Double loop over chambers in stations(1..) 4 and 5 to make track candidates
 158   for (Int_t ich1 = 6; ich1 <= 7; ich1++) {
 159     for (Int_t ich2 = 8; ich2 <= 9; ich2++) {
 160
 161       // Make segments in the station
 162       segments = MakeSegmentsBetweenChambers(clusterStore, ich1, ich2);
 163
 164       // Loop over segments
 165       for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++)
 166       {
 167         AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
 168
 169         // Transform segments to tracks and put them at the end of fRecTracksPtr
 170         iCurrentTrack = fRecTracksPtr->GetLast()+1;
 171         track = new ((*fRecTracksPtr)[iCurrentTrack]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]),GetRecoParam()->GetBendingVertexDispersion());
 172         fNRecTracks++;
 173
 174         // Look for compatible cluster(s) in the second chamber of station 5
 175         if (GetRecoParam()->MakeTrackCandidatesFast())
 176           clusterFound = FollowLinearTrackInChamber(*track, clusterStore, 17-ich2);
 177         else clusterFound = FollowTrackInChamber(*track, clusterStore, 17-ich2);
 178
 179         // skip the original track in case it has been removed
 180         if (GetRecoParam()->TrackAllTracks() && clusterFound) iCurrentTrack++;
 181
 182         // loop over every new tracks
 183         nCurrentTracks = fRecTracksPtr->GetLast()+1;
 184         while (iCurrentTrack < nCurrentTracks) {
 185           track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iCurrentTrack);
 186
 187           // Look for compatible cluster(s) in the second chamber of station 4
 188           if (GetRecoParam()->MakeTrackCandidatesFast())
 189             FollowLinearTrackInChamber(*track, clusterStore, 13-ich1);
 190           else FollowTrackInChamber(*track, clusterStore, 13-ich1);
 191
 192           iCurrentTrack++;
 193         }
 194
 195       }
 196
 197       // delete the array of segments
 198       delete segments;
 199     }
 200   }
 201
 202   fRecTracksPtr->Compress(); // this is essential before checking tracks
 203
 204   // Keep all different tracks or only the best ones as required
 205   if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks();
 206   else RemoveDoubleTracks();
 207
 208   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
 209
 210 }
 211
 212   //__________________________________________________________________________
 213 void AliMUONTrackReconstructor::FollowTracks(AliMUONVClusterStore& clusterStore)
 214 {
 215   /// Follow tracks in stations(1..) 3, 2 and 1
 216   AliDebug(1,"Enter FollowTracks");
 217
 218   AliMUONTrack *track, *nextTrack;
 219   AliMUONTrackParam *trackParam, *nextTrackParam;
 220   Int_t currentNRecTracks;
 221
 222   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() *
 223                        GetRecoParam()->GetSigmaCutForTracking();
 224
 225   for (Int_t station = 2; station >= 0; station--) {
 226
 227     // Save the actual number of reconstructed track in case of
 228     // tracks are added or suppressed during the tracking procedure
 229     // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
 230     currentNRecTracks = fNRecTracks;
 231
 232     for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
 233       AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
 234
 235       track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
 236
 237       // Fit the track:
 238       // Do not take into account the multiple scattering to speed up the fit
 239       // Calculate the track parameter covariance matrix
 240       // If there is no cluster out of station 4 or 5 then use the vertex to better constrain the fit
 241       if (((AliMUONTrackParam*) track->GetTrackParamAtCluster()->First())->GetClusterPtr()->GetChamberId() > 5)
 242         Fit(*track, kFALSE, kTRUE, kTRUE);
 243       else Fit(*track, kFALSE, kFALSE, kTRUE);
 244
 245       // remove track with absolute bending momentum out of limits
 246       // or if the normalized chi2 is too high
 247       Double_t bendingMomentum = TMath::Abs(1. / ((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First())->GetInverseBendingMomentum());
 248       if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
 249           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum() ||
 250           track->GetNormalizedChi2() > sigmaCut2) {
 251         fRecTracksPtr->Remove(track);
 252         fNRecTracks--;
 253         continue;
 254       }
 255
 256       // save parameters from fit into smoothed parameters to complete track afterward
 257       if (GetRecoParam()->ComplementTracks()) {
 258
 259         if (station==2) { // save track parameters on stations 4 and 5
 260
 261           // extrapolate track parameters and covariances at each cluster
 262           track->UpdateCovTrackParamAtCluster();
 263
 264           // save them
 265           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 266           while (trackParam) {
 267             trackParam->SetSmoothParameters(trackParam->GetParameters());
 268             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 269             trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
 270           }
 271
 272         } else { // or save track parameters on last station only
 273
 274           trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 275           if (trackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) {
 276
 277             // save parameters from fit
 278             trackParam->SetSmoothParameters(trackParam->GetParameters());
 279             trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 280
 281             // save parameters extrapolated to the second chamber of the same station if it has been hit
 282             nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
 283             if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) {
 284
 285               // reset parameters and covariances
 286               nextTrackParam->SetParameters(trackParam->GetParameters());
 287               nextTrackParam->SetZ(trackParam->GetZ());
 288               nextTrackParam->SetCovariances(trackParam->GetCovariances());
 289
 290               // extrapolate them to the z of the corresponding cluster
 291               AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ());
 292
 293               // save them
 294               nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
 295               nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
 296
 297             }
 298
 299           }
 300
 301         }
 302
 303       }
 304
 305       // Look for compatible cluster(s) in station(0..) "station"
 306       if (!FollowTrackInStation(*track, clusterStore, station)) {
 307
 308         // Try to recover track if required
 309         if (GetRecoParam()->RecoverTracks()) {
 310
 311           // work on a copy of the track if this station is not required
 312           // to keep the case where no cluster is reconstructed as a possible candidate
 313           if (!GetRecoParam()->RequestStation(station)) {
 314             track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*track);
 315             fNRecTracks++;
 316           }
 317
 318           // try to recover
 319           if (!RecoverTrack(*track, clusterStore, station)) {
 320             // remove track if no cluster found
 321             fRecTracksPtr->Remove(track);
 322             fNRecTracks--;
 323           }
 324
 325         } else if (GetRecoParam()->RequestStation(station)) {
 326           // remove track if no cluster found
 327           fRecTracksPtr->Remove(track);
 328           fNRecTracks--;
 329         }
 330
 331       }
 332
 333     }
 334
 335     // Compress fRecTracksPtr for the next step
 336     fRecTracksPtr->Compress();
 337
 338     // Keep only the best tracks if required
 339     if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
 340
 341   }
 342
 343   // Last fit of track candidates with all stations
 344   // Take into account the multiple scattering and remove bad tracks
 345   Int_t trackIndex = -1;
 346   track = (AliMUONTrack*) fRecTracksPtr->First();
 347   while (track) {
 348
 349     trackIndex++;
 350     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
 351
 352     Fit(*track, kTRUE, kFALSE, kTRUE);
 353
 354     // Printout for debuging
 355     if (AliLog::GetGlobalDebugLevel() >= 3) {
 356       cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
 357       track->RecursiveDump();
 358     }
 359
 360     // Remove the track if the normalized chi2 is too high
 361     if (track->GetNormalizedChi2() > sigmaCut2) {
 362       fRecTracksPtr->Remove(track);
 363       fNRecTracks--;
 364       track = nextTrack;
 365       continue;
 366     }
 367
 368     // save parameters from fit into smoothed parameters to complete track afterward
 369     if (GetRecoParam()->ComplementTracks()) {
 370
 371       trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
 372       if (trackParam->GetClusterPtr()->GetChamberId() < 2) {
 373
 374         // save parameters from fit
 375         trackParam->SetSmoothParameters(trackParam->GetParameters());
 376         trackParam->SetSmoothCovariances(trackParam->GetCovariances());
 377
 378         // save parameters extrapolated to the second chamber of the same station if it has been hit
 379         nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam);
 380         if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2) {
 381
 382           // reset parameters and covariances
 383           nextTrackParam->SetParameters(trackParam->GetParameters());
 384           nextTrackParam->SetZ(trackParam->GetZ());
 385           nextTrackParam->SetCovariances(trackParam->GetCovariances());
 386
 387           // extrapolate them to the z of the corresponding cluster
 388           AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ());
 389
 390           // save them
 391           nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters());
 392           nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances());
 393
 394         }
 395
 396       }
 397
 398     }
 399
 400     track = nextTrack;
 401
 402   }
 403
 404   fRecTracksPtr->Compress();
 405
 406 }
 407
 408   //__________________________________________________________________________
 409 Bool_t AliMUONTrackReconstructor::FollowTrackInChamber(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextChamber)
 410 {
 411   /// Follow trackCandidate in chamber(0..) nextChamber and search for compatible cluster(s)
 412   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 413   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 414   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 415   ///         Remove the obsolete "trackCandidate" at the end.
 416   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 417   AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1));
 418
 419   Double_t chi2WithOneCluster = 1.e10;
 420   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 421                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 422   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
 423   Bool_t foundOneCluster = kFALSE;
 424   AliMUONTrack *newTrack = 0x0;
 425   AliMUONVCluster *cluster;
 426   AliMUONTrackParam extrapTrackParam;
 427   AliMUONTrackParam extrapTrackParamAtCh;
 428   AliMUONTrackParam extrapTrackParamAtCluster;
 429   AliMUONTrackParam bestTrackParamAtCluster;
 430
 431   // Get track parameters according to the propagation direction
 432   if (nextChamber > 7) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 433   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 434
 435   // Printout for debuging
 436   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 437     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 438     extrapTrackParamAtCh.GetParameters().Print();
 439     extrapTrackParamAtCh.GetCovariances().Print();
 440   }
 441
 442   // Add MCS effect
 443   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 444
 445   // Add MCS in the missing chamber(s) if any
 446   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
 447   while (currentChamber > nextChamber + 1) {
 448     // extrapolation to the missing chamber
 449     currentChamber--;
 450     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber));
 451     // add MCS effect
 452     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 453   }
 454
 455   //Extrapolate trackCandidate to chamber
 456   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(nextChamber));
 457
 458   // Printout for debuging
 459   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 460     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(nextChamber)<<":"<<endl;
 461     extrapTrackParamAtCh.GetParameters().Print();
 462     extrapTrackParamAtCh.GetCovariances().Print();
 463   }
 464
 465   // Printout for debuging
 466   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 467     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << nextChamber+1 << endl;
 468   }
 469
 470   // Ask the clustering to reconstruct new clusters around the track position in the current chamber
 471   // except for station 4 and 5 that are already entirely clusterized
 472   if (GetRecoParam()->CombineClusterTrackReco()) {
 473     if (nextChamber < 6) AskForNewClustersInChamber(extrapTrackParamAtCh, clusterStore, nextChamber);
 474   }
 475
 476   // Create iterators to loop over clusters in both chambers
 477   TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber));
 478
 479   // look for cluster in chamber
 480   while ( ( cluster = static_cast<AliMUONVCluster*>(next()) ) ) {
 481
 482     // try to add the current cluster fast
 483     if (!TryOneClusterFast(extrapTrackParamAtCh, cluster)) continue;
 484
 485     // try to add the current cluster accuratly
 486     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, cluster, extrapTrackParamAtCluster);
 487
 488     // if good chi2 then consider to add cluster
 489     if (chi2WithOneCluster < maxChi2WithOneCluster) {
 490       foundOneCluster = kTRUE;
 491
 492       // Printout for debuging
 493       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 494         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << nextChamber+1
 495         << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 496         cluster->Print();
 497       }
 498
 499       if (GetRecoParam()->TrackAllTracks()) {
 500         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 501         newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 502         UpdateTrack(*newTrack,extrapTrackParamAtCluster);
 503         fNRecTracks++;
 504
 505         // Printout for debuging
 506         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 507           cout << "FollowTrackInStation: added one cluster in chamber(1..): " << nextChamber+1 << endl;
 508           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 509         }
 510
 511       } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
 512         // keep track of the best single cluster except if a couple of clusters has already been found
 513         bestChi2WithOneCluster = chi2WithOneCluster;
 514         bestTrackParamAtCluster = extrapTrackParamAtCluster;
 515       }
 516
 517     }
 518
 519   }
 520
 521   // fill out the best track if required else clean up the fRecTracksPtr array
 522   if (!GetRecoParam()->TrackAllTracks()) {
 523     if (foundOneCluster) {
 524       UpdateTrack(trackCandidate,bestTrackParamAtCluster);
 525
 526       // Printout for debuging
 527       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 528         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl;
 529         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 530       }
 531
 532     } else return kFALSE;
 533
 534   } else if (foundOneCluster) {
 535
 536     // remove obsolete track
 537     fRecTracksPtr->Remove(&trackCandidate);
 538     fNRecTracks--;
 539
 540   } else return kFALSE;
 541
 542   return kTRUE;
 543
 544 }
 545
 546   //__________________________________________________________________________
 547 Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
 548 {
 549   /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s)
 550   /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
 551   /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
 552   ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
 553   ///         Remove the obsolete "trackCandidate" at the end.
 554   /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority.
 555   AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
 556
 557   // Order the chamber according to the propagation direction (tracking starts with chamber 2):
 558   // - nextStation == station(1...) 5 => forward propagation
 559   // - nextStation < station(1...) 5 => backward propagation
 560   Int_t ch1, ch2;
 561   if (nextStation==4) {
 562     ch1 = 2*nextStation+1;
 563     ch2 = 2*nextStation;
 564   } else {
 565     ch1 = 2*nextStation;
 566     ch2 = 2*nextStation+1;
 567   }
 568
 569   Double_t chi2WithOneCluster = 1.e10;
 570   Double_t chi2WithTwoClusters = 1.e10;
 571   Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
 572                                         GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
 573   Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() *
 574                                          GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2
 575   Double_t bestChi2WithOneCluster = maxChi2WithOneCluster;
 576   Double_t bestChi2WithTwoClusters = maxChi2WithTwoClusters;
 577   Bool_t foundOneCluster = kFALSE;
 578   Bool_t foundTwoClusters = kFALSE;
 579   AliMUONTrack *newTrack = 0x0;
 580   AliMUONVCluster *clusterCh1, *clusterCh2;
 581   AliMUONTrackParam extrapTrackParam;
 582   AliMUONTrackParam extrapTrackParamAtCh;
 583   AliMUONTrackParam extrapTrackParamAtCluster1;
 584   AliMUONTrackParam extrapTrackParamAtCluster2;
 585   AliMUONTrackParam bestTrackParamAtCluster1;
 586   AliMUONTrackParam bestTrackParamAtCluster2;
 587
 588   Int_t nClusters = clusterStore.GetSize();
 589   Bool_t *clusterCh1Used = new Bool_t[nClusters];
 590   for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE;
 591   Int_t iCluster1;
 592
 593   // Get track parameters according to the propagation direction
 594   if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last();
 595   else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First();
 596
 597   // Printout for debuging
 598   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 599     cout<<endl<<"Track parameters and covariances at first cluster:"<<endl;
 600     extrapTrackParamAtCh.GetParameters().Print();
 601     extrapTrackParamAtCh.GetCovariances().Print();
 602   }
 603
 604   // Add MCS effect
 605   AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 606
 607   // Add MCS in the missing chamber(s) if any
 608   Int_t currentChamber = extrapTrackParamAtCh.GetClusterPtr()->GetChamberId();
 609   while (ch1 < ch2 && currentChamber > ch2 + 1) {
 610     // extrapolation to the missing chamber
 611     currentChamber--;
 612     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber));
 613     // add MCS effect
 614     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 615   }
 616
 617   //Extrapolate trackCandidate to chamber "ch2"
 618   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2));
 619
 620   // Printout for debuging
 621   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 622     cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch2)<<":"<<endl;
 623     extrapTrackParamAtCh.GetParameters().Print();
 624     extrapTrackParamAtCh.GetCovariances().Print();
 625   }
 626
 627   // Printout for debuging
 628   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 629     cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl;
 630   }
 631
 632   // Ask the clustering to reconstruct new clusters around the track position in the current station
 633   // except for station 4 and 5 that are already entirely clusterized
 634   if (GetRecoParam()->CombineClusterTrackReco()) {
 635     if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation);
 636   }
 637
 638   // Create iterators to loop over clusters in both chambers
 639   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
 640   TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2));
 641
 642   // look for candidates in chamber 2
 643   while ( ( clusterCh2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
 644
 645     // try to add the current cluster fast
 646     if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue;
 647
 648     // try to add the current cluster accuratly
 649     chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2);
 650
 651     // if good chi2 then try to attach a cluster in the other chamber too
 652     if (chi2WithOneCluster < maxChi2WithOneCluster) {
 653       Bool_t foundSecondCluster = kFALSE;
 654
 655       // Printout for debuging
 656       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 657         cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1
 658              << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 659         clusterCh2->Print();
 660         cout << "                      look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl;
 661       }
 662
 663       // add MCS effect for next step
 664       AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(),1.);
 665
 666       // copy new track parameters for next step
 667       extrapTrackParam = extrapTrackParamAtCluster2;
 668
 669       //Extrapolate track parameters to chamber "ch1"
 670       AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1));
 671
 672       // reset cluster iterator of chamber 1
 673       nextInCh1.Reset();
 674       iCluster1 = -1;
 675
 676       // look for second candidates in chamber 1
 677       while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 678         iCluster1++;
 679
 680         // try to add the current cluster fast
 681         if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue;
 682
 683         // try to add the current cluster accuratly
 684         chi2WithTwoClusters = TryTwoClusters(extrapTrackParamAtCluster2, clusterCh1, extrapTrackParamAtCluster1);
 685
 686         // if good chi2 then create a new track by adding the 2 clusters to the "trackCandidate"
 687         if (chi2WithTwoClusters < maxChi2WithTwoClusters) {
 688           foundSecondCluster = kTRUE;
 689           foundTwoClusters = kTRUE;
 690
 691           // Printout for debuging
 692           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 693             cout << "FollowTrackInStation: found second cluster in chamber(1..): " << ch1+1
 694                  << " (Global Chi2 = " << chi2WithTwoClusters << ")" << endl;
 695             clusterCh1->Print();
 696           }
 697
 698           if (GetRecoParam()->TrackAllTracks()) {
 699             // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters
 700             newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 701             UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2);
 702             fNRecTracks++;
 703
 704             // Tag clusterCh1 as used
 705             clusterCh1Used[iCluster1] = kTRUE;
 706
 707             // Printout for debuging
 708             if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 709               cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl;
 710               if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 711             }
 712
 713           } else if (chi2WithTwoClusters < bestChi2WithTwoClusters) {
 714             // keep track of the best couple of clusters
 715             bestChi2WithTwoClusters = chi2WithTwoClusters;
 716             bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 717             bestTrackParamAtCluster2 = extrapTrackParamAtCluster2;
 718           }
 719
 720         }
 721
 722       }
 723
 724       // if no clusterCh1 found then consider to add clusterCh2 only
 725       if (!foundSecondCluster) {
 726         foundOneCluster = kTRUE;
 727
 728         if (GetRecoParam()->TrackAllTracks()) {
 729           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 730           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 731           UpdateTrack(*newTrack,extrapTrackParamAtCluster2);
 732           fNRecTracks++;
 733
 734           // Printout for debuging
 735           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 736             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl;
 737             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 738           }
 739
 740         } else if (!foundTwoClusters && chi2WithOneCluster < bestChi2WithOneCluster) {
 741           // keep track of the best single cluster except if a couple of clusters has already been found
 742           bestChi2WithOneCluster = chi2WithOneCluster;
 743           bestTrackParamAtCluster1 = extrapTrackParamAtCluster2;
 744         }
 745
 746       }
 747
 748     }
 749
 750   }
 751
 752   // look for candidates in chamber 1 not already attached to a track
 753   // if we want to keep all possible tracks or if no good couple of clusters has been found
 754   if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) {
 755
 756     // add MCS effect for next step
 757     AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.);
 758
 759     //Extrapolate trackCandidate to chamber "ch1"
 760     AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1));
 761
 762     // Printout for debuging
 763     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
 764       cout<<endl<<"Track parameters and covariances at first cluster extrapolated to z = "<<AliMUONConstants::DefaultChamberZ(ch1)<<":"<<endl;
 765       extrapTrackParamAtCh.GetParameters().Print();
 766       extrapTrackParamAtCh.GetCovariances().Print();
 767     }
 768
 769     // Printout for debuging
 770     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 771       cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl;
 772     }
 773
 774     // reset cluster iterator of chamber 1
 775     nextInCh1.Reset();
 776     iCluster1 = -1;
 777
 778     // look for second candidates in chamber 1
 779     while ( ( clusterCh1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
 780       iCluster1++;
 781
 782       if (clusterCh1Used[iCluster1]) continue; // Skip cluster already used
 783
 784       // try to add the current cluster fast
 785       if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue;
 786
 787       // try to add the current cluster accuratly
 788       chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1);
 789
 790       // if good chi2 then consider to add clusterCh1
 791       // We do not try to attach a cluster in the other chamber too since it has already been done above
 792       if (chi2WithOneCluster < maxChi2WithOneCluster) {
 793         foundOneCluster = kTRUE;
 794
 795         // Printout for debuging
 796         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 797           cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1
 798                << " (Chi2 = " << chi2WithOneCluster << ")" << endl;
 799           clusterCh1->Print();
 800         }
 801
 802         if (GetRecoParam()->TrackAllTracks()) {
 803           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster
 804           newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate);
 805           UpdateTrack(*newTrack,extrapTrackParamAtCluster1);
 806           fNRecTracks++;
 807
 808           // Printout for debuging
 809           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 810             cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl;
 811             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 812           }
 813
 814         } else if (chi2WithOneCluster < bestChi2WithOneCluster) {
 815           // keep track of the best single cluster except if a couple of clusters has already been found
 816           bestChi2WithOneCluster = chi2WithOneCluster;
 817           bestTrackParamAtCluster1 = extrapTrackParamAtCluster1;
 818         }
 819
 820       }
 821
 822     }
 823
 824   }
 825
 826   // fill out the best track if required else clean up the fRecTracksPtr array
 827   if (!GetRecoParam()->TrackAllTracks()) {
 828     if (foundTwoClusters) {
 829       UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2);
 830
 831       // Printout for debuging
 832       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 833         cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl;
 834         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 835       }
 836
 837     } else if (foundOneCluster) {
 838       UpdateTrack(trackCandidate,bestTrackParamAtCluster1);
 839
 840       // Printout for debuging
 841       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
 842         cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl;
 843         if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
 844       }
 845
 846     } else {
 847       delete [] clusterCh1Used;
 848       return kFALSE;
 849     }
 850
 851   } else if (foundOneCluster || foundTwoClusters) {
 852
 853     // remove obsolete track
 854     fRecTracksPtr->Remove(&trackCandidate);
 855     fNRecTracks--;
 856
 857   } else {
 858     delete [] clusterCh1Used;
 859     return kFALSE;
 860   }
 861
 862   delete [] clusterCh1Used;
 863   return kTRUE;
 864
 865 }
 866
 867   //__________________________________________________________________________
 868 Double_t AliMUONTrackReconstructor::TryTwoClusters(const AliMUONTrackParam &trackParamAtCluster1, AliMUONVCluster* cluster2,
 869                                                    AliMUONTrackParam &trackParamAtCluster2)
 870 {
 871 /// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix):
 872 /// return the corresponding Chi2 accounting for covariances between the 2 clusters
 873 /// return trackParamAtCluster1 & 2
 874
 875   // extrapolate track parameters at the z position of the second cluster (no need to extrapolate the covariances)
 876   trackParamAtCluster2.SetParameters(trackParamAtCluster1.GetParameters());
 877   trackParamAtCluster2.SetZ(trackParamAtCluster1.GetZ());
 878   AliMUONTrackExtrap::ExtrapToZ(&trackParamAtCluster2, cluster2->GetZ());
 879
 880   // set pointer to cluster2 into trackParamAtCluster2
 881   trackParamAtCluster2.SetClusterPtr(cluster2);
 882
 883   // Set differences between track and the 2 clusters in the bending and non bending directions
 884   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
 885   TMatrixD dPos(4,1);
 886   dPos(0,0) = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor();
 887   dPos(1,0) = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor();
 888   dPos(2,0) = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor();
 889   dPos(3,0) = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor();
 890
 891   // Calculate the error matrix from the track parameter covariances at first cluster
 892   TMatrixD error(4,4);
 893   error.Zero();
 894   if (trackParamAtCluster1.CovariancesExist()) {
 895     // Save track parameters at first cluster
 896     TMatrixD paramAtCluster1Save(trackParamAtCluster1.GetParameters());
 897
 898     // Save track coordinates at second cluster
 899     Double_t nonBendingCoor2 = trackParamAtCluster2.GetNonBendingCoor();
 900     Double_t bendingCoor2    = trackParamAtCluster2.GetBendingCoor();
 901
 902     // copy track parameters at first cluster for jacobian calculation
 903     AliMUONTrackParam trackParam(trackParamAtCluster1);
 904
 905     // add MCS effect to the covariance matrix at first cluster
 906     AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
 907
 908     // Get the pointer to the parameter covariance matrix at first cluster
 909     const TMatrixD& kParamCov = trackParam.GetCovariances();
 910
 911     // Calculate the jacobian related to the transformation between track parameters
 912     // at first cluster and track coordinates at the 2 cluster z-positions
 913     TMatrixD jacob(4,5);
 914     jacob.Zero();
 915     // first derivative at the first cluster:
 916     jacob(0,0) = 1.; // dx1/dx
 917     jacob(1,2) = 1.; // dy1/dy
 918     // first derivative at the second cluster:
 919     TMatrixD dParam(5,1);
 920     for (Int_t i=0; i<5; i++) {
 921       // Skip jacobian calculation for parameters with no associated error
 922       if (kParamCov(i,i) == 0.) continue;
 923       // Small variation of parameter i only
 924       for (Int_t j=0; j<5; j++) {
 925         if (j==i) {
 926           dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
 927           if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramAtCluster1Save(4,0)); // variation always in the same direction
 928         } else dParam(j,0) = 0.;
 929       }
 930
 931       // Set new track parameters at first cluster
 932       trackParam.SetParameters(paramAtCluster1Save);
 933       trackParam.AddParameters(dParam);
 934       trackParam.SetZ(cluster1->GetZ());
 935
 936       // Extrapolate new track parameters to the z position of the second cluster
 937       AliMUONTrackExtrap::ExtrapToZ(&trackParam,cluster2->GetZ());
 938
 939       // Calculate the jacobian
 940       jacob(2,i) = (trackParam.GetNonBendingCoor() - nonBendingCoor2) / dParam(i,0); // dx2/dParami
 941       jacob(3,i) = (trackParam.GetBendingCoor()    - bendingCoor2   ) / dParam(i,0); // dy2/dParami
 942     }
 943
 944     // Calculate the error matrix
 945     TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov);
 946     error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob);
 947   }
 948
 949   // Add cluster resolution to the error matrix
 950   error(0,0) += cluster1->GetErrX2();
 951   error(1,1) += cluster1->GetErrY2();
 952   error(2,2) += cluster2->GetErrX2();
 953   error(3,3) += cluster2->GetErrY2();
 954
 955   // invert the error matrix for Chi2 calculation
 956   if (error.Determinant() != 0) {
 957     error.Invert();
 958   } else {
 959     AliWarning(" Determinant error=0");
 960     return 1.e10;
 961   }
 962
 963   // Compute the Chi2 value
 964   TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error);
 965   TMatrixD result(tmp2,TMatrixD::kMult,dPos);
 966
 967   return result(0,0);
 968
 969 }
 970
 971   //__________________________________________________________________________
 972 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster)
 973 {
 974   /// Add 1 cluster to the track candidate
 975   /// Update chi2 of the track
 976
 977   // Compute local chi2
 978   AliMUONVCluster* cluster = trackParamAtCluster.GetClusterPtr();
 979   Double_t deltaX = trackParamAtCluster.GetNonBendingCoor() - cluster->GetX();
 980   Double_t deltaY = trackParamAtCluster.GetBendingCoor() - cluster->GetY();
 981   Double_t localChi2 = deltaX*deltaX / cluster->GetErrX2() +
 982                        deltaY*deltaY / cluster->GetErrY2();
 983
 984   // Flag cluster as being not removable
 985   if (GetRecoParam()->RequestStation(cluster->GetChamberId()/2))
 986     trackParamAtCluster.SetRemovable(kFALSE);
 987   else trackParamAtCluster.SetRemovable(kTRUE);
 988   trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used
 989
 990   // Update the chi2 of the new track
 991   track.SetGlobalChi2(track.GetGlobalChi2() + localChi2);
 992
 993   // Update TrackParamAtCluster
 994   track.AddTrackParamAtCluster(trackParamAtCluster,*cluster);
 995
 996 }
 997
 998   //__________________________________________________________________________
 999 void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2)
1000 {
1001   /// Add 2 clusters to the track candidate
1002   /// Update track and local chi2
1003
1004   // Update local chi2 at first cluster
1005   AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr();
1006   Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX();
1007   Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY();
1008   Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() +
1009                                  deltaY*deltaY / cluster1->GetErrY2();
1010   trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1);
1011
1012   // Flag first cluster as being removable
1013   trackParamAtCluster1.SetRemovable(kTRUE);
1014
1015   // Update local chi2 at second cluster
1016   AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr();
1017   deltaX = trackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX();
1018   deltaY = trackParamAtCluster2.GetBendingCoor() - cluster2->GetY();
1019   Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() +
1020                                  deltaY*deltaY / cluster2->GetErrY2();
1021   trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2);
1022
1023   // Flag first cluster as being removable
1024   trackParamAtCluster2.SetRemovable(kTRUE);
1025
1026   // Update the chi2 of the new track
1027   track.SetGlobalChi2(track.GetGlobalChi2() + localChi2AtCluster1 + localChi2AtCluster2);
1028
1029   // Update TrackParamAtCluster
1030   track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1);
1031   track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2);
1032
1033 }
1034
1035   //__________________________________________________________________________
1036 Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation)
1037 {
1038   /// Try to recover the track candidate in the next station
1039   /// by removing the worst of the two clusters attached in the current station
1040   /// Return kTRUE if recovering succeeds
1041   AliDebug(1,"Enter RecoverTrack");
1042
1043   // Do not try to recover track until we have attached cluster(s) on station(1..) 3
1044   if (nextStation > 1) return kFALSE;
1045
1046   Int_t worstClusterNumber = -1;
1047   Double_t localChi2, worstLocalChi2 = -1.;
1048
1049   // Look for the cluster to remove
1050   for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) {
1051     AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber);
1052
1053     // check if current cluster is in the previous station
1054     if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break;
1055
1056     // check if current cluster is removable
1057     if (!trackParamAtCluster->IsRemovable()) return kFALSE;
1058
1059     // reset the current cluster as beig not removable if it is on a required station
1060     if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE);
1061
1062     // Pick up cluster with the worst chi2
1063     localChi2 = trackParamAtCluster->GetLocalChi2();
1064     if (localChi2 > worstLocalChi2) {
1065       worstLocalChi2 = localChi2;
1066       worstClusterNumber = clusterNumber;
1067     }
1068
1069   }
1070
1071   // check if worst cluster found
1072   if (worstClusterNumber < 0) return kFALSE;
1073
1074   // Remove the worst cluster
1075   trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber));
1076
1077   // Re-fit the track:
1078   // Do not take into account the multiple scattering to speed up the fit
1079   // Calculate the track parameter covariance matrix
1080   Fit(trackCandidate, kFALSE, kFALSE, kTRUE);
1081
1082   // Look for new cluster(s) in next station
1083   return FollowTrackInStation(trackCandidate,clusterStore,nextStation);
1084
1085 }
1086
1087   //__________________________________________________________________________
1088 void AliMUONTrackReconstructor::SetVertexErrXY2ForFit(AliMUONTrack &trackCandidate)
1089 {
1090   /// Compute the vertex resolution square from natural vertex dispersion and
1091   /// multiple scattering effets according to trackCandidate path in absorber
1092   /// It is necessary to account for multiple scattering effects here instead of during the fit of
1093   /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
1094   AliDebug(1,"Enter SetVertexForFit");
1095
1096   Double_t nonBendingReso2 = GetRecoParam()->GetNonBendingVertexDispersion() *
1097                              GetRecoParam()->GetNonBendingVertexDispersion();
1098   Double_t bendingReso2 = GetRecoParam()->GetBendingVertexDispersion() *
1099                           GetRecoParam()->GetBendingVertexDispersion();
1100
1101   // add multiple scattering effets
1102   AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtCluster()->First())));
1103   paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
1104   AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
1105   const TMatrixD& kParamCov = paramAtVertex.GetCovariances();
1106   nonBendingReso2 += kParamCov(0,0);
1107   bendingReso2 += kParamCov(2,2);
1108
1109   // Set the vertex resolution square
1110   trackCandidate.SetVertexErrXY2(nonBendingReso2,bendingReso2);
1111 }
1112
1113   //__________________________________________________________________________
1114 void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t fitWithVertex, Bool_t calcCov)
1115 {
1116   /// Fit the track
1117   /// w/wo multiple Coulomb scattering according to "includeMCS".
1118   /// w/wo constraining the vertex according to "fitWithVertex".
1119   /// calculating or not the covariance matrix according to "calcCov".
1120   AliDebug(1,"Enter Fit");
1121
1122   Double_t benC, errorParam, invBenP, nonBenC, x, y;
1123   AliMUONTrackParam *trackParam;
1124   Double_t arg[1], fedm, errdef, globalChi2;
1125   Int_t npari, nparx;
1126   Int_t status, covStatus;
1127
1128   // Instantiate gMinuit if not already done
1129   if (!gMinuit) gMinuit = new TMinuit(6);
1130   // Clear MINUIT parameters
1131   gMinuit->mncler();
1132   // Give the fitted track to MINUIT
1133   gMinuit->SetObjectFit(&track);
1134   if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
1135     // Define print level
1136     arg[0] = 1;
1137     gMinuit->mnexcm("SET PRI", arg, 1, status);
1138     // Print covariance matrix
1139     gMinuit->mnexcm("SHO COV", arg, 0, status);
1140   } else {
1141     arg[0] = -1;
1142     gMinuit->mnexcm("SET PRI", arg, 1, status);
1143   }
1144   // No warnings
1145   gMinuit->mnexcm("SET NOW", arg, 0, status);
1146   // Define strategy
1147   //arg[0] = 2;
1148   //gMinuit->mnexcm("SET STR", arg, 1, status);
1149
1150   // set flag w/wo multiple scattering according to "includeMCS"
1151   track.FitWithMCS(includeMCS);
1152   if (includeMCS) {
1153     // compute cluster weights only once
1154     if (!track.ComputeClusterWeights()) {
1155       AliWarning("cannot take into account the multiple scattering effects");
1156       track.FitWithMCS(kFALSE);
1157     }
1158   }
1159
1160   track.FitWithVertex(fitWithVertex);
1161   if (fitWithVertex) SetVertexErrXY2ForFit(track);
1162
1163   // Set fitting function
1164   gMinuit->SetFCN(TrackChi2);
1165
1166   // Set fitted parameters (!! The order is very important for the covariance matrix !!)
1167   // Mandatory limits to avoid NaN values of parameters
1168   trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First());
1169   Double_t maxIBM = 1. / GetRecoParam()->GetMinBendingMomentum();
1170   gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
1171   gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -1., 1., status);
1172   gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
1173   gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -1.5, 1.5, status);
1174   gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -maxIBM, maxIBM, status);
1175
1176   // minimization
1177   gMinuit->mnexcm("MIGRAD", arg, 0, status);
1178
1179   // Calculate the covariance matrix more accurately if required
1180   if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
1181
1182   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
1183   gMinuit->GetParameter(0, x, errorParam);
1184   trackParam->SetNonBendingCoor(x);
1185   gMinuit->GetParameter(1, nonBenC, errorParam);
1186   trackParam->SetNonBendingSlope(nonBenC);
1187   gMinuit->GetParameter(2, y, errorParam);
1188   trackParam->SetBendingCoor(y);
1189   gMinuit->GetParameter(3, benC, errorParam);
1190   trackParam->SetBendingSlope(benC);
1191   gMinuit->GetParameter(4, invBenP, errorParam);
1192   trackParam->SetInverseBendingMomentum(invBenP);
1193
1194   // global result of the fit
1195   gMinuit->mnstat(globalChi2, fedm, errdef, npari, nparx, covStatus);
1196   track.SetGlobalChi2(globalChi2);
1197
1198   // Get the covariance matrix if required
1199   if (calcCov) {
1200     // Covariance matrix according to HESSE status
1201     // If problem then keep only the diagonal terms (variances)
1202     Double_t matrix[5][5];
1203     gMinuit->mnemat(&matrix[0][0],5);
1204     if (covStatus == 3) trackParam->SetCovariances(matrix);
1205     else trackParam->SetVariances(matrix);
1206   } else trackParam->DeleteCovariances();
1207
1208 }
1209
1210   //__________________________________________________________________________
1211 void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/)
1212 {
1213   /// Return the "Chi2" to be minimized with Minuit for track fitting.
1214   /// Assumes that the trackParamAtCluster are sorted according to increasing Z.
1215   /// Track parameters at each cluster are updated accordingly.
1216   /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()".
1217   /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()".
1218
1219   AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
1220   AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtCluster()->First();
1221   Double_t dX, dY;
1222   chi2 = 0.; // initialize chi2
1223
1224   // update track parameters
1225   trackParamAtCluster->SetNonBendingCoor(param[0]);
1226   trackParamAtCluster->SetNonBendingSlope(param[1]);
1227   trackParamAtCluster->SetBendingCoor(param[2]);
1228   trackParamAtCluster->SetBendingSlope(param[3]);
1229   trackParamAtCluster->SetInverseBendingMomentum(param[4]);
1230   trackBeingFitted->UpdateTrackParamAtCluster();
1231
1232   // Take the vertex into account in the fit if required
1233   if (trackBeingFitted->FitWithVertex()) {
1234     Double_t nonBendingReso2,bendingReso2;
1235     trackBeingFitted->GetVertexErrXY2(nonBendingReso2,bendingReso2);
1236     if (nonBendingReso2 == 0. || bendingReso2 == 0.) chi2 += 1.e10;
1237     else {
1238       AliMUONTrackParam paramAtVertex(*trackParamAtCluster);
1239       AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.); // vextex position = (0,0,0)
1240       dX = paramAtVertex.GetNonBendingCoor();
1241       dY = paramAtVertex.GetBendingCoor();
1242       chi2 += dX * dX / nonBendingReso2 + dY * dY / bendingReso2;
1243     }
1244   }
1245
1246   // compute chi2 w/wo multiple scattering
1247   chi2 += trackBeingFitted->ComputeGlobalChi2(trackBeingFitted->FitWithMCS());
1248
1249 }
1250
1251   //__________________________________________________________________________
1252 void AliMUONTrackReconstructor::ComplementTracks(const AliMUONVClusterStore& clusterStore)
1253 {
1254   /// Complete tracks by adding missing clusters (if there is an overlap between
1255   /// two detection elements, the track may have two clusters in the same chamber)
1256   /// Re-fit track parameters and covariances
1257   AliDebug(1,"Enter ComplementTracks");
1258
1259   Int_t chamberId, detElemId;
1260   Double_t chi2OfCluster, bestChi2OfCluster;
1261   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() *
1262                        GetRecoParam()->GetSigmaCutForTracking();
1263   Bool_t foundOneCluster, trackModified;
1264   AliMUONVCluster* cluster;
1265   AliMUONTrackParam *trackParam, *nextTrackParam, copyOfTrackParam, trackParamAtCluster, bestTrackParamAtCluster;
1266
1267   AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First();
1268   while (track) {
1269     trackModified = kFALSE;
1270
1271     trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
1272     while (trackParam) {
1273       foundOneCluster = kFALSE;
1274       bestChi2OfCluster = 2. * sigmaCut2; // 2 because 2 quantities in chi2
1275       chamberId = trackParam->GetClusterPtr()->GetChamberId();
1276       detElemId = trackParam->GetClusterPtr()->GetDetElemId();
1277
1278       // prepare nextTrackParam before adding new cluster because of the sorting
1279       nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam);
1280
1281       // recover track parameters from local fit and put them into a copy of trackParam
1282       copyOfTrackParam.SetZ(trackParam->GetZ());
1283       copyOfTrackParam.SetParameters(trackParam->GetSmoothParameters());
1284       copyOfTrackParam.SetCovariances(trackParam->GetSmoothCovariances());
1285
1286       // Create iterators to loop over clusters in current chamber
1287       TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId));
1288
1289       // look for one second candidate in the same chamber
1290       while ( ( cluster = static_cast<AliMUONVCluster*>(nextInCh()) ) ) {
1291
1292         // look for a cluster in another detection element
1293         if (cluster->GetDetElemId() == detElemId) continue;
1294
1295         // try to add the current cluster fast
1296         if (!TryOneClusterFast(copyOfTrackParam, cluster)) continue;
1297
1298         // try to add the current cluster accurately
1299         chi2OfCluster = TryOneCluster(copyOfTrackParam, cluster, trackParamAtCluster);
1300
1301         // if better chi2 then prepare to add this cluster to the track
1302         if (chi2OfCluster < bestChi2OfCluster) {
1303           bestChi2OfCluster = chi2OfCluster;
1304           bestTrackParamAtCluster = trackParamAtCluster;
1305           foundOneCluster = kTRUE;
1306         }
1307
1308       }
1309
1310       // add new cluster if any
1311       if (foundOneCluster) {
1312
1313         // Printout for debuging
1314         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1315           cout << "ComplementTracks: found one cluster in chamber(1..): " << chamberId+1 << endl;
1316           bestTrackParamAtCluster.GetClusterPtr()->Print();
1317           cout<<endl<<"Track parameters and covariances at cluster:"<<endl;
1318           bestTrackParamAtCluster.GetParameters().Print();
1319           bestTrackParamAtCluster.GetCovariances().Print();
1320         }
1321
1322         trackParam->SetRemovable(kTRUE);
1323         bestTrackParamAtCluster.SetRemovable(kTRUE);
1324         track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr()));
1325         trackModified = kTRUE;
1326       }
1327
1328       trackParam = nextTrackParam;
1329     }
1330
1331     // re-fit track parameters if needed
1332     if (trackModified) Fit(*track, kTRUE, kFALSE, kTRUE);
1333
1334     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1335   }
1336
1337 }
1338
1339   //__________________________________________________________________________
1340 void AliMUONTrackReconstructor::ImproveTrack(AliMUONTrack &track)
1341 {
1342   /// Improve the given track by removing clusters with local chi2 highter than the defined cut
1343   /// Recompute track parameters and covariances at the remaining clusters
1344   AliDebug(1,"Enter ImproveTrack");
1345
1346   Double_t localChi2, worstLocalChi2;
1347   AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster;
1348   Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForImprovement() *
1349                        GetRecoParam()->GetSigmaCutForImprovement();
1350
1351   while (!track.IsImproved()) {
1352
1353     // identify removable clusters
1354     track.TagRemovableClusters(GetRecoParam()->RequestedStationMask());
1355
1356     // Update track parameters and covariances
1357     track.UpdateCovTrackParamAtCluster();
1358
1359     // Compute local chi2 of each clusters
1360     track.ComputeLocalChi2(kTRUE);
1361
1362     // Look for the cluster to remove
1363     worstTrackParamAtCluster = NULL;
1364     worstLocalChi2 = 0.;
1365     trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1366     while (trackParamAtCluster) {
1367
1368       // Pick up cluster with the worst chi2
1369       localChi2 = trackParamAtCluster->GetLocalChi2();
1370       if (localChi2 > worstLocalChi2) {
1371         worstLocalChi2 = localChi2;
1372         worstTrackParamAtCluster = trackParamAtCluster;
1373       }
1374
1375     trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(trackParamAtCluster);
1376     }
1377
1378     // Check if worst cluster found
1379     if (!worstTrackParamAtCluster) {
1380       AliWarning("Bad local chi2 values?");
1381       break;
1382     }
1383
1384     // Check whether the worst chi2 is under requirement or not
1385     if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2
1386       track.SetImproved(kTRUE);
1387       break;
1388     }
1389
1390     // if the worst cluster is not removable then stop improvement
1391     if (!worstTrackParamAtCluster->IsRemovable()) break;
1392
1393     // Remove the worst cluster
1394     track.RemoveTrackParamAtCluster(worstTrackParamAtCluster);
1395
1396     // Re-fit the track:
1397     // Take into account the multiple scattering
1398     // Calculate the track parameter covariance matrix
1399     Fit(track, kTRUE, kFALSE, kTRUE);
1400
1401     // Printout for debuging
1402     if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
1403       cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl;
1404     }
1405
1406   }
1407
1408 }
1409
1410   //__________________________________________________________________________
1411 void AliMUONTrackReconstructor::FinalizeTrack(AliMUONTrack &track)
1412 {
1413   /// Recompute track parameters and covariances at each attached cluster
1414   /// from those at the first one, if not already done
1415   AliDebug(1,"Enter FinalizeTrack");
1416   if (!track.IsImproved()) track.UpdateCovTrackParamAtCluster();
1417 }
1418
1419   //__________________________________________________________________________
1420 Bool_t AliMUONTrackReconstructor::RefitTrack(AliMUONTrack &track, Bool_t enableImprovement)
1421 {
1422   /// re-fit the given track
1423
1424   // check validity of the track
1425   if (track.GetNClusters() < 3) {
1426     AliWarning("the track does not contain enough clusters --> unable to refit");
1427     return kFALSE;
1428   }
1429
1430   // reset the seed (i.e. parameters at first cluster) before fitting
1431   AliMUONTrackParam* firstTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
1432   if (firstTrackParam->GetInverseBendingMomentum() == 0.) {
1433     AliWarning("track parameters at first chamber are not initialized --> unable to refit");
1434     return kFALSE;
1435   }
1436
1437   // compute track parameters at each cluster from parameters at the first one
1438   // necessary to compute multiple scattering effect during refitting
1439   track.UpdateTrackParamAtCluster();
1440
1441   // Re-fit the track:
1442   // Take into account the multiple scattering
1443   // Calculate the track parameter covariance matrix
1444   Fit(track, kTRUE, kFALSE, kTRUE);
1445
1446   // Improve the reconstructed tracks if required
1447   track.SetImproved(kFALSE);
1448   if (enableImprovement && GetRecoParam()->ImproveTracks()) ImproveTrack(track);
1449
1450   // Fill AliMUONTrack data members
1451   FinalizeTrack(track);
1452
1453   return kTRUE;
1454
1455 }
1456