MUON/AliMUONEventReconstructor.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 /*
  17 $Log$
  18 Revision 1.1.2.5  2000/06/16 07:00:26  gosset
  19 To remove problem in running RuleChecker
  20
  21 Revision 1.1.2.4  2000/06/12 08:00:07  morsch
  22 Dummy streamer to solve CINT compilation problem (to be investigated !)
  23
  24 Revision 1.1.2.3  2000/06/09 20:59:57  morsch
  25 Make includes consistent with new file structure.
  26
  27 Revision 1.1.2.2  2000/06/09 12:58:05  gosset
  28 Removed comment beginnings in Log sections of .cxx files
  29 Suppressed most violations of coding rules
  30
  31 Revision 1.1.2.1  2000/06/07 14:44:53  gosset
  32 Addition of files for track reconstruction in C++
  33 */
  34
  35 //__________________________________________________________________________
  36 //
  37 // MUON event reconstructor in ALICE
  38 //
  39 // This class contains as data:
  40 // * the parameters for the event reconstruction
  41 // * a pointer to the array of hits to be reconstructed (the event)
  42 // * a pointer to the array of segments made with these hits inside each station
  43 // * a pointer to the array of reconstructed tracks
  44 //
  45 // It contains as methods, among others:
  46 // * MakeEventToBeReconstructed to build the array of hits to be reconstructed
  47 // * MakeSegments to build the segments
  48 // * MakeTracks to build the tracks
  49 //__________________________________________________________________________
  50
  51 #include <iostream.h>
  52
  53 #include <TRandom.h>
  54 #include <TFile.h>
  55
  56 #include "AliCallf77.h"
  57 #include "AliMUONEventReconstructor.h"
  58 #include "AliMUON.h"
  59 #include "AliMUONHitForRec.h"
  60 #include "AliMUONSegment.h"
  61 #include "AliMUONHit.h"
  62 #include "AliMUONRawCluster.h"
  63 #include "AliMUONTrack.h"
  64 #include "AliMUONChamber.h"
  65 #include "AliMUONTrackHit.h"
  66 #include "AliRun.h"
  67
  68 #ifndef WIN32
  69 # define initfield initfield_
  70 # define reco_gufld reco_gufld_
  71 #else
  72 # define initfield INITFIELD
  73 # define reco_gufld RECO_GUFLD
  74 #endif
  75
  76 extern "C"
  77 {
  78 void type_of_call initfield();
  79 void type_of_call reco_gufld(Double_t *Coor, Double_t *Field);
  80 }
  81
  82 //************* Defaults parameters for reconstruction
  83 static const Double_t DefaultMinBendingMomentum = 3.0;
  84 static const Double_t DefaultMaxSigma2Distance = 16.0;
  85 static const Double_t DefaultBendingResolution = 0.01;
  86 static const Double_t DefaultNonBendingResolution = 0.144;
  87 static const Double_t DefaultChamberThicknessInX0 = 0.03;
  88 // Simple magnetic field:
  89 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
  90 // Length and Position from reco_muon.F, with opposite sign:
  91 // Length = ZMAGEND-ZCOIL
  92 // Position = (ZMAGEND+ZCOIL)/2
  93 // to be ajusted differently from real magnetic field ????
  94 static const Double_t DefaultSimpleBValue = 7.0;
  95 static const Double_t DefaultSimpleBLength = 428.0;
  96 static const Double_t DefaultSimpleBPosition = 1019.0;
  97 static const Int_t DefaultRecGeantHits = 0;
  98 static const Double_t DefaultEfficiency = 0.95;
  99
 100 static const Int_t DefaultPrintLevel = 0;
 101
 102 ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context
 103
 104   //__________________________________________________________________________
 105 AliMUONEventReconstructor::AliMUONEventReconstructor(void)
 106 {
 107   // Constructor for class AliMUONEventReconstructor
 108   SetReconstructionParametersToDefaults();
 109   // Memory allocation for the TClonesArray of hits for reconstruction
 110   // Is 10000 the right size ????
 111   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
 112   fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ????
 113   // Memory allocation for the TClonesArray's of segments in stations
 114   // Is 2000 the right size ????
 115   for (Int_t st = 0; st < MAX_MUON_TRACKING_STATIONS; st++) {
 116     fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
 117     fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
 118   }
 119   // Memory allocation for the TClonesArray of reconstructed tracks
 120   // Is 10 the right size ????
 121   fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
 122   fNRecTracks = 0; // really needed or GetEntriesFast sufficient ????
 123
 124   // Initialize magnetic field
 125   // using Fortran subroutine INITFIELD in "reco_muon.F".
 126   // Should rather use AliMagF ???? and remove prototyping ...
 127   initfield();
 128   // Impression de quelques valeurs
 129   Double_t coor[3], field[3];
 130   coor[0] = 50.0;
 131   coor[1] = 50.0;
 132   coor[2] = 950.0;
 133   reco_gufld(coor, field);
 134   cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl;
 135   cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl;
 136   coor[2] = -950.0;
 137   reco_gufld(coor, field);
 138   cout << "coor: " << coor[0] << ", " << coor[1] << ", " << coor[2] << endl;
 139   cout << "field: " << field[0] << ", " << field[1] << ", " << field[2] << endl;
 140   coor[2] = -950.0;
 141
 142   if (fPrintLevel >= 0) {
 143     cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();}
 144   return;
 145 }
 146
 147   //__________________________________________________________________________
 148 AliMUONEventReconstructor::~AliMUONEventReconstructor(void)
 149 {
 150   // Destructor for class AliMUONEventReconstructor
 151   delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
 152   for (Int_t st = 0; st < MAX_MUON_TRACKING_STATIONS; st++)
 153     delete fSegmentsPtr[st]; // Correct destruction of everything ????
 154   return;
 155 }
 156
 157   //__________________________________________________________________________
 158 void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
 159 {
 160   // Set reconstruction parameters to default values
 161   // Would be much more convenient with a structure (or class) ????
 162   fMinBendingMomentum = DefaultMinBendingMomentum;
 163   fMaxSigma2Distance = DefaultMaxSigma2Distance;
 164
 165   AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
 166   // ******** Parameters for making HitsForRec
 167   // minimum radius,
 168   // like in TRACKF_STAT:
 169   // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
 170   // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
 171   for (Int_t ch = 0; ch < MAX_MUON_TRACKING_CHAMBERS; ch++) {
 172     if (ch < 4) fRMin[ch] = TMath::Abs((&(MUON->Chamber(ch)))->Z()) *
 173                   2.0 * TMath::Pi() / 180.0;
 174     else fRMin[ch] = 30.0;
 175   }
 176   // maximum radius
 177   // like in TRACKF_STAT (10 degrees ????)
 178   fRMax[0] = fRMax[1] = 91.5;
 179   fRMax[2] = fRMax[3] = 122.5;
 180   fRMax[4] = fRMax[5] = 158.3;
 181   fRMax[6] = fRMax[7] = 260.0;
 182   fRMax[8] = fRMax[9] = 260.0;
 183
 184   // ******** Parameters for making segments
 185   // should be parametrized ????
 186   // according to interval between chambers in a station ????
 187   // Maximum distance in non bending plane
 188   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
 189   // SIGCUT*DYMAX(IZ)
 190   for (Int_t st = 0; st < MAX_MUON_TRACKING_STATIONS; st++)
 191     fSegmentMaxDistNonBending[st] = 5. * 0.22;
 192   // Maximum distance in bending plane
 193   // values from TRACKF_STAT corresponding to (J psi 20cm)
 194   fSegmentMaxDistBending[0] = 1.5;
 195   fSegmentMaxDistBending[1] = 1.5;
 196   fSegmentMaxDistBending[2] = 3.0;
 197   fSegmentMaxDistBending[3] = 6.0;
 198   fSegmentMaxDistBending[4] = 6.0;
 199
 200   fBendingResolution = DefaultBendingResolution;
 201   fNonBendingResolution = DefaultNonBendingResolution;
 202   fChamberThicknessInX0 = DefaultChamberThicknessInX0;
 203   fSimpleBValue = DefaultSimpleBValue;
 204   fSimpleBLength = DefaultSimpleBLength;
 205   fSimpleBPosition = DefaultSimpleBPosition;
 206   fRecGeantHits = DefaultRecGeantHits;
 207   fEfficiency = DefaultEfficiency;
 208   fPrintLevel = DefaultPrintLevel;
 209   return;
 210 }
 211
 212 //__________________________________________________________________________
 213 Double_t AliMUONEventReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum)
 214 {
 215   // Returns impact parameter at vertex in bending plane (cm),
 216   // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
 217   // using simple values for dipole magnetic field.
 218   // The sign is the sign of the charge.
 219   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 220           BendingMomentum);
 221 }
 222
 223 //__________________________________________________________________________
 224 Double_t AliMUONEventReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam)
 225 {
 226   // Returns signed bending momentum in bending plane (GeV/c),
 227   // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
 228   // using simple values for dipole magnetic field.
 229   // The sign is the sign of the charge.
 230   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
 231           ImpactParam);
 232 }
 233
 234 //__________________________________________________________________________
 235 void AliMUONEventReconstructor::SetBkgGeantFile(Text_t *BkgGeantFileName)
 236 {
 237   // Set background file ... for GEANT hits
 238   // Must be called after having loaded the firts signal event
 239   if (fPrintLevel >= 0) {
 240     cout << "Enter SetBkgGeantFile with BkgGeantFileName ``"
 241          << BkgGeantFileName << "''" << endl;}
 242   if (strlen(BkgGeantFileName)) {
 243     // BkgGeantFileName not empty: try to open the file
 244     if (fPrintLevel >= 2) {cout << "Before File(Bkg)" << endl; gDirectory->Dump();}
 245     fBkgGeantFile = new TFile(BkgGeantFileName);
 246     if (fPrintLevel >= 2) {cout << "After File(Bkg)" << endl; gDirectory->Dump();}
 247     if (fBkgGeantFile-> IsOpen()) {
 248       if (fPrintLevel >= 0) {
 249         cout << "Background for GEANT hits in file: ``" << BkgGeantFileName
 250              << "'' successfully opened" << endl;}
 251     }
 252     else {
 253       cout << "Background for GEANT hits in file: " << BkgGeantFileName << endl;
 254       cout << "NOT FOUND: EXIT" << endl;
 255       exit(0); // right instruction for exit ????
 256     }
 257     // Arrays for "particles" and "hits"
 258     fBkgGeantParticles = new TClonesArray("TParticle", 200);
 259     fBkgGeantHits = new TClonesArray("AliMUONHit", 2000);
 260     // Event number to -1 for initialization
 261     fBkgGeantEventNumber = -1;
 262     // Back to the signal file:
 263     // first signal event must have been loaded previously,
 264     // otherwise, Segmentation violation at the next instruction
 265     // How is it possible to do smething better ????
 266     ((gAlice->TreeK())->GetCurrentFile())->cd();
 267     if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
 268   }
 269   return;
 270 }
 271
 272 //__________________________________________________________________________
 273 void AliMUONEventReconstructor::NextBkgGeantEvent(void)
 274 {
 275   // Get next event in background file for GEANT hits
 276   // Goes back to event number 0 when end of file is reached
 277   char treeName[20];
 278   TBranch *branch;
 279   if (fPrintLevel >= 0) {
 280     cout << "Enter NextBkgGeantEvent" << endl;}
 281   // Clean previous event
 282   if(fBkgGeantTK) delete fBkgGeantTK;
 283   fBkgGeantTK = NULL;
 284   if(fBkgGeantParticles) fBkgGeantParticles->Clear();
 285   if(fBkgGeantTH) delete fBkgGeantTH;
 286   fBkgGeantTH = NULL;
 287   if(fBkgGeantHits) fBkgGeantHits->Clear();
 288   // Increment event number
 289   fBkgGeantEventNumber++;
 290   // Get access to Particles and Hits for event from background file
 291   if (fPrintLevel >= 2) {cout << "Before cd(Bkg)" << endl; gDirectory->Dump();}
 292   fBkgGeantFile->cd();
 293   if (fPrintLevel >= 2) {cout << "After cd(Bkg)" << endl; gDirectory->Dump();}
 294   // Particles: TreeK for event and branch "Particles"
 295   sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
 296   fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
 297   if (!fBkgGeantTK) {
 298     if (fPrintLevel >= 0) {
 299       cout << "Cannot find Kine Tree for background event: " <<
 300         fBkgGeantEventNumber << endl;
 301       cout << "Goes back to event 0" << endl;
 302     }
 303     fBkgGeantEventNumber = 0;
 304     sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
 305     fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
 306     if (!fBkgGeantTK) {
 307       cout << "ERROR: cannot find Kine Tree for background event: " <<
 308         fBkgGeantEventNumber << endl;
 309       exit(0);
 310     }
 311   }
 312   if (fBkgGeantTK)
 313     fBkgGeantTK->SetBranchAddress("Particles", &fBkgGeantParticles);
 314   fBkgGeantTK->GetEvent(0); // why event 0 ???? necessary ????
 315   // Hits: TreeH for event and branch "MUON"
 316   sprintf(treeName, "TreeH%d", fBkgGeantEventNumber);
 317   fBkgGeantTH = (TTree*)gDirectory->Get(treeName);
 318   if (!fBkgGeantTH) {
 319     cout << "ERROR: cannot find Hits Tree for background event: " <<
 320       fBkgGeantEventNumber << endl;
 321       exit(0);
 322   }
 323   if (fBkgGeantTH && fBkgGeantHits) {
 324     branch = fBkgGeantTH->GetBranch("MUON");
 325     if (branch) branch->SetAddress(&fBkgGeantHits);
 326   }
 327   fBkgGeantTH->GetEntries(); // necessary ????
 328   // Back to the signal file
 329   ((gAlice->TreeK())->GetCurrentFile())->cd();
 330   if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
 331   return;
 332 }
 333
 334 //__________________________________________________________________________
 335 void AliMUONEventReconstructor::EventReconstruct(void)
 336 {
 337   // To reconstruct one event
 338   if (fPrintLevel >= 1) cout << "enter EventReconstruct" << endl;
 339   MakeEventToBeReconstructed();
 340   MakeSegments();
 341   MakeTracks();
 342   return;
 343 }
 344
 345   //__________________________________________________________________________
 346 void AliMUONEventReconstructor::ResetHitsForRec(void)
 347 {
 348   // To reset the array and the number of HitsForRec,
 349   // and also the number of HitsForRec
 350   // and the index of the first HitForRec per chamber
 351   if (fHitsForRecPtr) fHitsForRecPtr->Clear();
 352   fNHitsForRec = 0;
 353   for (Int_t ch = 0; ch < MAX_MUON_TRACKING_CHAMBERS; ch++)
 354     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
 355   return;
 356 }
 357
 358   //__________________________________________________________________________
 359 void AliMUONEventReconstructor::ResetSegments(void)
 360 {
 361   // To reset the TClonesArray of segments and the number of Segments
 362   // for all stations
 363   for (Int_t st = 0; st < MAX_MUON_TRACKING_STATIONS; st++) {
 364     if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
 365     fNSegments[st] = 0;
 366   }
 367   return;
 368 }
 369
 370   //__________________________________________________________________________
 371 void AliMUONEventReconstructor::ResetTracks(void)
 372 {
 373   // To reset the TClonesArray of reconstructed tracks
 374   if (fRecTracksPtr) fRecTracksPtr->Clear();
 375   fNRecTracks = 0;
 376   return;
 377 }
 378
 379   //__________________________________________________________________________
 380 void AliMUONEventReconstructor::MakeEventToBeReconstructed(void)
 381 {
 382   // To make the list of hits to be reconstructed,
 383   // either from the GEANT hits or from the raw clusters
 384   // according to the parameter set for the reconstructor
 385   if (fPrintLevel >= 1) cout << "enter MakeEventToBeReconstructed" << endl;
 386   ResetHitsForRec();
 387   if (fRecGeantHits == 1) {
 388     // Reconstruction from GEANT hits
 389     // Back to the signal file
 390     ((gAlice->TreeK())->GetCurrentFile())->cd();
 391     // Signal hits
 392     // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 393     // Security on MUON ????
 394     AddHitsForRecFromGEANT(gAlice->TreeH());
 395     // Background hits
 396     AddHitsForRecFromBkgGEANT(fBkgGeantTH, fBkgGeantHits);
 397     // Sort HitsForRec in increasing order with respect to chamber number
 398     SortHitsForRecWithIncreasingChamber();
 399   }
 400   else {
 401     // Reconstruction from raw clusters
 402     // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 403     // Security on MUON ????
 404     // TreeR assumed to be be "prepared" in calling function
 405     // by "MUON->GetTreeR(nev)" ????
 406     TTree *TR = gAlice->TreeR();
 407     AddHitsForRecFromRawClusters(TR);
 408     // No sorting: it is done automatically in the previous function
 409   }
 410   if (fPrintLevel >= 10) {
 411     cout << "end of MakeEventToBeReconstructed" << endl;
 412     cout << "NHitsForRec: " << fNHitsForRec << endl;
 413     for (Int_t ch = 0; ch < MAX_MUON_TRACKING_CHAMBERS; ch++) {
 414       cout << "chamber(0...): " << ch
 415            << "  NHitsForRec: " << fNHitsForRecPerChamber[ch]
 416            << "  index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch]
 417            << endl;
 418       for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
 419            hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
 420            hit++) {
 421         cout << "HitForRec index(0...): " << hit << endl;
 422         ((*fHitsForRecPtr)[hit])->Dump();
 423       }
 424     }
 425   }
 426   return;
 427 }
 428
 429   //__________________________________________________________________________
 430 void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH)
 431 {
 432   // To add to the list of hits for reconstruction
 433   // the GEANT signal hits from a hit tree TH.
 434   if (fPrintLevel >= 2)
 435     cout << "enter AddHitsForRecFromGEANT with TH: " << TH << endl;
 436   if (TH == NULL) return;
 437   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 438   // Security on MUON ????
 439   // See whether it could be the same for signal and background ????
 440   // Loop over tracks in tree
 441   Int_t ntracks = (Int_t) TH->GetEntries();
 442   if (fPrintLevel >= 2)
 443     cout << "ntracks: " << ntracks << endl;
 444   for (Int_t track = 0; track < ntracks; track++) {
 445     gAlice->ResetHits();
 446     TH->GetEvent(track);
 447     // Loop over hits
 448     Int_t hit = 0;
 449     for (AliMUONHit* mHit = (AliMUONHit*) MUON->FirstHit(-1);
 450          mHit;
 451          mHit = (AliMUONHit*) MUON->NextHit(), hit++) {
 452       NewHitForRecFromGEANT(mHit,track, hit, 1);
 453     } // end of hit loop
 454   } // end of track loop
 455   return;
 456 }
 457
 458   //__________________________________________________________________________
 459 void AliMUONEventReconstructor::AddHitsForRecFromBkgGEANT(TTree *TH, TClonesArray *Hits)
 460 {
 461   // To add to the list of hits for reconstruction
 462   // the GEANT background hits from a hit tree TH and a pointer Hits to a hit list.
 463   // How to have only one function "AddHitsForRecFromGEANT" ????
 464   if (fPrintLevel >= 2)
 465     cout << "enter AddHitsForRecFromBkgGEANT with TH: " << TH << endl;
 466   if (TH == NULL) return;
 467   // Loop over tracks in tree
 468   Int_t ntracks = (Int_t) TH->GetEntries();
 469   if (fPrintLevel >= 2)
 470     cout << "ntracks: " << ntracks << endl;
 471   for (Int_t track = 0; track < ntracks; track++) {
 472     if (Hits) Hits->Clear();
 473     TH->GetEvent(track);
 474     // Loop over hits
 475     for (Int_t hit = 0; hit < Hits->GetEntriesFast(); hit++) {
 476       NewHitForRecFromGEANT((AliMUONHit*) (*Hits)[hit], track, hit, 0);
 477     } // end of hit loop
 478   } // end of track loop
 479   return;
 480 }
 481
 482   //__________________________________________________________________________
 483 AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* Hit, Int_t TrackNumber, Int_t HitNumber, Int_t Signal)
 484 {
 485   // To make a new hit for reconstruction from a GEANT hit pointed to by "Hit",
 486   // with hit number "HitNumber" in the track numbered "TrackNumber",
 487   // either from signal ("Signal" = 1) or background ("Signal" = 0) event.
 488   // Selects hits in tracking (not trigger) chambers.
 489   // Takes into account the efficiency (fEfficiency)
 490   // and the smearing from resolution (fBendingResolution and fNonBendingResolution).
 491   // Adds a condition on the radius between RMin and RMax
 492   // to better simulate the real chambers.
 493   // Returns the pointer to the new hit for reconstruction,
 494   // or NULL in case of inefficiency or non tracking chamber or bad radius.
 495   // No condition on at most 20 cm from a muon from a resonance
 496   // like in Fortran TRACKF_STAT.
 497   AliMUONHitForRec* hitForRec;
 498   Double_t bendCoor, nonBendCoor, radius;
 499   Int_t chamber = Hit->fChamber - 1; // chamber(0...)
 500   // only in tracking chambers (fChamber starts at 1)
 501   if (chamber >= MAX_MUON_TRACKING_CHAMBERS) return NULL;
 502   // only if hit is efficient (keep track for checking ????)
 503   if (gRandom->Rndm() > fEfficiency) return NULL;
 504   // only if radius between RMin and RMax
 505   bendCoor = Hit->fY;
 506   nonBendCoor = Hit->fX;
 507   radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor));
 508   if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
 509   // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's
 510   hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit);
 511   fNHitsForRec++;
 512   // add smearing from resolution
 513   hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution));
 514   hitForRec->SetNonBendingCoor(nonBendCoor
 515                                + gRandom->Gaus(0., fNonBendingResolution));
 516   // more information into HitForRec
 517   //  resolution: angular effect to be added here ????
 518   hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
 519   hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
 520   //  GEANT track info
 521   hitForRec->SetHitNumber(HitNumber);
 522   hitForRec->SetTHTrack(TrackNumber);
 523   hitForRec->SetGeantSignal(Signal);
 524   if (fPrintLevel >= 10) {
 525     cout << "track: " << TrackNumber << " hit: " << HitNumber << endl;
 526     Hit->Dump();
 527     cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 528     hitForRec->Dump();}
 529   return hitForRec;
 530 }
 531
 532   //__________________________________________________________________________
 533 void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber()
 534 {
 535   // Sort HitsForRec's in increasing order with respect to chamber number.
 536   // Uses the function "Compare".
 537   // Update the information for HitsForRec per chamber too.
 538   Int_t ch, nhits, prevch;
 539   fHitsForRecPtr->Sort();
 540   for (ch = 0; ch < MAX_MUON_TRACKING_CHAMBERS; ch++) {
 541     fNHitsForRecPerChamber[ch] = 0;
 542     fIndexOfFirstHitForRecPerChamber[ch] = 0;
 543   }
 544   prevch = 0; // previous chamber
 545   nhits = 0; // number of hits in current chamber
 546   // Loop over HitsForRec
 547   for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
 548     // chamber number (0...)
 549     ch = ((AliMUONHitForRec*)  ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
 550     // increment number of hits in current chamber
 551     (fNHitsForRecPerChamber[ch])++;
 552     // update index of first HitForRec in current chamber
 553     // if chamber number different from previous one
 554     if (ch != prevch) {
 555       fIndexOfFirstHitForRecPerChamber[ch] = hit;
 556       prevch = ch;
 557     }
 558   }
 559   return;
 560 }
 561
 562 //   //__________________________________________________________________________
 563 // void AliMUONEventReconstructor::AddHitsForRecFromCathodeCorrelations(TTree* TC)
 564 // {
 565 //   // OLD VERSION WHEN ONE ONE WAS USING SO CALLED CATHODE CORRELATIONS
 566 //   // To add to the list of hits for reconstruction
 567 //   // the (cathode correlated) raw clusters
 568 //   // No condition added, like in Fortran TRACKF_STAT,
 569 //   // on the radius between RMin and RMax.
 570 //   AliMUONHitForRec *hitForRec;
 571 //   if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromCathodeCorrelations" << endl;
 572 //   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 573 //   // Security on MUON ????
 574 //   // Loop over tracking chambers
 575 //   for (Int_t ch = 0; ch < MAX_MUON_TRACKING_CHAMBERS; ch++) {
 576 //     // number of HitsForRec to 0 for the chamber
 577 //     fNHitsForRecPerChamber[ch] = 0;
 578 //     // index of first HitForRec for the chamber
 579 //     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
 580 //     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
 581 //     TClonesArray *reconst_hits  = MUON->ReconstHitsAddress(ch);
 582 //     MUON->ResetReconstHits();
 583 //     TC->GetEvent();
 584 //     Int_t ncor = (Int_t)reconst_hits->GetEntries();
 585 //     // Loop over (cathode correlated) raw clusters
 586 //     for (Int_t cor = 0; cor < ncor; cor++) {
 587 //       AliMUONReconstHit * mCor =
 588 //      (AliMUONReconstHit*) reconst_hits->UncheckedAt(cor);
 589 //       // new AliMUONHitForRec from (cathode correlated) raw cluster
 590 //       // and increment number of AliMUONHitForRec's (total and in chamber)
 591 //       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(mCor);
 592 //       fNHitsForRec++;
 593 //       (fNHitsForRecPerChamber[ch])++;
 594 //       // more information into HitForRec
 595 //       hitForRec->SetChamberNumber(ch);
 596 //       hitForRec->SetHitNumber(cor);
 597 //       // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
 598 //       // could (should) be more exact from chamber geometry ????
 599 //       hitForRec->SetZ(-(&(MUON->Chamber(ch)))->Z());
 600 //       if (fPrintLevel >= 10) {
 601 //      cout << "chamber (0...): " << ch <<
 602 //        " cathcorrel (0...): " << cor << endl;
 603 //      mCor->Dump();
 604 //      cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 605 //      hitForRec->Dump();}
 606 //     } // end of cluster loop
 607 //   } // end of chamber loop
 608 //   return;
 609 // }
 610
 611   //__________________________________________________________________________
 612 void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
 613 {
 614   // To add to the list of hits for reconstruction all the raw clusters
 615   // No condition added, like in Fortran TRACKF_STAT,
 616   // on the radius between RMin and RMax.
 617   AliMUONHitForRec *hitForRec;
 618   AliMUONRawCluster *clus;
 619   Int_t iclus, nclus;
 620   TClonesArray *rawclusters;
 621   if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromRawClusters" << endl;
 622   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 623   // Security on MUON ????
 624   // Loop over tracking chambers
 625   for (Int_t ch = 0; ch < MAX_MUON_TRACKING_CHAMBERS; ch++) {
 626     // number of HitsForRec to 0 for the chamber
 627     fNHitsForRecPerChamber[ch] = 0;
 628     // index of first HitForRec for the chamber
 629     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
 630     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
 631     rawclusters = MUON->RawClustAddress(ch);
 632     MUON->ResetRawClusters();
 633     TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ????
 634     nclus = (Int_t) (rawclusters->GetEntries());
 635     // Loop over (cathode correlated) raw clusters
 636     for (iclus = 0; iclus < nclus; iclus++) {
 637       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
 638       // new AliMUONHitForRec from raw cluster
 639       // and increment number of AliMUONHitForRec's (total and in chamber)
 640       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
 641       fNHitsForRec++;
 642       (fNHitsForRecPerChamber[ch])++;
 643       // more information into HitForRec
 644       //  resolution: info should be already in raw cluster and taken from it ????
 645       hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
 646       hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
 647       //  original raw cluster
 648       hitForRec->SetChamberNumber(ch);
 649       hitForRec->SetHitNumber(iclus);
 650       // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
 651       // could (should) be more exact from chamber geometry ????
 652       hitForRec->SetZ(-(&(MUON->Chamber(ch)))->Z());
 653       if (fPrintLevel >= 10) {
 654         cout << "chamber (0...): " << ch <<
 655           " raw cluster (0...): " << iclus << endl;
 656         clus->Dump();
 657         cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
 658         hitForRec->Dump();}
 659     } // end of cluster loop
 660   } // end of chamber loop
 661   return;
 662 }
 663
 664   //__________________________________________________________________________
 665 void AliMUONEventReconstructor::MakeSegments(void)
 666 {
 667   // To make the list of segments in all stations,
 668   // from the list of hits to be reconstructed
 669   if (fPrintLevel >= 1) cout << "enter MakeSegments" << endl;
 670   ResetSegments();
 671   // Loop over stations
 672   for (Int_t st = 0; st < MAX_MUON_TRACKING_STATIONS; st++)
 673     MakeSegmentsPerStation(st);
 674   if (fPrintLevel >= 10) {
 675     cout << "end of MakeSegments" << endl;
 676     for (Int_t st = 0; st < MAX_MUON_TRACKING_STATIONS; st++) {
 677       cout << "station(0...): " << st
 678            << "  Segments: " << fNSegments[st]
 679            << endl;
 680       for (Int_t seg = 0;
 681            seg < fNSegments[st];
 682            seg++) {
 683         cout << "Segment index(0...): " << seg << endl;
 684         ((*fSegmentsPtr[st])[seg])->Dump();
 685       }
 686     }
 687   }
 688   return;
 689 }
 690
 691   //__________________________________________________________________________
 692 void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
 693 {
 694   // To make the list of segments in station number "Station" (0...)
 695   // from the list of hits to be reconstructed.
 696   // Updates "fNSegments"[Station].
 697   // Segments in stations 4 and 5 are sorted
 698   // according to increasing absolute value of "impact parameter"
 699   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
 700   AliMUONSegment *segment;
 701   Bool_t last2st;
 702   Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
 703     impactParam, maxImpactParam;
 704   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 705   if (fPrintLevel >= 1)
 706     cout << "enter MakeSegmentsPerStation (0...) " << Station << endl;
 707   // first and second chambers (0...) in the station
 708   Int_t ch1 = 2 * Station;
 709   Int_t ch2 = ch1 + 1;
 710   // variable true for stations downstream of the dipole:
 711   // Station(0..4) equal to 3 or 4
 712   if ((Station == 3) || (Station == 4)) {
 713     last2st = kTRUE;
 714     // maximum impact parameter (cm) according to fMinBendingMomentum...
 715     maxImpactParam =
 716       TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
 717   }
 718   else last2st = kFALSE;
 719   // extrapolation factor from Z of first chamber to Z of second chamber
 720   // dZ to be changed to take into account fine structure of chambers ????
 721   Double_t extrapFact =
 722     (&(MUON->Chamber(ch2)))->Z() / (&(MUON->Chamber(ch1)))->Z();
 723   // index for current segment
 724   Int_t segmentIndex = 0;
 725   // Loop over HitsForRec in the first chamber of the station
 726   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
 727        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
 728        hit1++) {
 729     // pointer to the HitForRec
 730     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
 731     // extrapolation,
 732     // on the straight line joining the HitForRec to the vertex (0,0,0),
 733     // to the Z of the second chamber of the station
 734     extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
 735     extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
 736     // Loop over HitsForRec in the second chamber of the station
 737     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
 738          hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
 739          hit2++) {
 740       // pointer to the HitForRec
 741       hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
 742       // absolute values of distances, in bending and non bending planes,
 743       // between the HitForRec in the second chamber
 744       // and the previous extrapolation
 745       distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
 746       distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
 747       if (last2st) {
 748         // bending slope
 749         bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
 750           (hit1Ptr->GetZ() - hit2Ptr->GetZ());
 751         // absolute value of impact parameter
 752         impactParam =
 753           TMath::Abs(hit1Ptr->GetBendingCoor() - hit2Ptr->GetZ() * bendingSlope);
 754       }
 755       // check for distances not too large,
 756       // and impact parameter not too big if stations downstream of the dipole.
 757       // Conditions "distBend" and "impactParam" correlated for these stations ????
 758       if ((distBend < fSegmentMaxDistBending[Station]) &&
 759           (distNonBend < fSegmentMaxDistNonBending[Station]) &&
 760           (!last2st || (impactParam < maxImpactParam))) {
 761         // make new segment
 762         segment = new ((*fSegmentsPtr[Station])[segmentIndex])
 763           AliMUONSegment(hit1Ptr, hit2Ptr);
 764         // update "link" to this segment from the hit in the first chamber
 765         if (hit1Ptr->GetNSegments() == 0)
 766           hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
 767         hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
 768         if (fPrintLevel >= 10) {
 769           cout << "segmentIndex(0...): " << segmentIndex
 770                << "  distBend: " << distBend
 771                << "  distNonBend: " << distNonBend
 772                << endl;
 773           segment->Dump();
 774           cout << "HitForRec in first chamber" << endl;
 775           hit1Ptr->Dump();
 776           cout << "HitForRec in second chamber" << endl;
 777           hit2Ptr->Dump();
 778         };
 779         // increment index for current segment
 780         segmentIndex++;
 781       }
 782     } //for (Int_t hit2
 783   } // for (Int_t hit1...
 784   fNSegments[Station] = segmentIndex;
 785   // Sorting according to "impact parameter" if station(1..5) 4 or 5,
 786   // i.e. Station(0..4) 3 or 4, using the function "Compare".
 787   // After this sorting, it is impossible to use
 788   // the "fNSegments" and "fIndexOfFirstSegment"
 789   // of the HitForRec in the first chamber to explore all segments formed with it.
 790   // Is this sorting really needed ????
 791   if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
 792   if (fPrintLevel >= 1) cout << "Station: " << Station << "  NSegments: "
 793                              << fNSegments[Station] << endl;
 794   return;
 795 }
 796
 797   //__________________________________________________________________________
 798 void AliMUONEventReconstructor::MakeTracks(void)
 799 {
 800   // To make the tracks,
 801   // from the list of segments and points in all stations
 802   if (fPrintLevel >= 1) cout << "enter MakeTracks" << endl;
 803   ResetTracks();
 804   // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
 805   MakeTrackCandidates();
 806   // Follow tracks in stations(1..) 3, 2 and 1
 807   FollowTracks();
 808   return;
 809 }
 810
 811   //__________________________________________________________________________
 812 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
 813 {
 814   // To make track candidates with two segments in stations(1..) 4 and 5,
 815   // the first segment being pointed to by "BegSegment".
 816   // Returns the number of such track candidates.
 817   Int_t endStation, iEndSegment, nbCan2Seg;
 818   AliMUONSegment *endSegment, *extrapSegment;
 819   AliMUONTrack *recTrack;
 820   Double_t mcsFactor;
 821   if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesWithTwoSegments" << endl;
 822   // Station for the end segment
 823   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 824   // multiple scattering factor corresponding to one chamber
 825   mcsFactor = 0.0136 /
 826     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 827   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 828   // linear extrapolation to end station
 829   extrapSegment =
 830     BegSegment->CreateSegmentFromLinearExtrapToStation(endStation, mcsFactor);
 831   // number of candidates with 2 segments to 0
 832   nbCan2Seg = 0;
 833   // Loop over segments in the end station
 834   for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
 835     // pointer to segment
 836     endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
 837     // test compatibility between current segment and "extrapSegment"
 838     if ((endSegment->
 839          NormalizedChi2WithSegment(extrapSegment,
 840                                    fMaxSigma2Distance)) <= 4.0) {
 841       // both segments compatible:
 842       // make track candidate from "begSegment" and "endSegment"
 843       if (fPrintLevel >= 2)
 844         cout << "TrackCandidate with Segment " << iEndSegment <<
 845           " in Station(0..) " << endStation << endl;
 846       // flag for both segments in one track:
 847       // to be done in track constructor ????
 848       BegSegment->SetInTrack(kTRUE);
 849       endSegment->SetInTrack(kTRUE);
 850       recTrack = new ((*fRecTracksPtr)[fNRecTracks])
 851         AliMUONTrack(BegSegment, endSegment, this);
 852       fNRecTracks++;
 853       if (fPrintLevel >= 10) recTrack->RecursiveDump();
 854       // increment number of track candidates with 2 segments
 855       nbCan2Seg++;
 856     }
 857   } // for (iEndSegment = 0;...
 858   delete extrapSegment; // should not delete HitForRec's it points to !!!!
 859   return nbCan2Seg;
 860 }
 861
 862   //__________________________________________________________________________
 863 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
 864 {
 865   // To make track candidates with one segment and one point
 866   // in stations(1..) 4 and 5,
 867   // the segment being pointed to by "BegSegment".
 868   Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
 869   AliMUONHitForRec *extrapHitForRec, *hit;
 870   AliMUONTrack *recTrack;
 871   Double_t mcsFactor;
 872   if (fPrintLevel >= 1)
 873     cout << "enter MakeTrackCandidatesWithOneSegmentAndOnePoint" << endl;
 874   // station for the end point
 875   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
 876   // multiple scattering factor corresponding to one chamber
 877   mcsFactor = 0.0136 /
 878     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
 879   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 880   // first and second chambers(0..) in the end station
 881   ch1 = 2 * endStation;
 882   ch2 = ch1 + 1;
 883   // number of candidates to 0
 884   nbCan1Seg1Hit = 0;
 885   // Loop over chambers of the end station
 886   for (ch = ch2; ch >= ch1; ch--) {
 887     // linear extrapolation to chamber
 888     extrapHitForRec =
 889       BegSegment->CreateHitForRecFromLinearExtrapToChamber(ch, mcsFactor);
 890     // limits for the hit index in the loop
 891     iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
 892     iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
 893     // Loop over HitForRec's in the chamber
 894     for (iHit = iHitMin; iHit < iHitMax; iHit++) {
 895       // pointer to HitForRec
 896       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
 897       // test compatibility between current HitForRec and "extrapHitForRec"
 898       if ((hit->
 899            NormalizedChi2WithHitForRec(extrapHitForRec,
 900                                        fMaxSigma2Distance)) <= 2.0) {
 901         // both HitForRec's compatible:
 902         // make track candidate from begSegment and current HitForRec
 903         if (fPrintLevel >= 2)
 904           cout << "TrackCandidate with HitForRec " << iHit <<
 905             " in Chamber(0..) " << ch << endl;
 906         // flag for beginning segments in one track:
 907         // to be done in track constructor ????
 908         BegSegment->SetInTrack(kTRUE);
 909         recTrack = new ((*fRecTracksPtr)[fNRecTracks])
 910           AliMUONTrack(BegSegment, hit, this);
 911         // the right place to eliminate "double counting" ???? how ????
 912         fNRecTracks++;
 913         if (fPrintLevel >= 10) recTrack->RecursiveDump();
 914         // increment number of track candidates
 915         nbCan1Seg1Hit++;
 916       }
 917     } // for (iHit = iHitMin;...
 918     delete extrapHitForRec;
 919   } // for (ch = ch2;...
 920   return nbCan1Seg1Hit;
 921 }
 922
 923   //__________________________________________________________________________
 924 void AliMUONEventReconstructor::MakeTrackCandidates(void)
 925 {
 926   // To make track candidates
 927   // with at least 3 aligned points in stations(1..) 4 and 5
 928   // (two Segment's or one Segment and one HitForRec)
 929   Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
 930   AliMUONSegment *begSegment;
 931   if (fPrintLevel >= 1) cout << "enter MakeTrackCandidates" << endl;
 932   // Loop over stations(1..) 5 and 4 for the beginning segment
 933   for (begStation = 4; begStation > 2; begStation--) {
 934     // Loop over segments in the beginning station
 935     for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
 936       // pointer to segment
 937       begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
 938       if (fPrintLevel >= 2)
 939         cout << "look for TrackCandidate's with Segment " << iBegSegment <<
 940           " in Station(0..) " << begStation << endl;
 941       // Look for track candidates with two segments,
 942       // "begSegment" and all compatible segments in other station.
 943       // Only for beginning station(1..) 5
 944       // because candidates with 2 segments have to looked for only once.
 945       if (begStation == 4)
 946         nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
 947       // Look for track candidates with one segments and one point,
 948       // "begSegment" and all compatible HitForRec's in other station.
 949       // Only if "begSegment" does not belong already to a track candidate.
 950       // Is that a too strong condition ????
 951       if (!(begSegment->GetInTrack()))
 952         nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
 953     } // for (iBegSegment = 0;...
 954   } // for (begStation = 4;...
 955   return;
 956 }
 957
 958   //__________________________________________________________________________
 959 void AliMUONEventReconstructor::FollowTracks(void)
 960 {
 961   // Follow tracks in stations(1..) 3, 2 and 1
 962   AliMUONHitForRec *bestHit, *extrapHit, *hit;
 963   AliMUONSegment *bestSegment, *extrapSegment, *segment;
 964   AliMUONTrack *track;
 965   AliMUONTrackParam *trackParam1, trackParam[2];
 966   Int_t ch, chBestHit, iHit, iSegment, station, trackIndex;
 967   Double_t bestChi2, chi2, dZ1, dZ2, maxSigma2Distance, mcsFactor;
 968   AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
 969   maxSigma2Distance = 4. * fMaxSigma2Distance; // sigma2cut increased by 4 !!!!
 970   if (fPrintLevel >= 2)
 971     cout << "enter FollowTracks" << endl;
 972   // Loop over track candidates
 973   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
 974     if (fPrintLevel >= 2)
 975       cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
 976     // function for each track candidate ????
 977     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
 978     // Fit track candidate from vertex at X = Y = 0
 979     track->Fit(track->GetTrackParamAtVertex(), 3);
 980     if (fPrintLevel >= 10) {
 981       cout << "FollowTracks: track candidate(0..): " << trackIndex
 982            << " after fit in stations(1..) 4 and 5" << endl;
 983       track->RecursiveDump();
 984     }
 985     // Loop over stations(1..) 3, 2 and 1
 986     // something SPECIAL for stations 2 and 1 for majority coincidence ????
 987     for (station = 2; station >= 0; station--) {
 988       // Track parameters at first track hit (smallest Z)
 989       trackParam1 = ((AliMUONTrackHit*)
 990                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
 991       // extrapolation to station
 992       trackParam1->ExtrapToStation(station, trackParam);
 993       extrapSegment = new AliMUONSegment(); //  empty segment
 994       // multiple scattering factor corresponding to one chamber
 995       // and momentum in bending plane (not total)
 996       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
 997       mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
 998       // Z difference from previous station
 999       dZ1 = (&(MUON->Chamber(2 * station)))->Z() -
1000         (&(MUON->Chamber(2 * station + 2)))->Z();
1001       // Z difference between the two previous stations
1002       dZ2 = (&(MUON->Chamber(2 * station + 2)))->Z() -
1003         (&(MUON->Chamber(2 * station + 4)))->Z();
1004       extrapSegment->SetBendingCoorReso2(fBendingResolution);
1005       extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution);
1006       extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2);
1007       bestChi2 = 5.0;
1008       bestSegment = NULL;
1009       if (fPrintLevel >= 10) {
1010         cout << "FollowTracks: track candidate(0..): " << trackIndex
1011              << " Look for segment in station(0..): " << station << endl;
1012       }
1013       // Loop over segments in station
1014       for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
1015         // Look for best compatible Segment in station
1016         // should consider all possibilities ????
1017         // multiple scattering ????
1018         // separation in 2 functions: Segment and HitForRec ????
1019         segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
1020         chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
1021         if (chi2 < bestChi2) {
1022           // update best Chi2 and Segment if better found
1023           bestSegment = segment;
1024           bestChi2 = chi2;
1025         }
1026       }
1027       if (bestSegment) {
1028         // best segment found: add it to track candidate
1029         track->AddSegment(bestSegment);
1030         // set track parameters at these two TrakHit's
1031         track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0]));
1032         track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1]));
1033         if (fPrintLevel >= 10) {
1034           cout << "FollowTracks: track candidate(0..): " << trackIndex
1035                << " Added segment in station(0..): " << station << endl;
1036           track->RecursiveDump();
1037         }
1038       }
1039       else {
1040         // No best segment found:
1041         // Look for best compatible HitForRec in station:
1042         // should consider all possibilities ????
1043         // multiple scattering ???? do about like for extrapSegment !!!!
1044         extrapHit = new AliMUONHitForRec(); //  empty hit
1045         bestChi2 = 3.0;
1046         bestHit = NULL;
1047         if (fPrintLevel >= 10) {
1048           cout << "FollowTracks: track candidate(0..): " << trackIndex
1049                << " Segment not found, look for hit in station(0..): " << station
1050                << endl;
1051         }
1052         // Loop over chambers of the station
1053         for (ch = 0; ch < 2; ch++) {
1054           // coordinates of extrapolated hit
1055           extrapHit->SetBendingCoor((&(trackParam[ch]))->GetBendingCoor());
1056           extrapHit->SetNonBendingCoor((&(trackParam[ch]))->GetNonBendingCoor());
1057           // resolutions from "extrapSegment"
1058           extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
1059           extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
1060           // Loop over hits in the chamber
1061           for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
1062                iHit < fIndexOfFirstHitForRecPerChamber[ch] +
1063                  fNHitsForRecPerChamber[ch];
1064                iHit++) {
1065             hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
1066             // condition for hit not already in segment ????
1067             chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
1068             if (chi2 < bestChi2) {
1069               // update best Chi2 and HitForRec if better found
1070               bestHit = hit;
1071               bestChi2 = chi2;
1072               chBestHit = ch;
1073             }
1074           }
1075         }
1076         if (bestHit) {
1077           // best hit found: add it to track candidate
1078           track->AddHitForRec(bestHit);
1079           // set track parameters at these two TrakHit's
1080           track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
1081                                     &(trackParam[chBestHit]));
1082           if (fPrintLevel >= 10) {
1083             cout << "FollowTracks: track candidate(0..): " << trackIndex
1084                  << " Added hit in station(0..): " << station << endl;
1085             track->RecursiveDump();
1086           }
1087         }
1088         else {
1089           fRecTracksPtr->RemoveAt(trackIndex); // Remove candidate
1090           break; // stop the search for this candidate:
1091           // exit from the loop over station
1092         }
1093       }
1094       // Sort track hits according to increasing Z
1095       track->GetTrackHitsPtr()->Sort();
1096       // Update track parameters at first track hit (smallest Z)
1097       trackParam1 = ((AliMUONTrackHit*)
1098                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1099       // Track fit from first track hit varying X and Y
1100       track->Fit(trackParam1, 5);
1101       if (fPrintLevel >= 10) {
1102         cout << "FollowTracks: track candidate(0..): " << trackIndex
1103              << " after fit from station(0..): " << station << " to 4" << endl;
1104         track->RecursiveDump();
1105       }
1106       delete extrapSegment;
1107     } // for (station = 2;...
1108   } // for (trackIndex = 0;...
1109   return;
1110 }
1111
1112 void AliMUONEventReconstructor::Streamer(TBuffer &R__b)
1113 {
1114     ;
1115 }