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