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