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