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