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