MUON/AliMUONTrackReconstructorK.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 /// \class AliMUONTrackReconstructorK
  19 /////////////////////////////////////
  20 ///
  21 /// MUON track reconstructor using the kalman method
  22 ///
  23 /// This class contains as data:
  24 /// - the parameters for the track reconstruction
  25 ///
  26 /// It contains as methods, among others:
  27 /// - MakeTracks to build the tracks
  28 ///
  29 ////////////////////////////////////
  30
  31 #include "AliMUONTrackReconstructorK.h"
  32 #include "AliMUONConstants.h"
  33 #include "AliMUONHitForRec.h"
  34 #include "AliMUONObjectPair.h"
  35 #include "AliMUONRawCluster.h"
  36 #include "AliMUONTrackK.h"
  37
  38 #include "AliLog.h"
  39
  40 #include <Riostream.h>
  41
  42 /// \cond CLASSIMP
  43 ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
  44 ClassImp(AliMUONConstants)
  45 /// \endcond
  46
  47 //__________________________________________________________________________
  48 AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(const Option_t* TrackMethod)
  49   : AliMUONVTrackReconstructor(),
  50     fTrackMethod(2) //tracking method (2-Kalman 3-Combination-Kalman/Clustering)
  51 {
  52   /// Constructor for class AliMUONTrackReconstructorK
  53
  54   if (strstr(TrackMethod,"Kalman")) {
  55     AliInfo(" *** Tracking with the Kalman filter *** ");
  56     fTrackMethod = 2;
  57   } else if (strstr(TrackMethod,"Combi")) {
  58     AliInfo(" *** Combined cluster / track finder ***");
  59     fTrackMethod = 3;
  60   } else AliFatal(Form("Tracking method %s not available",TrackMethod));
  61
  62   // Memory allocation for the TClonesArray of reconstructed tracks
  63   fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
  64 }
  65
  66 //__________________________________________________________________________
  67 AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK(void)
  68 {
  69   /// Destructor for class AliMUONTrackReconstructorK
  70   delete fRecTracksPtr;
  71 }
  72
  73 //__________________________________________________________________________
  74 void AliMUONTrackReconstructorK::MakeTracks(void)
  75 {
  76   /// To make the tracks from the list of segments and points in all stations
  77   AliDebug(1,"Enter MakeTracks");
  78   // The order may be important for the following Reset's
  79   //AZ ResetTracks();
  80   MakeTrackCandidates();
  81   if (fRecTracksPtr->GetEntriesFast() == 0) return;
  82   // Follow tracks in stations(1..) 3, 2 and 1
  83   FollowTracks();
  84   // Remove double tracks
  85   RemoveDoubleTracks();
  86   // Print out some track info (if necessary)
  87   EventDump();
  88   // Fill AliMUONTrack data members
  89   FillMUONTrack();
  90 }
  91
  92   //__________________________________________________________________________
  93 void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
  94 {
  95   /// To make initial tracks for Kalman filter from segments in stations(1..)  4 and 5
  96   Int_t istat, iseg;
  97   TClonesArray *segments;
  98   AliMUONObjectPair *segment;
  99   AliMUONTrackK *trackK;
 100
 101   AliDebug(1,"Enter MakeTrackCandidatesK");
 102
 103   AliMUONTrackK a(this, fHitsForRecPtr);
 104   // Loop over stations(1...) 5 and 4
 105   for (istat=4; istat>=3; istat--)
 106   {
 107     // Make segments in the station
 108     segments = MakeSegmentsInStation(istat);
 109     // Loop over segments in the station
 110     for (iseg=0; iseg<segments->GetEntriesFast(); iseg++)
 111     {
 112       // Transform segments to tracks
 113       segment = (AliMUONObjectPair*) ((*segments)[iseg]);
 114       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(*segment);
 115     } // for (iseg=0;...)
 116     delete segments;
 117   } // for (istat=4;...)
 118 }
 119
 120 //__________________________________________________________________________
 121 void AliMUONTrackReconstructorK::FollowTracks(void)
 122 {
 123   /// Follow tracks using Kalman filter
 124   Bool_t ok;
 125   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
 126   AliMUONTrackK *trackK;
 127   AliMUONHitForRec *hit;
 128   AliMUONRawCluster *clus;
 129   TClonesArray *rawclusters;
 130   clus = 0; rawclusters = 0;
 131
 132   Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
 133   Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
 134   Double_t zDipole1 = simpleBPosition + 0.5 * simpleBLength;
 135   Double_t zDipole2 = zDipole1 - simpleBLength;
 136
 137   // Print hits
 138   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
 139
 140   icand = -1;
 141   Int_t nSeeds;
 142   nSeeds = fNRecTracks; // starting number of seeds
 143   // Loop over track candidates
 144   while (icand < fNRecTracks-1) {
 145     icand ++;
 146     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
 147     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
 148     if (trackK->GetRecover() < 0) continue; // failed track
 149
 150     // Discard candidate which will produce the double track
 151     /*
 152     if (icand > 0) {
 153       ok = CheckCandidate(icand,nSeeds);
 154       if (!ok) {
 155         trackK->SetRecover(-1); // mark candidate to be removed
 156         continue;
 157       }
 158     }
 159     */
 160
 161     ok = kTRUE;
 162     if (trackK->GetRecover() == 0)
 163       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
 164     else
 165       hit = trackK->GetHitLastOk(); // hit where track stopped
 166
 167     if (hit) ichamBeg = hit->GetChamberNumber();
 168     if (ichamBeg == 8 && trackK->GetStation0() == 4) ichamBeg = 10;
 169
 170     ichamEnd = 0;
 171     // Check propagation direction
 172     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
 173     else if (trackK->GetTrackDir() < 0) {
 174       ichamEnd = 9; // forward propagation
 175       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
 176       if (ok) {
 177         ichamBeg = ichamEnd + 1;
 178         ichamEnd = 6; // backward propagation
 179         // Change weight matrix and zero fChi2 for backpropagation
 180         trackK->StartBack();
 181         trackK->SetTrackDir(1);
 182         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
 183         ichamBeg = ichamEnd;
 184         ichamEnd = 0;
 185       }
 186     } else {
 187       if (trackK->GetBPFlag()) {
 188         // backpropagation
 189         if (ichamBeg == 10) ichamEnd = 8;
 190         else {
 191           if (ichamBeg == 9) ichamBeg++;
 192           ichamEnd = 6;
 193           // Change weight matrix and zero fChi2 for backpropagation
 194           trackK->StartBack();
 195         }
 196         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
 197         ichamBeg = ichamEnd;
 198         ichamEnd = 0;
 199       }
 200     }
 201
 202     if (ok) {
 203       trackK->SetTrackDir(1);
 204       trackK->SetBPFlag(kFALSE);
 205       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
 206     }
 207     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
 208
 209     // Apply smoother
 210     if (trackK->GetRecover() >= 0) {
 211       ok = trackK->Smooth();
 212       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
 213     }
 214
 215     // Majority 3 of 4 in first 2 stations
 216     if (!ok) continue;
 217     chamBits = 0;
 218     Double_t chi2max = 0;
 219     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
 220       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
 221       chamBits |= BIT(hit->GetChamberNumber());
 222       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
 223     }
 224     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
 225       //trackK->Recover();
 226       trackK->SetRecover(-1); //mark candidate to be removed
 227       continue;
 228     }
 229     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
 230   } // while
 231
 232   for (Int_t i=0; i<fNRecTracks; i++) {
 233     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
 234     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
 235   }
 236
 237   // Compress TClonesArray
 238   fRecTracksPtr->Compress();
 239   fNRecTracks = fRecTracksPtr->GetEntriesFast();
 240   return;
 241 }
 242
 243 //__________________________________________________________________________
 244 Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const
 245 {
 246   /// Discards track candidate if it will produce the double track (having
 247   /// the same seed segment hits as hits of a good track found before)
 248   AliMUONTrackK *track1, *track2;
 249   AliMUONHitForRec *hit1, *hit2, *hit;
 250   AliMUONObjectPair *segment1, *segment2;
 251
 252   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
 253   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
 254   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
 255
 256   for (Int_t i=0; i<icand; i++) {
 257     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
 258     //if (track2->GetRecover() < 0) continue;
 259     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
 260     segment1 = track1->GetStartSegment();
 261     segment2 = track2->GetStartSegment();
 262     if (segment1->First()  == segment2->First() &&
 263         segment1->Second() == segment2->Second()) {
 264       return kFALSE;
 265     } else {
 266       Int_t nSame = 0;
 267       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
 268         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
 269         if (hit == hit1 || hit == hit2) {
 270           nSame++;
 271           if (nSame == 2) return kFALSE;
 272         }
 273       } // for (Int_t j=0;
 274     }
 275   } // for (Int_t i=0;
 276   return kTRUE;
 277 }
 278
 279   //__________________________________________________________________________
 280 void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
 281 {
 282   /// Removes double tracks (sharing more than half of their hits).
 283   /// Keeps the track with higher quality
 284   AliMUONTrackK *track1, *track2, *trackToKill;
 285
 286   // Sort tracks according to their quality
 287   fRecTracksPtr->Sort();
 288
 289   // Loop over first track of the pair
 290   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
 291   Int_t debug = track1->DebugLevel();
 292   while (track1) {
 293     // Loop over second track of the pair
 294     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
 295     while (track2) {
 296       // Check whether or not to keep track2
 297       if (!track2->KeepTrack(track1)) {
 298         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
 299           " " << track2->GetTrackQuality() << endl;
 300         trackToKill = track2;
 301         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
 302         trackToKill->Kill();
 303         fRecTracksPtr->Compress();
 304       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
 305     } // track2
 306     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
 307   } // track1
 308
 309   fNRecTracks = fRecTracksPtr->GetEntriesFast();
 310   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
 311 }
 312
 313   //__________________________________________________________________________
 314 void AliMUONTrackReconstructorK::FillMUONTrack()
 315 {
 316   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
 317   AliMUONTrackK *track;
 318   track = (AliMUONTrackK*) fRecTracksPtr->First();
 319   while (track) {
 320     track->FillMUONTrack();
 321     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
 322   }
 323   return;
 324 }
 325
 326
 327   //__________________________________________________________________________
 328 void AliMUONTrackReconstructorK::EventDump(void)
 329 {
 330   /// Dump reconstructed event (track parameters at vertex and at first hit),
 331   /// and the particle parameters
 332   AliDebug(1,"****** enter EventDump ******");
 333   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
 334
 335   Int_t debug =   ((AliMUONTrackK*) fRecTracksPtr->First())->DebugLevel();
 336   if (debug < 0) return;
 337
 338   for (Int_t i = 0; i < fNRecTracks; ++i) {
 339     AliMUONTrackK *track = (AliMUONTrackK*) fRecTracksPtr->UncheckedAt(i);
 340     track->Print("");
 341   }
 342 }