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