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 "AliMUONData.h"
  33 #include "AliMUONConstants.h"
  34 #include "AliMUONHitForRec.h"
  35 #include "AliMUONObjectPair.h"
  36 #include "AliMUONRawCluster.h"
  37 #include "AliMUONTrackK.h"
  38
  39 #include "AliLog.h"
  40
  41 #include <Riostream.h>
  42
  43 /// \cond CLASSIMP
  44 ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context
  45 ClassImp(AliMUONConstants)
  46 /// \endcond
  47
  48 //__________________________________________________________________________
  49 AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod)
  50   : AliMUONVTrackReconstructor(data),
  51     fTrackMethod(2), //tracking method (2-Kalman 3-Combination-Kalman/Clustering)
  52     fMuons(0)
  53 {
  54   /// Constructor for class AliMUONTrackReconstructorK
  55
  56   if (strstr(TrackMethod,"Kalman")) {
  57     cout << " *** Tracking with the Kalman filter *** " << endl;
  58     fTrackMethod = 2;
  59   } else if (strstr(TrackMethod,"Combi")) {
  60     cout << " *** Combined cluster / track finder ***" << endl;
  61     fTrackMethod = 3;
  62   } else AliFatal(Form("Tracking method %s not available",TrackMethod));
  63
  64   // Memory allocation for the TClonesArray of reconstructed tracks
  65   fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
  66 }
  67
  68   //__________________________________________________________________________
  69 AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK(void)
  70 {
  71   /// Destructor for class AliMUONTrackReconstructorK
  72   delete fRecTracksPtr;
  73 }
  74
  75   //__________________________________________________________________________
  76 void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
  77 {
  78   /// To add to the list of hits for reconstruction all the raw clusters
  79   TTree *treeR;
  80   AliMUONHitForRec *hitForRec;
  81   AliMUONRawCluster *clus;
  82   Int_t iclus, nclus;
  83   TClonesArray *rawclusters;
  84   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
  85
  86   treeR = fMUONData->TreeR();
  87   if (!treeR) {
  88     AliError("TreeR must be loaded");
  89     exit(0);
  90   }
  91
  92   if (fTrackMethod != 3) { //AZ
  93     fMUONData->SetTreeAddress("RC"); //AZ
  94     fMUONData->GetRawClusters(); // only one entry
  95   }
  96
  97   // Loop over tracking chambers
  98   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
  99     rawclusters =fMUONData->RawClusters(ch);
 100     nclus = (Int_t) (rawclusters->GetEntries());
 101     // Loop over (cathode correlated) raw clusters
 102     for (iclus = 0; iclus < nclus; iclus++) {
 103       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
 104       // new AliMUONHitForRec from raw cluster
 105       // and increment number of AliMUONHitForRec's (total and in chamber)
 106       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
 107       fNHitsForRec++;
 108       // more information into HitForRec
 109       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
 110       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
 111       //  original raw cluster
 112       hitForRec->SetChamberNumber(ch);
 113       hitForRec->SetHitNumber(iclus);
 114       // Z coordinate of the raw cluster (cm)
 115       hitForRec->SetZ(clus->GetZ(0));
 116       if (AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 3) {
 117         cout << "Chamber " << ch <<" raw cluster  " << iclus << " : " << endl;
 118         clus->Print("full");
 119         cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 120         hitForRec->Print("full");
 121       }
 122     } // end of cluster loop
 123   } // end of chamber loop
 124   SortHitsForRecWithIncreasingChamber();
 125
 126   AliDebug(1,"End of AddHitsForRecFromRawClusters");
 127     if (AliLog::GetGlobalDebugLevel() > 0) {
 128       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
 129       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
 130         AliDebug(1, Form("Chamber(0...): %d",ch));
 131         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
 132         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
 133         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
 134              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
 135              hit++) {
 136           AliDebug(1, Form("HitForRec index(0...): %d",hit));
 137           ((*fHitsForRecPtr)[hit])->Dump();
 138       }
 139     }
 140   }
 141
 142   return;
 143 }
 144
 145   //__________________________________________________________________________
 146 void AliMUONTrackReconstructorK::MakeTracks(void)
 147 {
 148   /// To make the tracks from the list of segments and points in all stations
 149   AliDebug(1,"Enter MakeTracks");
 150   // The order may be important for the following Reset's
 151   //AZ ResetTracks();
 152   MakeTrackCandidates();
 153   if (fRecTracksPtr->GetEntriesFast() == 0) return;
 154   // Follow tracks in stations(1..) 3, 2 and 1
 155   FollowTracks();
 156   // Remove double tracks
 157   RemoveDoubleTracks();
 158   // Fill AliMUONTrack data members
 159   FillMUONTrack();
 160 }
 161
 162   //__________________________________________________________________________
 163 void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
 164 {
 165   /// To make initial tracks for Kalman filter from segments in stations(1..)  4 and 5
 166   Int_t istat, iseg;
 167   TClonesArray *segments;
 168   AliMUONObjectPair *segment;
 169   AliMUONTrackK *trackK;
 170
 171   AliDebug(1,"Enter MakeTrackCandidatesK");
 172
 173   AliMUONTrackK a(this, fHitsForRecPtr);
 174   // Loop over stations(1...) 5 and 4
 175   for (istat=4; istat>=3; istat--) {
 176     // Make segments in the station
 177     segments = MakeSegmentsInStation(istat);
 178     // Loop over segments in the station
 179     for (iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
 180       // Transform segments to tracks
 181       segment = (AliMUONObjectPair*) ((*segments)[iseg]);
 182       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
 183     } // for (iseg=0;...)
 184   } // for (istat=4;...)
 185   return;
 186 }
 187
 188 //__________________________________________________________________________
 189 void AliMUONTrackReconstructorK::FollowTracks(void)
 190 {
 191   /// Follow tracks using Kalman filter
 192   Bool_t ok;
 193   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
 194   AliMUONTrackK *trackK;
 195   AliMUONHitForRec *hit;
 196   AliMUONRawCluster *clus;
 197   TClonesArray *rawclusters;
 198   clus = 0; rawclusters = 0;
 199
 200   Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
 201   Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
 202   Double_t zDipole1 = simpleBPosition + 0.5 * simpleBLength;
 203   Double_t zDipole2 = zDipole1 - simpleBLength;
 204
 205   // Print hits
 206   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
 207
 208   if (trackK->DebugLevel() > 0) {
 209     for (Int_t i1=0; i1<fNHitsForRec; i1++) {
 210       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
 211       printf(" Hit # %d %10.4f %10.4f %10.4f",
 212              hit->GetChamberNumber(), hit->GetBendingCoor(),
 213              hit->GetNonBendingCoor(), hit->GetZ());
 214
 215       // from raw clusters
 216       rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
 217       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
 218                                                            GetHitNumber());
 219       printf(" %d", clus->GetTrack(1));
 220       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
 221       else printf("\n");
 222
 223     }
 224   } // if (trackK->DebugLevel() > 0)
 225
 226   icand = -1;
 227   Int_t nSeeds;
 228   nSeeds = fNRecTracks; // starting number of seeds
 229   // Loop over track candidates
 230   while (icand < fNRecTracks-1) {
 231     icand ++;
 232     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
 233     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
 234     if (trackK->GetRecover() < 0) continue; // failed track
 235
 236     // Discard candidate which will produce the double track
 237     /*
 238     if (icand > 0) {
 239       ok = CheckCandidate(icand,nSeeds);
 240       if (!ok) {
 241         trackK->SetRecover(-1); // mark candidate to be removed
 242         continue;
 243       }
 244     }
 245     */
 246
 247     ok = kTRUE;
 248     if (trackK->GetRecover() == 0)
 249       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
 250     else
 251       hit = trackK->GetHitLastOk(); // hit where track stopped
 252
 253     if (hit) ichamBeg = hit->GetChamberNumber();
 254     ichamEnd = 0;
 255     // Check propagation direction
 256     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
 257     else if (trackK->GetTrackDir() < 0) {
 258       ichamEnd = 9; // forward propagation
 259       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
 260       if (ok) {
 261         ichamBeg = ichamEnd;
 262         ichamEnd = 6; // backward propagation
 263         // Change weight matrix and zero fChi2 for backpropagation
 264         trackK->StartBack();
 265         trackK->SetTrackDir(1);
 266         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
 267         ichamBeg = ichamEnd;
 268         ichamEnd = 0;
 269       }
 270     } else {
 271       if (trackK->GetBPFlag()) {
 272         // backpropagation
 273         ichamEnd = 6; // backward propagation
 274         // Change weight matrix and zero fChi2 for backpropagation
 275         trackK->StartBack();
 276         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
 277         ichamBeg = ichamEnd;
 278         ichamEnd = 0;
 279       }
 280     }
 281
 282     if (ok) {
 283       trackK->SetTrackDir(1);
 284       trackK->SetBPFlag(kFALSE);
 285       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
 286     }
 287     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
 288
 289     // Apply smoother
 290     if (trackK->GetRecover() >= 0) {
 291       ok = trackK->Smooth();
 292       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
 293     }
 294
 295     // Majority 3 of 4 in first 2 stations
 296     if (!ok) continue;
 297     chamBits = 0;
 298     Double_t chi2max = 0;
 299     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
 300       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
 301       chamBits |= BIT(hit->GetChamberNumber());
 302       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
 303     }
 304     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
 305       //trackK->Recover();
 306       trackK->SetRecover(-1); //mark candidate to be removed
 307       continue;
 308     }
 309     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
 310   } // while
 311
 312   for (Int_t i=0; i<fNRecTracks; i++) {
 313     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
 314     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
 315   }
 316
 317   // Compress TClonesArray
 318   fRecTracksPtr->Compress();
 319   fNRecTracks = fRecTracksPtr->GetEntriesFast();
 320   return;
 321 }
 322
 323 //__________________________________________________________________________
 324 Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const
 325 {
 326   /// Discards track candidate if it will produce the double track (having
 327   /// the same seed segment hits as hits of a good track found before)
 328   AliMUONTrackK *track1, *track2;
 329   AliMUONHitForRec *hit1, *hit2, *hit;
 330   AliMUONObjectPair *segment1, *segment2;
 331
 332   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
 333   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
 334   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
 335
 336   for (Int_t i=0; i<icand; i++) {
 337     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
 338     //if (track2->GetRecover() < 0) continue;
 339     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
 340     segment1 = track1->GetStartSegment();
 341     segment2 = track2->GetStartSegment();
 342     if (segment1->First()  == segment2->First() &&
 343         segment1->Second() == segment2->Second()) {
 344       return kFALSE;
 345     } else {
 346       Int_t nSame = 0;
 347       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
 348         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
 349         if (hit == hit1 || hit == hit2) {
 350           nSame++;
 351           if (nSame == 2) return kFALSE;
 352         }
 353       } // for (Int_t j=0;
 354     }
 355   } // for (Int_t i=0;
 356   return kTRUE;
 357 }
 358
 359   //__________________________________________________________________________
 360 void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
 361 {
 362   /// Removes double tracks (sharing more than half of their hits).
 363   /// Keeps the track with higher quality
 364   AliMUONTrackK *track1, *track2, *trackToKill;
 365
 366   // Sort tracks according to their quality
 367   fRecTracksPtr->Sort();
 368
 369   // Loop over first track of the pair
 370   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
 371   Int_t debug = track1->DebugLevel();
 372   while (track1) {
 373     // Loop over second track of the pair
 374     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
 375     while (track2) {
 376       // Check whether or not to keep track2
 377       if (!track2->KeepTrack(track1)) {
 378         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
 379           " " << track2->GetTrackQuality() << endl;
 380         trackToKill = track2;
 381         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
 382         trackToKill->Kill();
 383         fRecTracksPtr->Compress();
 384       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
 385     } // track2
 386     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
 387   } // track1
 388
 389   fNRecTracks = fRecTracksPtr->GetEntriesFast();
 390   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
 391 }
 392
 393   //__________________________________________________________________________
 394 void AliMUONTrackReconstructorK::FillMUONTrack()
 395 {
 396   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
 397   AliMUONTrackK *track;
 398   track = (AliMUONTrackK*) fRecTracksPtr->First();
 399   while (track) {
 400     track->FillMUONTrack();
 401     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
 402   }
 403   return;
 404 }
 405
 406   //__________________________________________________________________________
 407 void AliMUONTrackReconstructorK::EventDump(void)
 408 {
 409   /// Dump reconstructed event (track parameters at vertex and at first hit),
 410   /// and the particle parameters
 411   AliDebug(1,"****** enter EventDump ******");
 412   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
 413
 414   AliWarning("AliMUONTrackReconstructorK::EventDump: not implemented");
 415
 416   // informations about generated particles NO !!!!!!!!
 417
 418 //    for (Int_t iPart = 0; iPart < np; iPart++) {
 419 //      p = gAlice->Particle(iPart);
 420 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
 421 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
 422 //    }
 423   return;
 424 }
 425
 426