Using AliTracker (external mag field) instead of gAlice->Field (Nicolas & Christian)
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackReconstructor.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 in ALICE (class renamed from AliMUONEventReconstructor)
21 //
22 // This class contains as data:
23 // * the parameters for the track reconstruction
24 // * a pointer to the array of hits to be reconstructed (the event)
25 // * a pointer to the array of segments made with these hits inside each station
26 // * a pointer to the array of reconstructed tracks
27 //
28 // It contains as methods, among others:
29 // * MakeEventToBeReconstructed to build the array of hits to be reconstructed
30 // * MakeSegments to build the segments
31 // * MakeTracks to build the tracks
32 //
33 ////////////////////////////////////
34
35 #include <stdlib.h>
36 #include <Riostream.h>
37 #include <TDirectory.h>
38 #include <TFile.h>
39 #include <TMatrixD.h>
40
41 #include "AliMUONTrackReconstructor.h"
42 #include "AliMUON.h"
43 #include "AliMUONConstants.h"
44 #include "AliMUONHitForRec.h"
45 #include "AliMUONTriggerTrack.h"
46 #include "AliMUONTriggerCircuitNew.h"
47 #include "AliMUONRawCluster.h"
48 #include "AliMUONLocalTrigger.h"
49 #include "AliMUONGlobalTrigger.h"
50 #include "AliMUONSegment.h"
51 #include "AliMUONTrack.h"
52 #include "AliMUONTrackHit.h"
53 #include "AliMagF.h"
54 #include "AliRun.h" // for gAlice
55 #include "AliRunLoader.h"
56 #include "AliLoader.h"
57 #include "AliMUONTrackK.h" 
58 #include "AliLog.h"
59 #include "AliTracker.h"
60
61 //************* Defaults parameters for reconstruction
62 const Double_t AliMUONTrackReconstructor::fgkDefaultMinBendingMomentum = 3.0;
63 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxBendingMomentum = 3000.0;
64 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxChi2 = 100.0;
65 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxSigma2Distance = 16.0;
66 const Double_t AliMUONTrackReconstructor::fgkDefaultBendingResolution = 0.01;
67 const Double_t AliMUONTrackReconstructor::fgkDefaultNonBendingResolution = 0.144;
68 const Double_t AliMUONTrackReconstructor::fgkDefaultChamberThicknessInX0 = 0.03;
69 // Simple magnetic field:
70 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
71 // Length and Position from reco_muon.F, with opposite sign:
72 // Length = ZMAGEND-ZCOIL
73 // Position = (ZMAGEND+ZCOIL)/2
74 // to be ajusted differently from real magnetic field ????
75 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBValue = 7.0;
76 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBLength = 428.0;
77 const Double_t AliMUONTrackReconstructor::fgkDefaultSimpleBPosition = 1019.0;
78 const Double_t AliMUONTrackReconstructor::fgkDefaultEfficiency = 0.95;
79
80 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
81
82 //__________________________________________________________________________
83 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliLoader* loader, AliMUONData* data)
84   : TObject(),
85     fTrackMethod(1), //AZ - tracking method (1-default, 2-Kalman)
86     fMinBendingMomentum(fgkDefaultMinBendingMomentum),
87     fMaxBendingMomentum(fgkDefaultMaxBendingMomentum),
88     fMaxChi2(fgkDefaultMaxChi2),
89     fMaxSigma2Distance(fgkDefaultMaxSigma2Distance),
90     fBendingResolution(fgkDefaultBendingResolution),
91     fNonBendingResolution(fgkDefaultNonBendingResolution),
92     fChamberThicknessInX0(fgkDefaultChamberThicknessInX0),
93     fSimpleBValue(fgkDefaultSimpleBValue),
94     fSimpleBLength(fgkDefaultSimpleBLength),
95     fSimpleBPosition(fgkDefaultSimpleBPosition),
96     fEfficiency(fgkDefaultEfficiency),
97     fHitsForRecPtr(0x0),
98     fNHitsForRec(0),
99     fRecTracksPtr(0x0),
100     fNRecTracks(0),
101     fRecTrackHitsPtr(0x0),
102     fNRecTrackHits(0),
103     fMUONData(data),
104     fLoader(loader),
105     fMuons(0),
106     fTriggerTrack(new AliMUONTriggerTrack())
107 {
108   // Constructor for class AliMUONTrackReconstructor
109   SetReconstructionParametersToDefaults();
110
111   // Memory allocation for the TClonesArray of hits for reconstruction
112   // Is 10000 the right size ????
113   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
114
115   // Memory allocation for the TClonesArray's of segments in stations
116   // Is 2000 the right size ????
117   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
118     fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
119     fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
120   }
121   // Memory allocation for the TClonesArray of reconstructed tracks
122   // Is 10 the right size ????
123   fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
124
125   // Memory allocation for the TClonesArray of hits on reconstructed tracks
126   // Is 100 the right size ????
127   fRecTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 100);
128
129   const AliMagF* kField = AliTracker::GetFieldMap();
130   if (!kField) AliFatal("No field available");
131  // Sign of fSimpleBValue according to sign of Bx value at (50,50,-950).
132   Float_t b[3], x[3];
133   x[0] = 50.; x[1] = 50.; x[2] = -950.;
134   kField->Field(x,b);
135
136   fSimpleBValue    = TMath::Sign(fSimpleBValue,(Double_t) b[0]);
137   fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]);
138   // See how to get fSimple(BValue, BLength, BPosition)
139   // automatically calculated from the actual magnetic field ????
140
141   return;
142 }
143
144   //__________________________________________________________________________
145 AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void)
146 {
147   // Destructor for class AliMUONTrackReconstructor
148   delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
149   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
150     delete fSegmentsPtr[st]; // Correct destruction of everything ????
151
152   delete fTriggerTrack;
153 }
154   //__________________________________________________________________________
155 void AliMUONTrackReconstructor::SetReconstructionParametersToDefaults(void)
156 {
157   // Set reconstruction parameters to default values
158   // Would be much more convenient with a structure (or class) ????
159
160   // ******** Parameters for making HitsForRec
161   // minimum radius,
162   // like in TRACKF_STAT:
163   // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
164   // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
165   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
166     if (ch < 4) fRMin[ch] = TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)) *
167                   2.0 * TMath::Pi() / 180.0;
168     else fRMin[ch] = 30.0;
169     // maximum radius at 10 degrees and Z of chamber
170     fRMax[ch] = TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)) *
171                   10.0 * TMath::Pi() / 180.0;
172   }
173
174   // ******** Parameters for making segments
175   // should be parametrized ????
176   // according to interval between chambers in a station ????
177   // Maximum distance in non bending plane
178   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
179   // SIGCUT*DYMAX(IZ)
180   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
181     fSegmentMaxDistNonBending[st] = 5. * 0.22;
182   // Maximum distance in bending plane:
183   // values from TRACKF_STAT, corresponding to (J psi 20cm),
184   // scaled to the real distance between chambers in a station
185   fSegmentMaxDistBending[0] = TMath::Abs( 1.5 *
186                                           (AliMUONConstants::DefaultChamberZ(1) - AliMUONConstants::DefaultChamberZ(0)) / 20.0);
187   fSegmentMaxDistBending[1] = TMath::Abs( 1.5 *
188                                           (AliMUONConstants::DefaultChamberZ(3) - AliMUONConstants::DefaultChamberZ(2)) / 20.0);
189   fSegmentMaxDistBending[2] = TMath::Abs( 3.0 *
190                                           (AliMUONConstants::DefaultChamberZ(5) - AliMUONConstants::DefaultChamberZ(4)) / 20.0);
191   fSegmentMaxDistBending[3] = TMath::Abs( 6.0 *
192                                           (AliMUONConstants::DefaultChamberZ(7) - AliMUONConstants::DefaultChamberZ(6)) / 20.0);
193   fSegmentMaxDistBending[4] = TMath::Abs( 6.0 *
194                                           (AliMUONConstants::DefaultChamberZ(9) - AliMUONConstants::DefaultChamberZ(8)) / 20.0);
195
196   return;
197 }
198
199 //__________________________________________________________________________
200 Double_t AliMUONTrackReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const
201 {
202   // Returns impact parameter at vertex in bending plane (cm),
203   // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
204   // using simple values for dipole magnetic field.
205   // The sign of "BendingMomentum" is the sign of the charge.
206   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
207           BendingMomentum);
208 }
209
210 //__________________________________________________________________________
211 Double_t AliMUONTrackReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam) const
212 {
213   // Returns signed bending momentum in bending plane (GeV/c),
214   // the sign being the sign of the charge for particles moving forward in Z,
215   // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
216   // using simple values for dipole magnetic field.
217   return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
218           ImpactParam);
219 }
220
221 //__________________________________________________________________________
222 void AliMUONTrackReconstructor::EventReconstruct(void)
223 {
224   // To reconstruct one event
225   AliDebug(1,"Enter EventReconstruct");
226   ResetTracks(); //AZ
227   ResetTrackHits(); //AZ 
228   ResetSegments(); //AZ
229   ResetHitsForRec(); //AZ
230   MakeEventToBeReconstructed();
231   MakeSegments();
232   MakeTracks();
233   if (fMUONData->IsTriggerTrackBranchesInTree()) 
234     ValidateTracksWithTrigger(); 
235
236   // Add tracks to MUON data container 
237   for(Int_t i=0; i<GetNRecTracks(); i++) {
238     AliMUONTrack * track = (AliMUONTrack*) GetRecTracksPtr()->At(i);
239     fMUONData->AddRecTrack(*track);
240   }
241
242   return;
243 }
244
245 //__________________________________________________________________________
246 void AliMUONTrackReconstructor::EventReconstructTrigger(void)
247 {
248   // To reconstruct one event
249   AliDebug(1,"Enter EventReconstructTrigger");
250   MakeTriggerTracks();  
251   return;
252 }
253
254   //__________________________________________________________________________
255 void AliMUONTrackReconstructor::ResetHitsForRec(void)
256 {
257   // To reset the array and the number of HitsForRec,
258   // and also the number of HitsForRec
259   // and the index of the first HitForRec per chamber
260   if (fHitsForRecPtr) fHitsForRecPtr->Delete();
261   fNHitsForRec = 0;
262   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++)
263     fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
264   return;
265 }
266
267   //__________________________________________________________________________
268 void AliMUONTrackReconstructor::ResetSegments(void)
269 {
270   // To reset the TClonesArray of segments and the number of Segments
271   // for all stations
272   for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
273     if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
274     fNSegments[st] = 0;
275   }
276   return;
277 }
278
279   //__________________________________________________________________________
280 void AliMUONTrackReconstructor::ResetTracks(void)
281 {
282   // To reset the TClonesArray of reconstructed tracks
283   if (fRecTracksPtr) fRecTracksPtr->Delete();
284   // Delete in order that the Track destructors are called,
285   // hence the space for the TClonesArray of pointers to TrackHit's is freed
286   fNRecTracks = 0;
287   return;
288 }
289
290
291   //__________________________________________________________________________
292 void AliMUONTrackReconstructor::ResetTrackHits(void)
293 {
294   // To reset the TClonesArray of hits on reconstructed tracks
295   if (fRecTrackHitsPtr) fRecTrackHitsPtr->Clear();
296   fNRecTrackHits = 0;
297   return;
298 }
299
300   //__________________________________________________________________________
301 void AliMUONTrackReconstructor::MakeEventToBeReconstructed(void)
302 {
303   // To make the list of hits to be reconstructed,
304   // either from the track ref. hits or from the raw clusters
305   // according to the parameter set for the reconstructor
306
307   AliRunLoader *runLoader = fLoader->GetRunLoader();
308
309   AliDebug(1,"Enter MakeEventToBeReconstructed");
310   //AZ ResetHitsForRec();
311  
312   // Reconstruction from raw clusters
313   // AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
314   // Security on MUON ????
315   // TreeR assumed to be be "prepared" in calling function
316   // by "MUON->GetTreeR(nev)" ????
317   TTree *treeR = fLoader->TreeR();
318
319   //AZ? fMUONData->SetTreeAddress("RC");
320   AddHitsForRecFromRawClusters(treeR);
321   // No sorting: it is done automatically in the previous function
322   
323  
324   AliDebug(1,"End of MakeEventToBeReconstructed");
325     if (AliLog::GetGlobalDebugLevel() > 0) {
326       AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
327       for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
328         AliDebug(1, Form("Chamber(0...): %d",ch));
329         AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
330         AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
331         for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
332              hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
333              hit++) {
334           AliDebug(1, Form("HitForRec index(0...): %d",hit));
335           ((*fHitsForRecPtr)[hit])->Dump();
336       }
337     }
338   }
339   return;
340 }
341
342   //__________________________________________________________________________
343 void AliMUONTrackReconstructor::SortHitsForRecWithIncreasingChamber()
344 {
345   // Sort HitsForRec's in increasing order with respect to chamber number.
346   // Uses the function "Compare".
347   // Update the information for HitsForRec per chamber too.
348   Int_t ch, nhits, prevch;
349   fHitsForRecPtr->Sort();
350   for (ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
351     fNHitsForRecPerChamber[ch] = 0;
352     fIndexOfFirstHitForRecPerChamber[ch] = 0;
353   }
354   prevch = 0; // previous chamber
355   nhits = 0; // number of hits in current chamber
356   // Loop over HitsForRec
357   for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
358     // chamber number (0...)
359     ch = ((AliMUONHitForRec*)  ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
360     // increment number of hits in current chamber
361     (fNHitsForRecPerChamber[ch])++;
362     // update index of first HitForRec in current chamber
363     // if chamber number different from previous one
364     if (ch != prevch) {
365       fIndexOfFirstHitForRecPerChamber[ch] = hit;
366       prevch = ch;
367     }
368   }
369   return;
370 }
371
372   //__________________________________________________________________________
373 void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
374 {
375   // To add to the list of hits for reconstruction all the raw clusters
376   // No condition added, like in Fortran TRACKF_STAT,
377   // on the radius between RMin and RMax.
378   AliMUONHitForRec *hitForRec;
379   AliMUONRawCluster *clus;
380   Int_t iclus, nclus, nTRentries;
381   TClonesArray *rawclusters;
382   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
383
384   if (fTrackMethod != 3) { //AZ
385     fMUONData->SetTreeAddress("RC"); //AZ
386     nTRentries = Int_t(TR->GetEntries());
387     if (nTRentries != 1) {
388       AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
389       exit(0);
390     }
391     fLoader->TreeR()->GetEvent(0); // only one entry  
392   }
393
394   // Loop over tracking chambers
395   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
396     // number of HitsForRec to 0 for the chamber
397     fNHitsForRecPerChamber[ch] = 0;
398     // index of first HitForRec for the chamber
399     if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
400     else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
401     rawclusters =fMUONData->RawClusters(ch);
402     nclus = (Int_t) (rawclusters->GetEntries());
403     // Loop over (cathode correlated) raw clusters
404     for (iclus = 0; iclus < nclus; iclus++) {
405       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
406       // new AliMUONHitForRec from raw cluster
407       // and increment number of AliMUONHitForRec's (total and in chamber)
408       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
409       fNHitsForRec++;
410       (fNHitsForRecPerChamber[ch])++;
411       // more information into HitForRec
412       //  resolution: info should be already in raw cluster and taken from it ????
413       //hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
414       //hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
415       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
416       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
417       //  original raw cluster
418       hitForRec->SetChamberNumber(ch);
419       hitForRec->SetHitNumber(iclus);
420       // Z coordinate of the raw cluster (cm)
421       hitForRec->SetZ(clus->GetZ(0));
422       if (AliLog::GetGlobalDebugLevel() > 1) {
423         cout << "chamber (0...): " << ch <<
424           " raw cluster (0...): " << iclus << endl;
425         clus->Dump();
426         cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
427         hitForRec->Dump();}
428     } // end of cluster loop
429   } // end of chamber loop
430   SortHitsForRecWithIncreasingChamber(); 
431   return;
432 }
433
434   //__________________________________________________________________________
435 void AliMUONTrackReconstructor::MakeSegments(void)
436 {
437   // To make the list of segments in all stations,
438   // from the list of hits to be reconstructed
439   AliDebug(1,"Enter MakeSegments");
440   //AZ ResetSegments();
441   // Loop over stations
442   Int_t nb = (fTrackMethod != 1) ? 3 : 0; //AZ
443   for (Int_t st = nb; st < AliMUONConstants::NTrackingCh()/2; st++) 
444     MakeSegmentsPerStation(st); 
445   if (AliLog::GetGlobalDebugLevel() > 1) {
446     cout << "end of MakeSegments" << endl;
447     for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
448       cout << "station(0...): " << st
449            << "  Segments: " << fNSegments[st]
450            << endl;
451       for (Int_t seg = 0;
452            seg < fNSegments[st];
453            seg++) {
454         cout << "Segment index(0...): " << seg << endl;
455         ((*fSegmentsPtr[st])[seg])->Dump();
456       }
457     }
458   }
459   return;
460 }
461
462   //__________________________________________________________________________
463 void AliMUONTrackReconstructor::MakeSegmentsPerStation(Int_t Station)
464 {
465   // To make the list of segments in station number "Station" (0...)
466   // from the list of hits to be reconstructed.
467   // Updates "fNSegments"[Station].
468   // Segments in stations 4 and 5 are sorted
469   // according to increasing absolute value of "impact parameter"
470   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
471   AliMUONSegment *segment;
472   Bool_t last2st;
473   Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
474       impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
475   AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",Station));
476   // first and second chambers (0...) in the station
477   Int_t ch1 = 2 * Station;
478   Int_t ch2 = ch1 + 1;
479   // variable true for stations downstream of the dipole:
480   // Station(0..4) equal to 3 or 4
481   if ((Station == 3) || (Station == 4)) {
482     last2st = kTRUE;
483     // maximum impact parameter (cm) according to fMinBendingMomentum...
484     maxImpactParam =
485       TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
486     // minimum impact parameter (cm) according to fMaxBendingMomentum...
487     minImpactParam =
488       TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum));
489   }
490   else last2st = kFALSE;
491   // extrapolation factor from Z of first chamber to Z of second chamber
492   // dZ to be changed to take into account fine structure of chambers ????
493   Double_t extrapFact;
494   // index for current segment
495   Int_t segmentIndex = 0;
496   // Loop over HitsForRec in the first chamber of the station
497   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
498        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
499        hit1++) {
500     // pointer to the HitForRec
501     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
502     // extrapolation,
503     // on the straight line joining the HitForRec to the vertex (0,0,0),
504     // to the Z of the second chamber of the station
505     // Loop over HitsForRec in the second chamber of the station
506     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
507          hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
508          hit2++) {
509       // pointer to the HitForRec
510       hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
511       // absolute values of distances, in bending and non bending planes,
512       // between the HitForRec in the second chamber
513       // and the previous extrapolation
514       extrapFact = hit2Ptr->GetZ()/ hit1Ptr->GetZ();
515       extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
516       extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
517       distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
518       distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
519       if (last2st) {
520         // bending slope
521         if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0.0 ) {
522           bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
523             (hit1Ptr->GetZ() - hit2Ptr->GetZ());
524           // absolute value of impact parameter
525           impactParam =
526             TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope);
527          } 
528          else {
529            AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): impactParam set to maxImpactParam");
530            impactParam = maxImpactParam;   
531          }   
532       }
533       // check for distances not too large,
534       // and impact parameter not too big if stations downstream of the dipole.
535       // Conditions "distBend" and "impactParam" correlated for these stations ????
536       if ((distBend < fSegmentMaxDistBending[Station]) &&
537           (distNonBend < fSegmentMaxDistNonBending[Station]) &&
538           (!last2st || (impactParam < maxImpactParam)) &&
539           (!last2st || (impactParam > minImpactParam))) {
540         // make new segment
541         segment = new ((*fSegmentsPtr[Station])[segmentIndex])
542           AliMUONSegment(hit1Ptr, hit2Ptr);
543         // update "link" to this segment from the hit in the first chamber
544         if (hit1Ptr->GetNSegments() == 0)
545           hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
546         hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
547         if (AliLog::GetGlobalDebugLevel() > 1) {
548           cout << "segmentIndex(0...): " << segmentIndex
549                << "  distBend: " << distBend
550                << "  distNonBend: " << distNonBend
551                << endl;
552           segment->Dump();
553           cout << "HitForRec in first chamber" << endl;
554           hit1Ptr->Dump();
555           cout << "HitForRec in second chamber" << endl;
556           hit2Ptr->Dump();
557         };
558         // increment index for current segment
559         segmentIndex++;
560       }
561     } //for (Int_t hit2
562   } // for (Int_t hit1...
563   fNSegments[Station] = segmentIndex;
564   // Sorting according to "impact parameter" if station(1..5) 4 or 5,
565   // i.e. Station(0..4) 3 or 4, using the function "Compare".
566   // After this sorting, it is impossible to use
567   // the "fNSegments" and "fIndexOfFirstSegment"
568   // of the HitForRec in the first chamber to explore all segments formed with it.
569   // Is this sorting really needed ????
570   if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
571   AliDebug(1,Form("Station: %d  NSegments:  %d ", Station, fNSegments[Station]));
572   return;
573 }
574
575   //__________________________________________________________________________
576 void AliMUONTrackReconstructor::MakeTracks(void)
577 {
578   // To make the tracks,
579   // from the list of segments and points in all stations
580   AliDebug(1,"Enter MakeTracks");
581   // The order may be important for the following Reset's
582   //AZ ResetTracks();
583   //AZ ResetTrackHits();
584   if (fTrackMethod != 1) { //AZ - Kalman filter
585     MakeTrackCandidatesK();
586     if (fRecTracksPtr->GetEntriesFast() == 0) return;
587     // Follow tracks in stations(1..) 3, 2 and 1
588     FollowTracksK();
589     // Remove double tracks
590     RemoveDoubleTracksK();
591     // Propagate tracks to the vertex thru absorber
592     GoToVertex();
593     // Fill AliMUONTrack data members
594     FillMUONTrack();
595   } else { 
596     // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
597     MakeTrackCandidates();
598     // Follow tracks in stations(1..) 3, 2 and 1
599     FollowTracks();
600     // Remove double tracks
601     RemoveDoubleTracks();
602     UpdateTrackParamAtHit();
603     UpdateHitForRecAtHit();
604   }
605   return;
606 }
607
608   //__________________________________________________________________________
609 void AliMUONTrackReconstructor::ValidateTracksWithTrigger(void)
610 {
611   // Try to match track from tracking system with trigger track
612   AliMUONTrack *track;
613   TClonesArray *recTriggerTracks;
614   
615   fMUONData->SetTreeAddress("RL");
616   fMUONData->GetRecTriggerTracks();
617   recTriggerTracks = fMUONData->RecTriggerTracks();
618
619   track = (AliMUONTrack*) fRecTracksPtr->First();
620   while (track) {
621     track->MatchTriggerTrack(recTriggerTracks);
622     track = (AliMUONTrack*) fRecTracksPtr->After(track);
623   }
624
625 }
626
627   //__________________________________________________________________________
628 Bool_t AliMUONTrackReconstructor::MakeTriggerTracks(void)
629 {
630     // To make the trigger tracks from Local Trigger
631   AliDebug(1, "Enter MakeTriggerTracks");
632     
633     Int_t nTRentries;
634     Long_t gloTrigPat;
635     TClonesArray *localTrigger;
636     TClonesArray *globalTrigger;
637     AliMUONLocalTrigger *locTrg;
638     AliMUONGlobalTrigger *gloTrg;
639
640     TTree* treeR = fLoader->TreeR();
641     
642     // Loading MUON subsystem
643     AliMUON * pMUON = (AliMUON *) gAlice->GetDetector("MUON");
644
645     nTRentries = Int_t(treeR->GetEntries());
646      
647     treeR->GetEvent(0); // only one entry  
648
649     if (!(fMUONData->IsTriggerBranchesInTree())) {
650       AliWarning(Form("Trigger information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
651       return kFALSE;
652     }
653
654     fMUONData->SetTreeAddress("TC");
655     fMUONData->GetTrigger();
656
657     // global trigger for trigger pattern
658     gloTrigPat = 0;
659     globalTrigger = fMUONData->GlobalTrigger(); 
660     gloTrg = (AliMUONGlobalTrigger*)globalTrigger->UncheckedAt(0);
661  
662     if (gloTrg)
663       gloTrigPat = gloTrg->GetGlobalPattern();
664   
665
666     // local trigger for tracking 
667     localTrigger = fMUONData->LocalTrigger();    
668     Int_t nlocals = (Int_t) (localTrigger->GetEntries());
669
670     Float_t z11 = AliMUONConstants::DefaultChamberZ(10);
671     Float_t z21 = AliMUONConstants::DefaultChamberZ(12);
672
673     Float_t y11 = 0.;
674     Int_t stripX21 = 0;
675     Float_t y21 = 0.;
676     Float_t x11 = 0.;
677
678     for (Int_t i=0; i<nlocals; i++) { // loop on Local Trigger
679       locTrg = (AliMUONLocalTrigger*)localTrigger->UncheckedAt(i);      
680
681       AliDebug(1, "AliMUONTrackReconstructor::MakeTriggerTrack using NEW trigger \n");
682       AliMUONTriggerCircuitNew * circuit = 
683           &(pMUON->TriggerCircuitNew(locTrg->LoCircuit()-1)); // -1 !!!
684       y11 = circuit->GetY11Pos(locTrg->LoStripX()); 
685       stripX21 = locTrg->LoStripX()+locTrg->LoDev()+1;
686       y21 = circuit->GetY21Pos(stripX21);       
687       x11 = circuit->GetX11Pos(locTrg->LoStripY());
688       
689       AliDebug(1, Form(" MakeTriggerTrack %d %d %d %d %d %f %f %f \n",i,locTrg->LoCircuit(),
690                        locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11));
691       
692       Float_t thetax = TMath::ATan2( x11 , z11 );
693       Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) );
694       
695       fTriggerTrack->SetX11(x11);
696       fTriggerTrack->SetY11(y11);
697       fTriggerTrack->SetThetax(thetax);
698       fTriggerTrack->SetThetay(thetay);
699       fTriggerTrack->SetGTPattern(gloTrigPat);
700             
701       fMUONData->AddRecTriggerTrack(*fTriggerTrack);
702     } // end of loop on Local Trigger
703     return kTRUE;    
704 }
705
706   //__________________________________________________________________________
707 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
708 {
709   // To make track candidates with two segments in stations(1..) 4 and 5,
710   // the first segment being pointed to by "BegSegment".
711   // Returns the number of such track candidates.
712   Int_t endStation, iEndSegment, nbCan2Seg;
713   AliMUONSegment *endSegment;
714   AliMUONSegment *extrapSegment = NULL;
715   AliMUONTrack *recTrack;
716   Double_t mcsFactor;
717   AliDebug(1,"Enter MakeTrackCandidatesWithTwoSegments");
718   // Station for the end segment
719   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
720   // multiple scattering factor corresponding to one chamber
721   mcsFactor = 0.0136 /
722     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
723   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
724   // linear extrapolation to end station
725   // number of candidates with 2 segments to 0
726   nbCan2Seg = 0;
727   // Loop over segments in the end station
728   for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
729     // pointer to segment
730     endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
731     // test compatibility between current segment and "extrapSegment"
732     // 4 because 4 quantities in chi2
733     extrapSegment =
734       BegSegment->CreateSegmentFromLinearExtrapToStation(endSegment->GetZ(), mcsFactor);
735     if ((endSegment->
736          NormalizedChi2WithSegment(extrapSegment,
737                                    fMaxSigma2Distance)) <= 4.0) {
738       // both segments compatible:
739       // make track candidate from "begSegment" and "endSegment"
740       AliDebug(2,Form("TrackCandidate with Segment %d in Station(0..) %d", iEndSegment, endStation));
741       // flag for both segments in one track:
742       // to be done in track constructor ????
743       BegSegment->SetInTrack(kTRUE);
744       endSegment->SetInTrack(kTRUE);
745       recTrack = new ((*fRecTracksPtr)[fNRecTracks])
746         AliMUONTrack(BegSegment, endSegment, this);
747       fNRecTracks++;
748       if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
749       // increment number of track candidates with 2 segments
750       nbCan2Seg++;
751     }
752     delete extrapSegment; // should not delete HitForRec's it points to !!!!
753   } // for (iEndSegment = 0;...
754   return nbCan2Seg;
755 }
756
757   //__________________________________________________________________________
758 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
759 {
760   // To make track candidates with one segment and one point
761   // in stations(1..) 4 and 5,
762   // the segment being pointed to by "BegSegment".
763   Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
764   AliMUONHitForRec *extrapHitForRec= NULL;
765   AliMUONHitForRec *hit;
766   AliMUONTrack *recTrack;
767   Double_t mcsFactor;
768   AliDebug(1,"Enter MakeTrackCandidatesWithOneSegmentAndOnePoint");
769   // station for the end point
770   endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
771   // multiple scattering factor corresponding to one chamber
772   mcsFactor = 0.0136 /
773     GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
774   mcsFactor     = fChamberThicknessInX0 * mcsFactor * mcsFactor;
775   // first and second chambers(0..) in the end station
776   ch1 = 2 * endStation;
777   ch2 = ch1 + 1;
778   // number of candidates to 0
779   nbCan1Seg1Hit = 0;
780   // Loop over chambers of the end station
781   for (ch = ch2; ch >= ch1; ch--) {
782     // limits for the hit index in the loop
783     iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
784     iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
785     // Loop over HitForRec's in the chamber
786     for (iHit = iHitMin; iHit < iHitMax; iHit++) {
787       // pointer to HitForRec
788       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
789       // test compatibility between current HitForRec and "extrapHitForRec"
790       // 2 because 2 quantities in chi2
791       // linear extrapolation to chamber
792       extrapHitForRec =
793         BegSegment->CreateHitForRecFromLinearExtrapToChamber( hit->GetZ(), mcsFactor);
794       if ((hit->
795            NormalizedChi2WithHitForRec(extrapHitForRec,
796                                        fMaxSigma2Distance)) <= 2.0) {
797         // both HitForRec's compatible:
798         // make track candidate from begSegment and current HitForRec
799         AliDebug(2, Form("TrackCandidate with HitForRec  %d in Chamber(0..) %d", iHit, ch));
800         // flag for beginning segments in one track:
801         // to be done in track constructor ????
802         BegSegment->SetInTrack(kTRUE);
803         recTrack = new ((*fRecTracksPtr)[fNRecTracks])
804           AliMUONTrack(BegSegment, hit, this);
805         // the right place to eliminate "double counting" ???? how ????
806         fNRecTracks++;
807         if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
808         // increment number of track candidates
809         nbCan1Seg1Hit++;
810       }
811       delete extrapHitForRec;
812     } // for (iHit = iHitMin;...
813   } // for (ch = ch2;...
814   return nbCan1Seg1Hit;
815 }
816
817   //__________________________________________________________________________
818 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
819 {
820   // To make track candidates
821   // with at least 3 aligned points in stations(1..) 4 and 5
822   // (two Segment's or one Segment and one HitForRec)
823   Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
824   AliMUONSegment *begSegment;
825   AliDebug(1,"Enter MakeTrackCandidates");
826   // Loop over stations(1..) 5 and 4 for the beginning segment
827   for (begStation = 4; begStation > 2; begStation--) {
828     // Loop over segments in the beginning station
829     for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
830       // pointer to segment
831       begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
832       AliDebug(2,Form("Look for TrackCandidate's with Segment %d  in Station(0..) %d", iBegSegment, begStation));
833       // Look for track candidates with two segments,
834       // "begSegment" and all compatible segments in other station.
835       // Only for beginning station(1..) 5
836       // because candidates with 2 segments have to looked for only once.
837       if (begStation == 4)
838         nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
839       // Look for track candidates with one segment and one point,
840       // "begSegment" and all compatible HitForRec's in other station.
841       // Only if "begSegment" does not belong already to a track candidate.
842       // Is that a too strong condition ????
843       if (!(begSegment->GetInTrack()))
844         nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
845     } // for (iBegSegment = 0;...
846   } // for (begStation = 4;...
847   return;
848 }
849
850   //__________________________________________________________________________
851 void AliMUONTrackReconstructor::FollowTracks(void)
852 {
853   // Follow tracks in stations(1..) 3, 2 and 1
854   // too long: should be made more modular !!!!
855   AliMUONHitForRec *bestHit, *extrapHit, *hit;
856   AliMUONSegment *bestSegment, *extrapSegment, *segment;
857   AliMUONTrack *track, *nextTrack;
858   AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
859   // -1 to avoid compilation warnings
860   Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex; 
861   Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
862   Double_t bendingMomentum, chi2Norm = 0.;
863
864
865   // local maxSigma2Distance, for easy increase in testing
866   maxSigma2Distance = fMaxSigma2Distance;
867   AliDebug(2,"Enter FollowTracks");
868   // Loop over track candidates
869   track = (AliMUONTrack*) fRecTracksPtr->First();
870   trackIndex = -1;
871   while (track) {
872     // Follow function for each track candidate ????
873     trackIndex++;
874     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
875     AliDebug(2,Form("FollowTracks: track candidate(0..): %d", trackIndex));
876     // Fit track candidate
877     track->SetFitMCS(0); // without multiple Coulomb scattering
878     track->SetFitNParam(3); // with 3 parameters (X = Y = 0)
879     track->SetFitStart(0); // from parameters at vertex
880     track->Fit();
881     if (AliLog::GetGlobalDebugLevel()> 2) {
882       cout << "FollowTracks: track candidate(0..): " << trackIndex
883            << " after fit in stations(0..) 3 and 4" << endl;
884       track->RecursiveDump();
885     }
886     // Loop over stations(1..) 3, 2 and 1
887     // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
888     // otherwise: majority coincidence 2 !!!!
889     for (station = 2; station >= 0; station--) {
890       // Track parameters at first track hit (smallest Z)
891       trackParam1 = ((AliMUONTrackHit*)
892                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
893       // extrapolation to station
894       trackParam1->ExtrapToStation(station, trackParam);
895       extrapSegment = new AliMUONSegment(); //  empty segment
896       // multiple scattering factor corresponding to one chamber
897       // and momentum in bending plane (not total)
898       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
899       mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
900       // Z difference from previous station
901       dZ1 = AliMUONConstants::DefaultChamberZ(2 * station) -
902             AliMUONConstants::DefaultChamberZ(2 * station + 2);
903       // Z difference between the two previous stations
904       dZ2 = AliMUONConstants::DefaultChamberZ(2 * station + 2) -
905             AliMUONConstants::DefaultChamberZ(2 * station + 4);
906       // Z difference between the two chambers in the previous station
907       dZ3 = AliMUONConstants::DefaultChamberZ(2 * station) -
908             AliMUONConstants::DefaultChamberZ(2 * station + 1);
909       extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
910       extrapSegment->
911         SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
912       extrapSegment->UpdateFromStationTrackParam
913         (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
914          trackParam1->GetInverseBendingMomentum());
915       bestChi2 = 5.0;
916       bestSegment = NULL;
917       if (AliLog::GetGlobalDebugLevel() > 2) {
918         cout << "FollowTracks: track candidate(0..): " << trackIndex
919              << " Look for segment in station(0..): " << station << endl;
920       }
921
922       // Loop over segments in station
923       for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
924         // Look for best compatible Segment in station
925         // should consider all possibilities ????
926         // multiple scattering ????
927         // separation in 2 functions: Segment and HitForRec ????
928         segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
929         // correction of corrected segment (fBendingCoor and fNonBendingCoor)
930         // according to real Z value of "segment" and slopes of "extrapSegment"
931         (&(trackParam[0]))->ExtrapToZ(segment->GetZ());
932         (&(trackParam[1]))->ExtrapToZ(segment->GetZ());
933         extrapSegment->SetBendingCoor((&(trackParam[0]))->GetBendingCoor());
934         extrapSegment->SetNonBendingCoor((&(trackParam[0]))->GetNonBendingCoor());
935         extrapSegment->SetBendingSlope((&(trackParam[0]))->GetBendingSlope());
936         extrapSegment->SetNonBendingSlope((&(trackParam[0]))->GetNonBendingSlope());
937         chi2 = segment->
938           NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
939         if (chi2 < bestChi2) {
940           // update best Chi2 and Segment if better found
941           bestSegment = segment;
942           bestChi2 = chi2;
943         }
944       }
945       if (bestSegment) {
946         // best segment found: add it to track candidate
947         (&(trackParam[0]))->ExtrapToZ(bestSegment->GetZ());
948         (&(trackParam[1]))->ExtrapToZ(bestSegment->GetZ());
949         track->AddSegment(bestSegment);
950         // set track parameters at these two TrakHit's
951         track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0]));
952         track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1]));
953         AliDebug(3, Form("FollowTracks: track candidate(0..): %d  Added segment in station(0..): %d", trackIndex, station));
954         if (AliLog::GetGlobalDebugLevel()>2) track->RecursiveDump();
955       }
956       else {
957         // No best segment found:
958         // Look for best compatible HitForRec in station:
959         // should consider all possibilities ????
960         // multiple scattering ???? do about like for extrapSegment !!!!
961         extrapHit = new AliMUONHitForRec(); //  empty hit
962         bestChi2 = 3.0;
963         bestHit = NULL;
964         AliDebug(3, Form("FollowTracks: track candidate(0..): %d Segment not found, look for hit in station(0..): %d ", 
965                          trackIndex, station));
966         
967         // Loop over chambers of the station
968         for (chInStation = 0; chInStation < 2; chInStation++) {
969           ch = 2 * station + chInStation;
970           for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
971                iHit < fIndexOfFirstHitForRecPerChamber[ch] +
972                  fNHitsForRecPerChamber[ch];
973                iHit++) {
974             hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
975             // coordinates of extrapolated hit
976             (&(trackParam[chInStation]))->ExtrapToZ(hit->GetZ());
977             extrapHit->
978               SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
979             extrapHit->
980               SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
981             // resolutions from "extrapSegment"
982             extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
983             extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
984             // Loop over hits in the chamber
985             // condition for hit not already in segment ????
986             chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
987             if (chi2 < bestChi2) {
988               // update best Chi2 and HitForRec if better found
989               bestHit = hit;
990               bestChi2 = chi2;
991               chBestHit = chInStation;
992             }
993           }
994         }
995         if (bestHit) {
996           // best hit found: add it to track candidate
997           (&(trackParam[chBestHit]))->ExtrapToZ(bestHit->GetZ());
998           track->AddHitForRec(bestHit);
999           // set track parameters at this TrackHit
1000           track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
1001                                     &(trackParam[chBestHit]));
1002           if (AliLog::GetGlobalDebugLevel() > 2) {
1003             cout << "FollowTracks: track candidate(0..): " << trackIndex
1004                  << " Added hit in station(0..): " << station << endl;
1005             track->RecursiveDump();
1006           }
1007         }
1008         else {
1009           // Remove current track candidate
1010           // and corresponding TrackHit's, ...
1011           track->Remove();
1012           delete extrapSegment;
1013           delete extrapHit;
1014           break; // stop the search for this candidate:
1015           // exit from the loop over station
1016         }
1017         delete extrapHit;
1018       }
1019       delete extrapSegment;
1020       // Sort track hits according to increasing Z
1021       track->GetTrackHitsPtr()->Sort();
1022       // Update track parameters at first track hit (smallest Z)
1023       trackParam1 = ((AliMUONTrackHit*)
1024                      (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1025       bendingMomentum = 0.;
1026       if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
1027         bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
1028       // Track removed if bendingMomentum not in window [min, max]
1029       if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
1030         track->Remove();
1031         break; // stop the search for this candidate:
1032         // exit from the loop over station 
1033       }
1034       // Track fit
1035       // with multiple Coulomb scattering if all stations
1036       if (station == 0) track->SetFitMCS(1);
1037       // without multiple Coulomb scattering if not all stations
1038       else track->SetFitMCS(0);
1039       track->SetFitNParam(5);  // with 5 parameters (momentum and position)
1040       track->SetFitStart(1);  // from parameters at first hit
1041       track->Fit();
1042       Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
1043       if (numberOfDegFree > 0) {
1044         chi2Norm =  track->GetFitFMin() / numberOfDegFree;
1045       } else {
1046         chi2Norm = 1.e10;
1047       }
1048       // Track removed if normalized chi2 too high
1049       if (chi2Norm > fMaxChi2) {
1050         track->Remove();
1051         break; // stop the search for this candidate:
1052         // exit from the loop over station 
1053       }
1054       if (AliLog::GetGlobalDebugLevel() > 2) {
1055         cout << "FollowTracks: track candidate(0..): " << trackIndex
1056              << " after fit from station(0..): " << station << " to 4" << endl;
1057         track->RecursiveDump();
1058       }
1059       // Track extrapolation to the vertex through the absorber (Branson)
1060       // after going through the first station
1061       if (station == 0) {
1062         trackParamVertex = *trackParam1;
1063         (&trackParamVertex)->ExtrapToVertex(0.,0.,0.);
1064         track->SetTrackParamAtVertex(&trackParamVertex);
1065         if (AliLog::GetGlobalDebugLevel() > 0) {
1066           cout << "FollowTracks: track candidate(0..): " << trackIndex
1067                << " after extrapolation to vertex" << endl;
1068           track->RecursiveDump();
1069         }
1070       }
1071     } // for (station = 2;...
1072     // go really to next track
1073     track = nextTrack;
1074   } // while (track)
1075   // Compression of track array (necessary after Remove ????)
1076   fRecTracksPtr->Compress();
1077   return;
1078 }
1079
1080   //__________________________________________________________________________
1081 void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
1082 {
1083   // To remove double tracks.
1084   // Tracks are considered identical
1085   // if they have at least half of their hits in common.
1086   // Among two identical tracks, one keeps the track with the larger number of hits
1087   // or, if these numbers are equal, the track with the minimum Chi2.
1088   AliMUONTrack *track1, *track2, *trackToRemove;
1089   Bool_t identicalTracks;
1090   Int_t hitsInCommon, nHits1, nHits2;
1091   identicalTracks = kTRUE;
1092   while (identicalTracks) {
1093     identicalTracks = kFALSE;
1094     // Loop over first track of the pair
1095     track1 = (AliMUONTrack*) fRecTracksPtr->First();
1096     while (track1 && (!identicalTracks)) {
1097       nHits1 = track1->GetNTrackHits();
1098       // Loop over second track of the pair
1099       track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1100       while (track2 && (!identicalTracks)) {
1101         nHits2 = track2->GetNTrackHits();
1102         // number of hits in common between two tracks
1103         hitsInCommon = track1->HitsInCommon(track2);
1104         // check for identical tracks
1105         if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
1106           identicalTracks = kTRUE;
1107           // decide which track to remove
1108           if (nHits1 > nHits2) trackToRemove = track2;
1109           else if (nHits1 < nHits2) trackToRemove = track1;
1110           else if ((track1->GetFitFMin()) < (track2->GetFitFMin()))
1111             trackToRemove = track2;
1112           else trackToRemove = track1;
1113           // remove it
1114           trackToRemove->Remove();
1115         }
1116         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1117       } // track2
1118       track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1119     } // track1
1120   }
1121   return;
1122 }
1123
1124   //__________________________________________________________________________
1125 void AliMUONTrackReconstructor::UpdateTrackParamAtHit()
1126 {
1127   // Set track parameters after track fitting. Fill fTrackParamAtHit of AliMUONTrack's
1128   AliMUONTrack *track;
1129   AliMUONTrackHit *trackHit;
1130   AliMUONTrackParam *trackParam;
1131   track = (AliMUONTrack*) fRecTracksPtr->First();
1132   while (track) {
1133     trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->First();
1134     while (trackHit) {
1135       trackParam = trackHit->GetTrackParam();
1136       track->AddTrackParamAtHit(trackParam);
1137       trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->After(trackHit); 
1138     } // trackHit    
1139     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1140   } // track
1141   return;
1142 }
1143
1144   //__________________________________________________________________________
1145 void AliMUONTrackReconstructor::UpdateHitForRecAtHit()
1146 {
1147   // Set cluster parameterss after track fitting. Fill fHitForRecAtHit of AliMUONTrack's
1148   AliMUONTrack *track;
1149   AliMUONTrackHit *trackHit;
1150   AliMUONHitForRec *hitForRec;
1151   track = (AliMUONTrack*) fRecTracksPtr->First();
1152   while (track) {
1153     trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->First();
1154     while (trackHit) {
1155       hitForRec = trackHit->GetHitForRecPtr();
1156       track->AddHitForRecAtHit(hitForRec);
1157       trackHit = (AliMUONTrackHit*) (track->GetTrackHitsPtr())->After(trackHit); 
1158     } // trackHit    
1159     track = (AliMUONTrack*) fRecTracksPtr->After(track);
1160   } // track
1161   return;
1162 }
1163
1164   //__________________________________________________________________________
1165 void AliMUONTrackReconstructor::FillMUONTrack()
1166 {
1167   // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's
1168   AliMUONTrackK *track;
1169   track = (AliMUONTrackK*) fRecTracksPtr->First();
1170   while (track) {
1171     track->FillMUONTrack();
1172     track = (AliMUONTrackK*) fRecTracksPtr->After(track);
1173   } 
1174   return;
1175 }
1176
1177   //__________________________________________________________________________
1178 void AliMUONTrackReconstructor::EventDump(void)
1179 {
1180   // Dump reconstructed event (track parameters at vertex and at first hit),
1181   // and the particle parameters
1182
1183   AliMUONTrack *track;
1184   AliMUONTrackParam *trackParam, *trackParam1;
1185   Double_t bendingSlope, nonBendingSlope, pYZ;
1186   Double_t pX, pY, pZ, x, y, z, c;
1187   Int_t np, trackIndex, nTrackHits;
1188  
1189   AliDebug(1,"****** enter EventDump ******");
1190   AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); 
1191   
1192   fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1193   // Loop over reconstructed tracks
1194   for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1195     if (fTrackMethod != 1) continue; //AZ - skip the rest for now
1196     AliDebug(1, Form("track number: %d", trackIndex));
1197     // function for each track for modularity ????
1198     track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1199     nTrackHits = track->GetNTrackHits();
1200     AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1201     // track parameters at Vertex
1202     trackParam = track->GetTrackParamAtVertex();
1203     x = trackParam->GetNonBendingCoor();
1204     y = trackParam->GetBendingCoor();
1205     z = trackParam->GetZ();
1206     bendingSlope = trackParam->GetBendingSlope();
1207     nonBendingSlope = trackParam->GetNonBendingSlope();
1208     pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1209     pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1210     pX = pZ * nonBendingSlope;
1211     pY = pZ * bendingSlope;
1212     c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1213     AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1214                      z, x, y, pX, pY, pZ, c));
1215
1216     // track parameters at first hit
1217     trackParam1 = ((AliMUONTrackHit*)
1218                    (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1219     x = trackParam1->GetNonBendingCoor();
1220     y = trackParam1->GetBendingCoor();
1221     z = trackParam1->GetZ();
1222     bendingSlope = trackParam1->GetBendingSlope();
1223     nonBendingSlope = trackParam1->GetNonBendingSlope();
1224     pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1225     pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1226     pX = pZ * nonBendingSlope;
1227     pY = pZ * bendingSlope;
1228     c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1229     AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1230                      z, x, y, pX, pY, pZ, c));
1231   }
1232   // informations about generated particles NO !!!!!!!!
1233   
1234 //    for (Int_t iPart = 0; iPart < np; iPart++) {
1235 //      p = gAlice->Particle(iPart);
1236 //      printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1237 //         iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());    
1238 //    }
1239   return;
1240 }
1241
1242
1243 //__________________________________________________________________________
1244 void AliMUONTrackReconstructor::EventDumpTrigger(void)
1245 {
1246   // Dump reconstructed trigger event 
1247   // and the particle parameters
1248     
1249   AliMUONTriggerTrack *triggertrack ;
1250   Int_t nTriggerTracks = fMUONData->RecTriggerTracks()->GetEntriesFast();
1251  
1252   AliDebug(1, "****** enter EventDumpTrigger ******");
1253   AliDebug(1, Form("Number of Reconstructed tracks : %d ",  nTriggerTracks));
1254   
1255   // Loop over reconstructed tracks
1256   for (Int_t trackIndex = 0; trackIndex < nTriggerTracks; trackIndex++) {
1257     triggertrack = (AliMUONTriggerTrack*)fMUONData->RecTriggerTracks()->At(trackIndex);
1258       printf(" trigger track number %i x11=%f y11=%f thetax=%f thetay=%f \n",
1259              trackIndex,
1260              triggertrack->GetX11(),triggertrack->GetY11(),
1261              triggertrack->GetThetax(),triggertrack->GetThetay());      
1262   } 
1263 }
1264
1265 //__________________________________________________________________________
1266 void AliMUONTrackReconstructor::MakeTrackCandidatesK(void)
1267 {
1268   // To make initial tracks for Kalman filter from the list of segments
1269   Int_t istat, iseg;
1270   AliMUONSegment *segment;
1271   AliMUONTrackK *trackK;
1272
1273   AliDebug(1,"Enter MakeTrackCandidatesK");
1274
1275   AliMUONTrackK a(this, fHitsForRecPtr);
1276   // Loop over stations(1...) 5 and 4
1277   for (istat=4; istat>=3; istat--) {
1278     // Loop over segments in the station
1279     for (iseg=0; iseg<fNSegments[istat]; iseg++) {
1280       // Transform segments to tracks and evaluate covariance matrix
1281       segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
1282       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
1283     } // for (iseg=0;...)
1284   } // for (istat=4;...)
1285   return;
1286 }
1287
1288 //__________________________________________________________________________
1289 void AliMUONTrackReconstructor::FollowTracksK(void)
1290 {
1291   // Follow tracks using Kalman filter
1292   Bool_t ok;
1293   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
1294   Double_t zDipole1, zDipole2;
1295   AliMUONTrackK *trackK;
1296   AliMUONHitForRec *hit;
1297   AliMUONRawCluster *clus;
1298   TClonesArray *rawclusters;
1299   clus = 0; rawclusters = 0;
1300
1301   zDipole1 = GetSimpleBPosition() + GetSimpleBLength()/2;
1302   zDipole2 = zDipole1 - GetSimpleBLength();
1303
1304   // Print hits
1305   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
1306
1307   if (trackK->DebugLevel() > 0) {
1308     for (Int_t i1=0; i1<fNHitsForRec; i1++) {
1309       hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
1310       printf(" Hit # %d %10.4f %10.4f %10.4f",
1311              hit->GetChamberNumber(), hit->GetBendingCoor(),
1312              hit->GetNonBendingCoor(), hit->GetZ());
1313  
1314       // from raw clusters
1315       rawclusters = fMUONData->RawClusters(hit->GetChamberNumber());
1316       clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
1317                                                            GetHitNumber());
1318       printf(" %d", clus->GetTrack(1));
1319       if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2));
1320       else printf("\n");
1321      
1322     }
1323   } // if (trackK->DebugLevel() > 0)
1324
1325   icand = -1;
1326   Int_t nSeeds;
1327   nSeeds = fNRecTracks; // starting number of seeds
1328   // Loop over track candidates
1329   while (icand < fNRecTracks-1) {
1330     icand ++;
1331     if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl;
1332     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1333     if (trackK->GetRecover() < 0) continue; // failed track
1334
1335     // Discard candidate which will produce the double track
1336     /*
1337     if (icand > 0) {
1338       ok = CheckCandidateK(icand,nSeeds);
1339       if (!ok) {
1340         trackK->SetRecover(-1); // mark candidate to be removed
1341         continue;
1342       }
1343     }
1344     */
1345
1346     ok = kTRUE;
1347     if (trackK->GetRecover() == 0) 
1348       hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit
1349     else 
1350       hit = trackK->GetHitLastOk(); // hit where track stopped
1351
1352     if (hit) ichamBeg = hit->GetChamberNumber();
1353     ichamEnd = 0;
1354     // Check propagation direction
1355     if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); }
1356     else if (trackK->GetTrackDir() < 0) {
1357       ichamEnd = 9; // forward propagation
1358       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1359       if (ok) {
1360         ichamBeg = ichamEnd;
1361         ichamEnd = 6; // backward propagation
1362         // Change weight matrix and zero fChi2 for backpropagation
1363         trackK->StartBack();
1364         trackK->SetTrackDir(1);
1365         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1366         ichamBeg = ichamEnd;
1367         ichamEnd = 0;
1368       }
1369     } else {
1370       if (trackK->GetBPFlag()) {
1371         // backpropagation
1372         ichamEnd = 6; // backward propagation
1373         // Change weight matrix and zero fChi2 for backpropagation
1374         trackK->StartBack();
1375         ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1376         ichamBeg = ichamEnd;
1377         ichamEnd = 0;
1378       }
1379     }
1380
1381     if (ok) {
1382       trackK->SetTrackDir(1);
1383       trackK->SetBPFlag(kFALSE);
1384       ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1385     }
1386     if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed
1387
1388     // Apply smoother
1389     if (trackK->GetRecover() >= 0) {
1390       ok = trackK->Smooth();
1391       if (!ok) trackK->SetRecover(-1); // mark candidate to be removed
1392     }
1393
1394     // Majority 3 of 4 in first 2 stations
1395     if (!ok) continue;
1396     chamBits = 0;
1397     Double_t chi2max = 0;
1398     for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
1399       hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i];
1400       chamBits |= BIT(hit->GetChamberNumber());
1401       if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i);
1402     }
1403     if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) {
1404       //trackK->Recover();
1405       trackK->SetRecover(-1); //mark candidate to be removed
1406       continue;
1407     }
1408     if (ok) trackK->SetTrackQuality(0); // compute "track quality"
1409   } // while
1410
1411   for (Int_t i=0; i<fNRecTracks; i++) {
1412     trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1413     if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
1414   }
1415
1416   // Compress TClonesArray
1417   fRecTracksPtr->Compress();
1418   fNRecTracks = fRecTracksPtr->GetEntriesFast();
1419   return;
1420 }
1421
1422 //__________________________________________________________________________
1423 Bool_t AliMUONTrackReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds) const
1424 {
1425   // Discards track candidate if it will produce the double track (having
1426   // the same seed segment hits as hits of a good track found before)
1427   AliMUONTrackK *track1, *track2;
1428   AliMUONHitForRec *hit1, *hit2, *hit;
1429
1430   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1431   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
1432   hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit
1433
1434   for (Int_t i=0; i<icand; i++) {
1435     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1436     //if (track2->GetRecover() < 0) continue;
1437     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
1438
1439     if (track1->GetStartSegment() == track2->GetStartSegment()) {
1440       return kFALSE;
1441     } else {
1442       Int_t nSame = 0;
1443       for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
1444         hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j];
1445         if (hit == hit1 || hit == hit2) {
1446           nSame++;
1447           if (nSame == 2) return kFALSE;
1448         }
1449       } // for (Int_t j=0;
1450     }
1451   } // for (Int_t i=0;
1452   return kTRUE;
1453 }
1454
1455 //__________________________________________________________________________
1456 void AliMUONTrackReconstructor::RemoveDoubleTracksK(void)
1457 {
1458   // Removes double tracks (sharing more than half of their hits). Keeps
1459   // the track with higher quality
1460   AliMUONTrackK *track1, *track2, *trackToKill;
1461
1462   // Sort tracks according to their quality
1463   fRecTracksPtr->Sort();
1464
1465   // Loop over first track of the pair
1466   track1 = (AliMUONTrackK*) fRecTracksPtr->First();
1467   Int_t debug = track1->DebugLevel();
1468   while (track1) {
1469     // Loop over second track of the pair
1470     track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1471     while (track2) {
1472       // Check whether or not to keep track2
1473       if (!track2->KeepTrack(track1)) {
1474         if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
1475           " " << track2->GetTrackQuality() << endl;
1476         trackToKill = track2;
1477         track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1478         trackToKill->Kill();
1479         fRecTracksPtr->Compress();
1480       } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1481     } // track2
1482     track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1483   } // track1
1484
1485   fNRecTracks = fRecTracksPtr->GetEntriesFast();
1486   if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl;
1487 }
1488
1489 //__________________________________________________________________________
1490 void AliMUONTrackReconstructor::GoToVertex(void)
1491 {
1492   // Propagates track to the vertex thru absorber
1493   // (using Branson correction for now)
1494
1495   Double_t zVertex;
1496   zVertex = 0;
1497   for (Int_t i=0; i<fNRecTracks; i++) {
1498     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
1499     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
1500     //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
1501     ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber
1502   }
1503 }
1504
1505 //__________________________________________________________________________
1506 void AliMUONTrackReconstructor::SetTrackMethod(Int_t iTrackMethod)
1507 {
1508   // Set track method and recreate track container if necessary
1509   
1510   fTrackMethod = TMath::Min (iTrackMethod, 3);
1511   fTrackMethod = TMath::Max (fTrackMethod, 1);
1512   if (fTrackMethod != 1) {
1513     if (fRecTracksPtr) delete fRecTracksPtr;
1514     fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10);
1515     if (fTrackMethod == 2) cout << " *** Tracking with the Kalman filter *** " << endl;
1516     else cout << " *** Combined cluster / track finder ***" << endl;
1517   } else cout << " *** Original tracking *** " << endl;
1518
1519 }