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