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