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