1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.31 2002/10/23 07:24:56 alibrary
19 Introducing Riostream.h
21 Revision 1.30 2002/10/14 14:57:29 hristov
22 Merging the VirtualMC branch to the main development branch (HEAD)
24 Revision 1.28.8.1 2002/10/11 06:56:47 hristov
25 Updating VirtualMC to v3-09-02
27 Revision 1.29 2002/09/20 13:32:26 cussonno
28 Minor bugs in the definition of the bending impact parameter corrected (thanks to A. Zinchenko)
30 Revision 1.28 2001/08/31 06:50:11 gosset
31 Different handling of TreeR for track reconstruction from raw clusters
33 Revision 1.27 2001/07/27 13:03:12 hristov
34 Default Branch split level set to 99
36 Revision 1.26 2001/05/03 08:11:31 hristov
37 stdlib.h included to define exit()
39 Revision 1.25 2001/04/25 14:50:42 gosset
40 Corrections to violations of coding conventions
42 Revision 1.24 2001/03/30 09:37:58 gosset
43 Initialisations of pointers... for GEANT background events in the constructor
45 Revision 1.23 2001/01/26 21:44:45 morsch
46 Use access functions to AliMUONDigit, ... member data.
48 Revision 1.22 2001/01/26 20:00:53 hristov
49 Major upgrade of AliRoot code
50 Revision 1.20 2001/01/08 11:01:02 gosset
51 Modifications used for addendum to Dimuon TDR (JP Cussonneau):
52 *. MaxBendingMomentum to make both a segment and a track (default 500)
53 *. MaxChi2 per degree of freedom to make a track (default 100)
54 *. MinBendingMomentum used also to make a track
55 and not only a segment (default 3)
56 *. wider roads for track search in stations 1 to 3
57 *. extrapolation to actual Z instead of Z(chamber) in FollowTracks
59 - limits on parameters X and Y (+/-500)
60 - covariance matrices in double precision
61 - normalization of covariance matrices before inversion
62 - suppression of Minuit printouts
63 *. correction against memory leak (delete extrapHit) in FollowTracks
64 *. RMax to 10 degrees with Z(chamber) instead of fixed values;
65 RMin and Rmax cuts suppressed in NewHitForRecFromGEANT,
66 because useless with realistic geometry
68 Revision 1.19 2000/11/20 11:24:10 gosset
69 New package for reconstructed tracks (A. Gheata):
70 * tree output in the file "tree_reco.root"
71 * display events and make histograms from this file
73 Revision 1.18 2000/10/26 12:47:03 gosset
74 Real distance between chambers of each station taken into account
75 for the reconstruction parameters "fSegmentMaxDistBending[5]"
77 Revision 1.17 2000/10/24 09:26:20 gosset
80 Revision 1.16 2000/10/24 09:22:35 gosset
81 Method AddHitsForRecFromRawClusters: real Z of raw cluster and not Z of chamber
83 Revision 1.15 2000/10/12 15:17:30 gosset
84 Sign of fSimpleBValue corrected: sign ox Bx and not Bz (thanks to Galina)
86 Revision 1.14 2000/10/04 18:21:26 morsch
89 Revision 1.13 2000/10/02 21:28:09 fca
90 Removal of useless dependecies via forward declarations
92 Revision 1.12 2000/10/02 16:58:29 egangler
93 Cleaning of the code :
96 -> some useless includes removed or replaced by "class" statement
98 Revision 1.11 2000/09/22 09:16:33 hristov
101 Revision 1.10 2000/09/19 09:49:50 gosset
102 AliMUONEventReconstructor package
103 * track extrapolation independent from reco_muon.F, use of AliMagF...
104 * possibility to use new magnetic field (automatic from generated root file)
106 Revision 1.9 2000/07/20 12:45:27 gosset
107 New "EventReconstructor..." structure,
108 hopefully more adapted to tree/streamer.
109 "AliMUONEventReconstructor::RemoveDoubleTracks"
110 to keep only one track among similar ones.
112 Revision 1.8 2000/07/18 16:04:06 gosset
113 AliMUONEventReconstructor package:
114 * a few minor modifications and more comments
116 * right sign for Z of raw clusters
117 * right loop over chambers inside station
118 * symmetrized covariance matrix for measurements (TrackChi2MCS)
119 * right sign of charge in extrapolation (ExtrapToZ)
120 * right zEndAbsorber for Branson correction below 3 degrees
121 * use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
122 * no parameter for AliMUONTrack::Fit() but more fit parameters in Track object
124 Revision 1.7 2000/07/03 12:28:06 gosset
125 Printout at the right place after extrapolation to vertex
127 Revision 1.6 2000/06/30 12:01:06 gosset
128 Correction for hit search in the right chamber (JPC)
130 Revision 1.5 2000/06/30 10:15:48 gosset
131 Changes to EventReconstructor...:
132 precision fit with multiple Coulomb scattering;
133 extrapolation to vertex with Branson correction in absorber (JPC)
135 Revision 1.4 2000/06/27 14:11:36 gosset
136 Corrections against violations of coding conventions
138 Revision 1.3 2000/06/16 07:27:08 gosset
139 To remove problem in running RuleChecker, like in MUON-dev
141 Revision 1.1.2.5 2000/06/16 07:00:26 gosset
142 To remove problem in running RuleChecker
144 Revision 1.1.2.4 2000/06/12 08:00:07 morsch
145 Dummy streamer to solve CINT compilation problem (to be investigated !)
147 Revision 1.1.2.3 2000/06/09 20:59:57 morsch
148 Make includes consistent with new file structure.
150 Revision 1.1.2.2 2000/06/09 12:58:05 gosset
151 Removed comment beginnings in Log sections of .cxx files
152 Suppressed most violations of coding rules
154 Revision 1.1.2.1 2000/06/07 14:44:53 gosset
155 Addition of files for track reconstruction in C++
158 ////////////////////////////////////
160 // MUON event reconstructor in ALICE
162 // This class contains as data:
163 // * the parameters for the event reconstruction
164 // * a pointer to the array of hits to be reconstructed (the event)
165 // * a pointer to the array of segments made with these hits inside each station
166 // * a pointer to the array of reconstructed tracks
168 // It contains as methods, among others:
169 // * MakeEventToBeReconstructed to build the array of hits to be reconstructed
170 // * MakeSegments to build the segments
171 // * MakeTracks to build the tracks
173 ////////////////////////////////////
175 #include <Riostream.h> // for cout
176 #include <stdlib.h> // for exit()
181 #include "AliMUONChamber.h"
182 #include "AliMUONEventReconstructor.h"
183 #include "AliMUONHitForRec.h"
184 #include "AliMUONRawCluster.h"
185 #include "AliMUONRecoEvent.h"
186 #include "AliMUONSegment.h"
187 #include "AliMUONTrack.h"
188 #include "AliMUONTrackHit.h"
190 #include "AliRun.h" // for gAlice
191 #include "AliMUONTrackK.h" //AZ
192 #include <TMatrixD.h> //AZ
194 //************* Defaults parameters for reconstruction
195 static const Double_t kDefaultMinBendingMomentum = 3.0;
196 static const Double_t kDefaultMaxBendingMomentum = 500.0;
197 static const Double_t kDefaultMaxChi2 = 100.0;
198 static const Double_t kDefaultMaxSigma2Distance = 16.0;
199 static const Double_t kDefaultBendingResolution = 0.01;
200 static const Double_t kDefaultNonBendingResolution = 0.144;
201 static const Double_t kDefaultChamberThicknessInX0 = 0.03;
202 // Simple magnetic field:
203 // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
204 // Length and Position from reco_muon.F, with opposite sign:
205 // Length = ZMAGEND-ZCOIL
206 // Position = (ZMAGEND+ZCOIL)/2
207 // to be ajusted differently from real magnetic field ????
208 static const Double_t kDefaultSimpleBValue = 7.0;
209 static const Double_t kDefaultSimpleBLength = 428.0;
210 static const Double_t kDefaultSimpleBPosition = 1019.0;
211 static const Int_t kDefaultRecGeantHits = 0;
212 static const Double_t kDefaultEfficiency = 0.95;
214 static const Int_t kDefaultPrintLevel = 0;
216 ClassImp(AliMUONEventReconstructor) // Class implementation in ROOT context
218 //__________________________________________________________________________
219 AliMUONEventReconstructor::AliMUONEventReconstructor(void)
221 // Constructor for class AliMUONEventReconstructor
222 SetReconstructionParametersToDefaults();
223 fTrackMethod = 1; //AZ - tracking method (1-default, 2-Kalman)
224 // Memory allocation for the TClonesArray of hits for reconstruction
225 // Is 10000 the right size ????
226 fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
227 fNHitsForRec = 0; // really needed or GetEntriesFast sufficient ????
228 // Memory allocation for the TClonesArray's of segments in stations
229 // Is 2000 the right size ????
230 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
231 fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
232 fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
234 // Memory allocation for the TClonesArray of reconstructed tracks
235 // Is 10 the right size ????
236 fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
237 fNRecTracks = 0; // really needed or GetEntriesFast sufficient ????
238 // Memory allocation for the TClonesArray of hits on reconstructed tracks
239 // Is 100 the right size ????
240 fRecTrackHitsPtr = new TClonesArray("AliMUONTrack", 100);
241 fNRecTrackHits = 0; // really needed or GetEntriesFast sufficient ????
243 // Sign of fSimpleBValue according to sign of Bx value at (50,50,950).
245 x[0] = 50.; x[1] = 50.; x[2] = 950.;
246 gAlice->Field()->Field(x, b);
247 fSimpleBValue = TMath::Sign(fSimpleBValue,(Double_t) b[0]);
248 // See how to get fSimple(BValue, BLength, BPosition)
249 // automatically calculated from the actual magnetic field ????
251 if (fPrintLevel >= 0) {
252 cout << "AliMUONEventReconstructor constructed with defaults" << endl; Dump();
253 cout << endl << "Magnetic field from root file:" << endl;
254 gAlice->Field()->Dump();
258 // Initializions for GEANT background events
261 fBkgGeantParticles = 0;
264 fBkgGeantEventNumber = -1;
266 // Initialize to 0 pointers to RecoEvent, tree and tree file
274 AliMUONEventReconstructor::AliMUONEventReconstructor (const AliMUONEventReconstructor& Reconstructor)
276 // Dummy copy constructor
279 AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEventReconstructor& Reconstructor)
281 // Dummy assignment operator
285 //__________________________________________________________________________
286 AliMUONEventReconstructor::~AliMUONEventReconstructor(void)
288 // Destructor for class AliMUONEventReconstructor
293 // if (fEventTree) delete fEventTree;
294 if (fRecoEvent) delete fRecoEvent;
295 delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
296 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
297 delete fSegmentsPtr[st]; // Correct destruction of everything ????
301 //__________________________________________________________________________
302 void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
304 // Set reconstruction parameters to default values
305 // Would be much more convenient with a structure (or class) ????
306 fMinBendingMomentum = kDefaultMinBendingMomentum;
307 fMaxBendingMomentum = kDefaultMaxBendingMomentum;
308 fMaxChi2 = kDefaultMaxChi2;
309 fMaxSigma2Distance = kDefaultMaxSigma2Distance;
311 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
312 // ******** Parameters for making HitsForRec
314 // like in TRACKF_STAT:
315 // 2 degrees for stations 1 and 2, or ch(0...) from 0 to 3;
316 // 30 cm for stations 3 to 5, or ch(0...) from 4 to 9
317 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
318 if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
319 2.0 * TMath::Pi() / 180.0;
320 else fRMin[ch] = 30.0;
321 // maximum radius at 10 degrees and Z of chamber
322 fRMax[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
323 10.0 * TMath::Pi() / 180.0;
326 // ******** Parameters for making segments
327 // should be parametrized ????
328 // according to interval between chambers in a station ????
329 // Maximum distance in non bending plane
330 // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
332 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
333 fSegmentMaxDistNonBending[st] = 5. * 0.22;
334 // Maximum distance in bending plane:
335 // values from TRACKF_STAT, corresponding to (J psi 20cm),
336 // scaled to the real distance between chambers in a station
337 fSegmentMaxDistBending[0] = 1.5 *
338 ((&(pMUON->Chamber(1)))->Z() - (&(pMUON->Chamber(0)))->Z()) / 20.0;
339 fSegmentMaxDistBending[1] = 1.5 *
340 ((&(pMUON->Chamber(3)))->Z() - (&(pMUON->Chamber(2)))->Z()) / 20.0;
341 fSegmentMaxDistBending[2] = 3.0 *
342 ((&(pMUON->Chamber(5)))->Z() - (&(pMUON->Chamber(4)))->Z()) / 20.0;
343 fSegmentMaxDistBending[3] = 6.0 *
344 ((&(pMUON->Chamber(7)))->Z() - (&(pMUON->Chamber(6)))->Z()) / 20.0;
345 fSegmentMaxDistBending[4] = 6.0 *
346 ((&(pMUON->Chamber(9)))->Z() - (&(pMUON->Chamber(8)))->Z()) / 20.0;
348 fBendingResolution = kDefaultBendingResolution;
349 fNonBendingResolution = kDefaultNonBendingResolution;
350 fChamberThicknessInX0 = kDefaultChamberThicknessInX0;
351 fSimpleBValue = kDefaultSimpleBValue;
352 fSimpleBLength = kDefaultSimpleBLength;
353 fSimpleBPosition = kDefaultSimpleBPosition;
354 fRecGeantHits = kDefaultRecGeantHits;
355 fEfficiency = kDefaultEfficiency;
356 fPrintLevel = kDefaultPrintLevel;
360 //__________________________________________________________________________
361 Double_t AliMUONEventReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const
363 // Returns impact parameter at vertex in bending plane (cm),
364 // from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
365 // using simple values for dipole magnetic field.
366 // The sign of "BendingMomentum" is the sign of the charge.
367 return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
371 //__________________________________________________________________________
372 Double_t AliMUONEventReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam) const
374 // Returns signed bending momentum in bending plane (GeV/c),
375 // the sign being the sign of the charge for particles moving forward in Z,
376 // from the impact parameter "ImpactParam" at vertex in bending plane (cm),
377 // using simple values for dipole magnetic field.
378 return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
382 //__________________________________________________________________________
383 void AliMUONEventReconstructor::SetBkgGeantFile(Text_t *BkgGeantFileName)
385 // Set background file ... for GEANT hits
386 // Must be called after having loaded the firts signal event
387 if (fPrintLevel >= 0) {
388 cout << "Enter SetBkgGeantFile with BkgGeantFileName ``"
389 << BkgGeantFileName << "''" << endl;}
390 if (strlen(BkgGeantFileName)) {
391 // BkgGeantFileName not empty: try to open the file
392 if (fPrintLevel >= 2) {cout << "Before File(Bkg)" << endl; gDirectory->Dump();}
393 fBkgGeantFile = new TFile(BkgGeantFileName);
394 if (fPrintLevel >= 2) {cout << "After File(Bkg)" << endl; gDirectory->Dump();}
395 if (fBkgGeantFile-> IsOpen()) {
396 if (fPrintLevel >= 0) {
397 cout << "Background for GEANT hits in file: ``" << BkgGeantFileName
398 << "'' successfully opened" << endl;}
401 cout << "Background for GEANT hits in file: " << BkgGeantFileName << endl;
402 cout << "NOT FOUND: EXIT" << endl;
403 exit(0); // right instruction for exit ????
405 // Arrays for "particles" and "hits"
406 fBkgGeantParticles = new TClonesArray("TParticle", 200);
407 fBkgGeantHits = new TClonesArray("AliMUONHit", 2000);
408 // Event number to -1 for initialization
409 fBkgGeantEventNumber = -1;
410 // Back to the signal file:
411 // first signal event must have been loaded previously,
412 // otherwise, Segmentation violation at the next instruction
413 // How is it possible to do smething better ????
414 ((gAlice->TreeK())->GetCurrentFile())->cd();
415 if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
420 //__________________________________________________________________________
421 void AliMUONEventReconstructor::NextBkgGeantEvent(void)
423 // Get next event in background file for GEANT hits
424 // Goes back to event number 0 when end of file is reached
427 if (fPrintLevel >= 0) {
428 cout << "Enter NextBkgGeantEvent" << endl;}
429 // Clean previous event
430 if(fBkgGeantTK) delete fBkgGeantTK;
432 if(fBkgGeantParticles) fBkgGeantParticles->Clear();
433 if(fBkgGeantTH) delete fBkgGeantTH;
435 if(fBkgGeantHits) fBkgGeantHits->Clear();
436 // Increment event number
437 fBkgGeantEventNumber++;
438 // Get access to Particles and Hits for event from background file
439 if (fPrintLevel >= 2) {cout << "Before cd(Bkg)" << endl; gDirectory->Dump();}
441 if (fPrintLevel >= 2) {cout << "After cd(Bkg)" << endl; gDirectory->Dump();}
442 // Particles: TreeK for event and branch "Particles"
443 sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
444 fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
446 if (fPrintLevel >= 0) {
447 cout << "Cannot find Kine Tree for background event: " <<
448 fBkgGeantEventNumber << endl;
449 cout << "Goes back to event 0" << endl;
451 fBkgGeantEventNumber = 0;
452 sprintf(treeName, "TreeK%d", fBkgGeantEventNumber);
453 fBkgGeantTK = (TTree*)gDirectory->Get(treeName);
455 cout << "ERROR: cannot find Kine Tree for background event: " <<
456 fBkgGeantEventNumber << endl;
461 fBkgGeantTK->SetBranchAddress("Particles", &fBkgGeantParticles);
462 fBkgGeantTK->GetEvent(0); // why event 0 ???? necessary ????
463 // Hits: TreeH for event and branch "MUON"
464 sprintf(treeName, "TreeH%d", fBkgGeantEventNumber);
465 fBkgGeantTH = (TTree*)gDirectory->Get(treeName);
467 cout << "ERROR: cannot find Hits Tree for background event: " <<
468 fBkgGeantEventNumber << endl;
471 if (fBkgGeantTH && fBkgGeantHits) {
472 branch = fBkgGeantTH->GetBranch("MUON");
473 if (branch) branch->SetAddress(&fBkgGeantHits);
475 fBkgGeantTH->GetEntries(); // necessary ????
476 // Back to the signal file
477 ((gAlice->TreeK())->GetCurrentFile())->cd();
478 if (fPrintLevel >= 2) {cout << "After cd(gAlice)" << endl; gDirectory->Dump();}
482 //__________________________________________________________________________
483 void AliMUONEventReconstructor::EventReconstruct(void)
485 // To reconstruct one event
486 if (fPrintLevel >= 1) cout << "enter EventReconstruct" << endl;
487 MakeEventToBeReconstructed();
493 //__________________________________________________________________________
494 void AliMUONEventReconstructor::ResetHitsForRec(void)
496 // To reset the array and the number of HitsForRec,
497 // and also the number of HitsForRec
498 // and the index of the first HitForRec per chamber
499 if (fHitsForRecPtr) fHitsForRecPtr->Clear();
501 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++)
502 fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0;
506 //__________________________________________________________________________
507 void AliMUONEventReconstructor::ResetSegments(void)
509 // To reset the TClonesArray of segments and the number of Segments
511 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
512 if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
518 //__________________________________________________________________________
519 void AliMUONEventReconstructor::ResetTracks(void)
521 // To reset the TClonesArray of reconstructed tracks
522 if (fRecTracksPtr) fRecTracksPtr->Delete();
523 // Delete in order that the Track destructors are called,
524 // hence the space for the TClonesArray of pointers to TrackHit's is freed
529 //__________________________________________________________________________
530 void AliMUONEventReconstructor::ResetTrackHits(void)
532 // To reset the TClonesArray of hits on reconstructed tracks
533 if (fRecTrackHitsPtr) fRecTrackHitsPtr->Clear();
538 //__________________________________________________________________________
539 void AliMUONEventReconstructor::MakeEventToBeReconstructed(void)
541 // To make the list of hits to be reconstructed,
542 // either from the GEANT hits or from the raw clusters
543 // according to the parameter set for the reconstructor
544 if (fPrintLevel >= 1) cout << "enter MakeEventToBeReconstructed" << endl;
546 if (fRecGeantHits == 1) {
547 // Reconstruction from GEANT hits
548 // Back to the signal file
549 ((gAlice->TreeK())->GetCurrentFile())->cd();
551 // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
552 // Security on MUON ????
553 AddHitsForRecFromGEANT(gAlice->TreeH());
555 AddHitsForRecFromBkgGEANT(fBkgGeantTH, fBkgGeantHits);
556 // Sort HitsForRec in increasing order with respect to chamber number
557 SortHitsForRecWithIncreasingChamber();
560 // Reconstruction from raw clusters
561 // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
562 // Security on MUON ????
563 // TreeR assumed to be be "prepared" in calling function
564 // by "MUON->GetTreeR(nev)" ????
565 TTree *treeR = gAlice->TreeR();
566 AddHitsForRecFromRawClusters(treeR);
567 // No sorting: it is done automatically in the previous function
569 if (fPrintLevel >= 10) {
570 cout << "end of MakeEventToBeReconstructed" << endl;
571 cout << "NHitsForRec: " << fNHitsForRec << endl;
572 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
573 cout << "chamber(0...): " << ch
574 << " NHitsForRec: " << fNHitsForRecPerChamber[ch]
575 << " index(first HitForRec): " << fIndexOfFirstHitForRecPerChamber[ch]
577 for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
578 hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
580 cout << "HitForRec index(0...): " << hit << endl;
581 ((*fHitsForRecPtr)[hit])->Dump();
588 //__________________________________________________________________________
589 void AliMUONEventReconstructor::AddHitsForRecFromGEANT(TTree *TH)
591 // To add to the list of hits for reconstruction
592 // the GEANT signal hits from a hit tree TH.
593 Int_t hitBits, chamBits; //AZ
594 if (fPrintLevel >= 2)
595 cout << "enter AddHitsForRecFromGEANT with TH: " << TH << endl;
596 if (TH == NULL) return;
597 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
598 // Security on MUON ????
599 // See whether it could be the same for signal and background ????
600 // Loop over tracks in tree
601 Int_t ntracks = (Int_t) TH->GetEntries();
602 if (fPrintLevel >= 2)
603 cout << "ntracks: " << ntracks << endl;
605 for (Int_t track = 0; track < ntracks; track++) {
612 Int_t itrack = track; //AZ
613 for (AliMUONHit* mHit = (AliMUONHit*) pMUON->FirstHit(-1);
615 mHit = (AliMUONHit*) pMUON->NextHit(), hit++) {
616 Int_t ipart = TMath::Abs ((Int_t) mHit->Particle()); //AZ
617 //itrack = mHit->Track(); //AZ
618 //AZNewHitForRecFromGEANT(mHit,track, hit, 1);
619 if (NewHitForRecFromGEANT(mHit,track, hit, 1) && ipart == 13
620 //if (NewHitForRecFromGEANT(mHit,itrack-1, hit, 1) && ipart == 13
621 && itrack <= 2) chamBits |= BIT(mHit->Chamber()-1); //AZ - set bit
623 if (chamBits&3 && chamBits>>2&3 && chamBits>>4&3 && chamBits>>6&3 &&
624 chamBits>>8&3 && ((chamBits>>6&3)==3 || (chamBits>>8&3)==3))
626 //if (chamBits&3 && chamBits>>2&3 && chamBits>>4&3 && chamBits>>6&3 &&
627 // chamBits>>8&3 && ((chamBits>>6&3)==3 || (chamBits>>8&3)==3) &&
628 // ((chamBits&3)==3 || (chamBits>>2&3)==3)) fMuons += 1;
629 } // end of track loop
633 //__________________________________________________________________________
634 void AliMUONEventReconstructor::AddHitsForRecFromBkgGEANT(TTree *TH, TClonesArray *Hits)
636 // To add to the list of hits for reconstruction
637 // the GEANT background hits from a hit tree TH and a pointer Hits to a hit list.
638 // How to have only one function "AddHitsForRecFromGEANT" ????
639 if (fPrintLevel >= 2)
640 cout << "enter AddHitsForRecFromBkgGEANT with TH: " << TH << endl;
641 if (TH == NULL) return;
642 // Loop over tracks in tree
643 Int_t ntracks = (Int_t) TH->GetEntries();
644 if (fPrintLevel >= 2)
645 cout << "ntracks: " << ntracks << endl;
646 for (Int_t track = 0; track < ntracks; track++) {
647 if (Hits) Hits->Clear();
650 for (Int_t hit = 0; hit < Hits->GetEntriesFast(); hit++) {
651 NewHitForRecFromGEANT((AliMUONHit*) (*Hits)[hit], track, hit, 0);
653 } // end of track loop
657 //__________________________________________________________________________
658 AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* Hit, Int_t TrackNumber, Int_t HitNumber, Int_t Signal)
660 // To make a new hit for reconstruction from a GEANT hit pointed to by "Hit",
661 // with hit number "HitNumber" in the track numbered "TrackNumber",
662 // either from signal ("Signal" = 1) or background ("Signal" = 0) event.
663 // Selects hits in tracking (not trigger) chambers.
664 // Takes into account the efficiency (fEfficiency)
665 // and the smearing from resolution (fBendingResolution and fNonBendingResolution).
666 // Adds a condition on the radius between RMin and RMax
667 // to better simulate the real chambers.
668 // Returns the pointer to the new hit for reconstruction,
669 // or NULL in case of inefficiency or non tracking chamber or bad radius.
670 // No condition on at most 20 cm from a muon from a resonance
671 // like in Fortran TRACKF_STAT.
672 AliMUONHitForRec* hitForRec;
673 Double_t bendCoor, nonBendCoor, radius;
674 Int_t chamber = Hit->Chamber() - 1; // chamber(0...)
675 // only in tracking chambers (fChamber starts at 1)
676 if (chamber >= kMaxMuonTrackingChambers) return NULL;
677 // only if hit is efficient (keep track for checking ????)
678 if (gRandom->Rndm() > fEfficiency) return NULL;
679 // only if radius between RMin and RMax
681 nonBendCoor = Hit->X();
682 radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor));
683 // This cut is not needed with a realistic chamber geometry !!!!
684 // if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
685 // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's
686 hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit);
688 // add smearing from resolution
689 hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution));
690 hitForRec->SetNonBendingCoor(nonBendCoor
691 + gRandom->Gaus(0., fNonBendingResolution));
692 // // !!!! without smearing
693 // hitForRec->SetBendingCoor(bendCoor);
694 // hitForRec->SetNonBendingCoor(nonBendCoor);
695 // more information into HitForRec
696 // resolution: angular effect to be added here ????
697 hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
698 hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
700 hitForRec->SetHitNumber(HitNumber);
701 hitForRec->SetTHTrack(TrackNumber);
702 hitForRec->SetGeantSignal(Signal);
703 if (fPrintLevel >= 10) {
704 cout << "track: " << TrackNumber << " hit: " << HitNumber << endl;
706 cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
711 //__________________________________________________________________________
712 void AliMUONEventReconstructor::SortHitsForRecWithIncreasingChamber()
714 // Sort HitsForRec's in increasing order with respect to chamber number.
715 // Uses the function "Compare".
716 // Update the information for HitsForRec per chamber too.
717 Int_t ch, nhits, prevch;
718 fHitsForRecPtr->Sort();
719 for (ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
720 fNHitsForRecPerChamber[ch] = 0;
721 fIndexOfFirstHitForRecPerChamber[ch] = 0;
723 prevch = 0; // previous chamber
724 nhits = 0; // number of hits in current chamber
725 // Loop over HitsForRec
726 for (Int_t hit = 0; hit < fNHitsForRec; hit++) {
727 // chamber number (0...)
728 ch = ((AliMUONHitForRec*) ((*fHitsForRecPtr)[hit]))->GetChamberNumber();
729 // increment number of hits in current chamber
730 (fNHitsForRecPerChamber[ch])++;
731 // update index of first HitForRec in current chamber
732 // if chamber number different from previous one
734 fIndexOfFirstHitForRecPerChamber[ch] = hit;
741 // //__________________________________________________________________________
742 // void AliMUONEventReconstructor::AddHitsForRecFromCathodeCorrelations(TTree* TC)
744 // // OLD VERSION WHEN ONE ONE WAS USING SO CALLED CATHODE CORRELATIONS
745 // // To add to the list of hits for reconstruction
746 // // the (cathode correlated) raw clusters
747 // // No condition added, like in Fortran TRACKF_STAT,
748 // // on the radius between RMin and RMax.
749 // AliMUONHitForRec *hitForRec;
750 // if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromCathodeCorrelations" << endl;
751 // AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
752 // // Security on MUON ????
753 // // Loop over tracking chambers
754 // for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
755 // // number of HitsForRec to 0 for the chamber
756 // fNHitsForRecPerChamber[ch] = 0;
757 // // index of first HitForRec for the chamber
758 // if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
759 // else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
760 // TClonesArray *reconst_hits = MUON->ReconstHitsAddress(ch);
761 // MUON->ResetReconstHits();
763 // Int_t ncor = (Int_t)reconst_hits->GetEntries();
764 // // Loop over (cathode correlated) raw clusters
765 // for (Int_t cor = 0; cor < ncor; cor++) {
766 // AliMUONReconstHit * mCor =
767 // (AliMUONReconstHit*) reconst_hits->UncheckedAt(cor);
768 // // new AliMUONHitForRec from (cathode correlated) raw cluster
769 // // and increment number of AliMUONHitForRec's (total and in chamber)
770 // hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(mCor);
772 // (fNHitsForRecPerChamber[ch])++;
773 // // more information into HitForRec
774 // hitForRec->SetChamberNumber(ch);
775 // hitForRec->SetHitNumber(cor);
776 // // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
777 // // could (should) be more exact from chamber geometry ????
778 // hitForRec->SetZ(-(&(MUON->Chamber(ch)))->Z());
779 // if (fPrintLevel >= 10) {
780 // cout << "chamber (0...): " << ch <<
781 // " cathcorrel (0...): " << cor << endl;
783 // cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
784 // hitForRec->Dump();}
785 // } // end of cluster loop
786 // } // end of chamber loop
790 //__________________________________________________________________________
791 void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
793 // To add to the list of hits for reconstruction all the raw clusters
794 // No condition added, like in Fortran TRACKF_STAT,
795 // on the radius between RMin and RMax.
796 AliMUONHitForRec *hitForRec;
797 AliMUONRawCluster *clus;
798 Int_t iclus, nclus, nTRentries;
799 TClonesArray *rawclusters;
800 if (fPrintLevel >= 1) cout << "enter AddHitsForRecFromRawClusters" << endl;
801 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
802 // Security on MUON ????
803 pMUON->ResetRawClusters();
804 nTRentries = Int_t(TR->GetEntries());
805 if (nTRentries != 1) {
806 cout << "Error in AliMUONEventReconstructor::AddHitsForRecFromRawClusters"
808 cout << "nTRentries = " << nTRentries << " not equal to 1" << endl;
811 TR->GetEvent(0); // only one entry
812 // Loop over tracking chambers
813 for (Int_t ch = 0; ch < kMaxMuonTrackingChambers; ch++) {
814 // number of HitsForRec to 0 for the chamber
815 fNHitsForRecPerChamber[ch] = 0;
816 // index of first HitForRec for the chamber
817 if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
818 else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
819 rawclusters = pMUON->RawClustAddress(ch);
820 // pMUON->ResetRawClusters();
821 // TR->GetEvent((Int_t) (TR->GetEntries()) - 1); // to be checked ????
822 nclus = (Int_t) (rawclusters->GetEntries());
823 // Loop over (cathode correlated) raw clusters
824 for (iclus = 0; iclus < nclus; iclus++) {
825 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
826 // new AliMUONHitForRec from raw cluster
827 // and increment number of AliMUONHitForRec's (total and in chamber)
828 hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
830 (fNHitsForRecPerChamber[ch])++;
831 // more information into HitForRec
832 // resolution: info should be already in raw cluster and taken from it ????
833 hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
834 hitForRec->SetNonBendingReso2(fNonBendingResolution * fNonBendingResolution);
835 // original raw cluster
836 hitForRec->SetChamberNumber(ch);
837 hitForRec->SetHitNumber(iclus);
838 // Z coordinate of the raw cluster (cm)
839 hitForRec->SetZ(clus->fZ[0]);
840 if (fPrintLevel >= 10) {
841 cout << "chamber (0...): " << ch <<
842 " raw cluster (0...): " << iclus << endl;
844 cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
846 } // end of cluster loop
847 } // end of chamber loop
848 SortHitsForRecWithIncreasingChamber(); //AZ
852 //__________________________________________________________________________
853 void AliMUONEventReconstructor::MakeSegments(void)
855 // To make the list of segments in all stations,
856 // from the list of hits to be reconstructed
857 if (fPrintLevel >= 1) cout << "enter MakeSegments" << endl;
859 // Loop over stations
860 Int_t nb = (fTrackMethod == 2) ? 3 : 0; //AZ
861 //AZ for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
862 for (Int_t st = nb; st < kMaxMuonTrackingStations; st++) //AZ
863 MakeSegmentsPerStation(st);
864 if (fPrintLevel >= 10) {
865 cout << "end of MakeSegments" << endl;
866 for (Int_t st = 0; st < kMaxMuonTrackingStations; st++) {
867 cout << "station(0...): " << st
868 << " Segments: " << fNSegments[st]
871 seg < fNSegments[st];
873 cout << "Segment index(0...): " << seg << endl;
874 ((*fSegmentsPtr[st])[seg])->Dump();
881 //__________________________________________________________________________
882 void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
884 // To make the list of segments in station number "Station" (0...)
885 // from the list of hits to be reconstructed.
886 // Updates "fNSegments"[Station].
887 // Segments in stations 4 and 5 are sorted
888 // according to increasing absolute value of "impact parameter"
889 AliMUONHitForRec *hit1Ptr, *hit2Ptr;
890 AliMUONSegment *segment;
892 Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
893 impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
894 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
895 if (fPrintLevel >= 1)
896 cout << "enter MakeSegmentsPerStation (0...) " << Station << endl;
897 // first and second chambers (0...) in the station
898 Int_t ch1 = 2 * Station;
900 // variable true for stations downstream of the dipole:
901 // Station(0..4) equal to 3 or 4
902 if ((Station == 3) || (Station == 4)) {
904 // maximum impact parameter (cm) according to fMinBendingMomentum...
906 TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
907 // minimum impact parameter (cm) according to fMaxBendingMomentum...
909 TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum));
911 else last2st = kFALSE;
912 // extrapolation factor from Z of first chamber to Z of second chamber
913 // dZ to be changed to take into account fine structure of chambers ????
914 Double_t extrapFact =
915 (&(pMUON->Chamber(ch2)))->Z() / (&(pMUON->Chamber(ch1)))->Z();
916 // index for current segment
917 Int_t segmentIndex = 0;
918 // Loop over HitsForRec in the first chamber of the station
919 for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
920 hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
922 // pointer to the HitForRec
923 hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
925 // on the straight line joining the HitForRec to the vertex (0,0,0),
926 // to the Z of the second chamber of the station
927 extBendCoor = extrapFact * hit1Ptr->GetBendingCoor();
928 extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
929 // Loop over HitsForRec in the second chamber of the station
930 for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
931 hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
933 // pointer to the HitForRec
934 hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]);
935 // absolute values of distances, in bending and non bending planes,
936 // between the HitForRec in the second chamber
937 // and the previous extrapolation
938 distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
939 distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
942 bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
943 (hit1Ptr->GetZ() - hit2Ptr->GetZ());
944 // absolute value of impact parameter
946 TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope);
948 // check for distances not too large,
949 // and impact parameter not too big if stations downstream of the dipole.
950 // Conditions "distBend" and "impactParam" correlated for these stations ????
951 if ((distBend < fSegmentMaxDistBending[Station]) &&
952 (distNonBend < fSegmentMaxDistNonBending[Station]) &&
953 (!last2st || (impactParam < maxImpactParam)) &&
954 (!last2st || (impactParam > minImpactParam))) {
956 segment = new ((*fSegmentsPtr[Station])[segmentIndex])
957 AliMUONSegment(hit1Ptr, hit2Ptr);
958 // update "link" to this segment from the hit in the first chamber
959 if (hit1Ptr->GetNSegments() == 0)
960 hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
961 hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
962 if (fPrintLevel >= 10) {
963 cout << "segmentIndex(0...): " << segmentIndex
964 << " distBend: " << distBend
965 << " distNonBend: " << distNonBend
968 cout << "HitForRec in first chamber" << endl;
970 cout << "HitForRec in second chamber" << endl;
973 // increment index for current segment
977 } // for (Int_t hit1...
978 fNSegments[Station] = segmentIndex;
979 // Sorting according to "impact parameter" if station(1..5) 4 or 5,
980 // i.e. Station(0..4) 3 or 4, using the function "Compare".
981 // After this sorting, it is impossible to use
982 // the "fNSegments" and "fIndexOfFirstSegment"
983 // of the HitForRec in the first chamber to explore all segments formed with it.
984 // Is this sorting really needed ????
985 if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
986 if (fPrintLevel >= 1) cout << "Station: " << Station << " NSegments: "
987 << fNSegments[Station] << endl;
991 //__________________________________________________________________________
992 void AliMUONEventReconstructor::MakeTracks(void)
994 // To make the tracks,
995 // from the list of segments and points in all stations
996 if (fPrintLevel >= 1) cout << "enter MakeTracks" << endl;
997 // The order may be important for the following Reset's
1000 if (fTrackMethod == 2) { //AZ - Kalman filter
1001 MakeTrackCandidatesK();
1002 // Follow tracks in stations(1..) 3, 2 and 1
1004 // Remove double tracks
1005 RemoveDoubleTracksK();
1006 // Propagate tracks to the vertex thru absorber
1009 // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
1010 MakeTrackCandidates();
1011 // Follow tracks in stations(1..) 3, 2 and 1
1013 // Remove double tracks
1014 RemoveDoubleTracks();
1019 //__________________________________________________________________________
1020 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
1022 // To make track candidates with two segments in stations(1..) 4 and 5,
1023 // the first segment being pointed to by "BegSegment".
1024 // Returns the number of such track candidates.
1025 Int_t endStation, iEndSegment, nbCan2Seg;
1026 AliMUONSegment *endSegment, *extrapSegment;
1027 AliMUONTrack *recTrack;
1029 if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesWithTwoSegments" << endl;
1030 // Station for the end segment
1031 endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
1032 // multiple scattering factor corresponding to one chamber
1033 mcsFactor = 0.0136 /
1034 GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
1035 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
1036 // linear extrapolation to end station
1038 BegSegment->CreateSegmentFromLinearExtrapToStation(endStation, mcsFactor);
1039 // number of candidates with 2 segments to 0
1041 // Loop over segments in the end station
1042 for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
1043 // pointer to segment
1044 endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
1045 // test compatibility between current segment and "extrapSegment"
1046 // 4 because 4 quantities in chi2
1048 NormalizedChi2WithSegment(extrapSegment,
1049 fMaxSigma2Distance)) <= 4.0) {
1050 // both segments compatible:
1051 // make track candidate from "begSegment" and "endSegment"
1052 if (fPrintLevel >= 2)
1053 cout << "TrackCandidate with Segment " << iEndSegment <<
1054 " in Station(0..) " << endStation << endl;
1055 // flag for both segments in one track:
1056 // to be done in track constructor ????
1057 BegSegment->SetInTrack(kTRUE);
1058 endSegment->SetInTrack(kTRUE);
1059 recTrack = new ((*fRecTracksPtr)[fNRecTracks])
1060 AliMUONTrack(BegSegment, endSegment, this);
1062 if (fPrintLevel >= 10) recTrack->RecursiveDump();
1063 // increment number of track candidates with 2 segments
1066 } // for (iEndSegment = 0;...
1067 delete extrapSegment; // should not delete HitForRec's it points to !!!!
1071 //__________________________________________________________________________
1072 Int_t AliMUONEventReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
1074 // To make track candidates with one segment and one point
1075 // in stations(1..) 4 and 5,
1076 // the segment being pointed to by "BegSegment".
1077 Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
1078 AliMUONHitForRec *extrapHitForRec, *hit;
1079 AliMUONTrack *recTrack;
1081 if (fPrintLevel >= 1)
1082 cout << "enter MakeTrackCandidatesWithOneSegmentAndOnePoint" << endl;
1083 // station for the end point
1084 endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
1085 // multiple scattering factor corresponding to one chamber
1086 mcsFactor = 0.0136 /
1087 GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
1088 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
1089 // first and second chambers(0..) in the end station
1090 ch1 = 2 * endStation;
1092 // number of candidates to 0
1094 // Loop over chambers of the end station
1095 for (ch = ch2; ch >= ch1; ch--) {
1096 // linear extrapolation to chamber
1098 BegSegment->CreateHitForRecFromLinearExtrapToChamber(ch, mcsFactor);
1099 // limits for the hit index in the loop
1100 iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
1101 iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
1102 // Loop over HitForRec's in the chamber
1103 for (iHit = iHitMin; iHit < iHitMax; iHit++) {
1104 // pointer to HitForRec
1105 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
1106 // test compatibility between current HitForRec and "extrapHitForRec"
1107 // 2 because 2 quantities in chi2
1109 NormalizedChi2WithHitForRec(extrapHitForRec,
1110 fMaxSigma2Distance)) <= 2.0) {
1111 // both HitForRec's compatible:
1112 // make track candidate from begSegment and current HitForRec
1113 if (fPrintLevel >= 2)
1114 cout << "TrackCandidate with HitForRec " << iHit <<
1115 " in Chamber(0..) " << ch << endl;
1116 // flag for beginning segments in one track:
1117 // to be done in track constructor ????
1118 BegSegment->SetInTrack(kTRUE);
1119 recTrack = new ((*fRecTracksPtr)[fNRecTracks])
1120 AliMUONTrack(BegSegment, hit, this);
1121 // the right place to eliminate "double counting" ???? how ????
1123 if (fPrintLevel >= 10) recTrack->RecursiveDump();
1124 // increment number of track candidates
1127 } // for (iHit = iHitMin;...
1128 delete extrapHitForRec;
1129 } // for (ch = ch2;...
1130 return nbCan1Seg1Hit;
1133 //__________________________________________________________________________
1134 void AliMUONEventReconstructor::MakeTrackCandidates(void)
1136 // To make track candidates
1137 // with at least 3 aligned points in stations(1..) 4 and 5
1138 // (two Segment's or one Segment and one HitForRec)
1139 Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
1140 AliMUONSegment *begSegment;
1141 if (fPrintLevel >= 1) cout << "enter MakeTrackCandidates" << endl;
1142 // Loop over stations(1..) 5 and 4 for the beginning segment
1143 for (begStation = 4; begStation > 2; begStation--) {
1144 // Loop over segments in the beginning station
1145 for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
1146 // pointer to segment
1147 begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
1148 if (fPrintLevel >= 2)
1149 cout << "look for TrackCandidate's with Segment " << iBegSegment <<
1150 " in Station(0..) " << begStation << endl;
1151 // Look for track candidates with two segments,
1152 // "begSegment" and all compatible segments in other station.
1153 // Only for beginning station(1..) 5
1154 // because candidates with 2 segments have to looked for only once.
1155 if (begStation == 4)
1156 nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
1157 // Look for track candidates with one segment and one point,
1158 // "begSegment" and all compatible HitForRec's in other station.
1159 // Only if "begSegment" does not belong already to a track candidate.
1160 // Is that a too strong condition ????
1161 if (!(begSegment->GetInTrack()))
1162 nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
1163 } // for (iBegSegment = 0;...
1164 } // for (begStation = 4;...
1168 //__________________________________________________________________________
1169 void AliMUONEventReconstructor::FollowTracks(void)
1171 // Follow tracks in stations(1..) 3, 2 and 1
1172 // too long: should be made more modular !!!!
1173 AliMUONHitForRec *bestHit, *extrapHit, *extrapCorrHit, *hit;
1174 AliMUONSegment *bestSegment, *extrapSegment, *extrapCorrSegment, *segment;
1175 AliMUONTrack *track, *nextTrack;
1176 AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
1177 // -1 to avoid compilation warnings
1178 Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
1179 Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
1180 Double_t bendingMomentum, chi2Norm = 0.;
1181 AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
1182 // local maxSigma2Distance, for easy increase in testing
1183 maxSigma2Distance = fMaxSigma2Distance;
1184 if (fPrintLevel >= 2)
1185 cout << "enter FollowTracks" << endl;
1186 // Loop over track candidates
1187 track = (AliMUONTrack*) fRecTracksPtr->First();
1190 // Follow function for each track candidate ????
1192 nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
1193 if (fPrintLevel >= 2)
1194 cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
1195 // Fit track candidate
1196 track->SetFitMCS(0); // without multiple Coulomb scattering
1197 track->SetFitNParam(3); // with 3 parameters (X = Y = 0)
1198 track->SetFitStart(0); // from parameters at vertex
1200 if (fPrintLevel >= 10) {
1201 cout << "FollowTracks: track candidate(0..): " << trackIndex
1202 << " after fit in stations(0..) 3 and 4" << endl;
1203 track->RecursiveDump();
1205 // Loop over stations(1..) 3, 2 and 1
1206 // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
1207 // otherwise: majority coincidence 2 !!!!
1208 for (station = 2; station >= 0; station--) {
1209 // Track parameters at first track hit (smallest Z)
1210 trackParam1 = ((AliMUONTrackHit*)
1211 (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1212 // extrapolation to station
1213 trackParam1->ExtrapToStation(station, trackParam);
1214 extrapSegment = new AliMUONSegment(); // empty segment
1215 extrapCorrSegment = new AliMUONSegment(); // empty corrected segment
1216 // multiple scattering factor corresponding to one chamber
1217 // and momentum in bending plane (not total)
1218 mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
1219 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
1220 // Z difference from previous station
1221 dZ1 = (&(pMUON->Chamber(2 * station)))->Z() -
1222 (&(pMUON->Chamber(2 * station + 2)))->Z();
1223 // Z difference between the two previous stations
1224 dZ2 = (&(pMUON->Chamber(2 * station + 2)))->Z() -
1225 (&(pMUON->Chamber(2 * station + 4)))->Z();
1226 // Z difference between the two chambers in the previous station
1227 dZ3 = (&(pMUON->Chamber(2 * station)))->Z() -
1228 (&(pMUON->Chamber(2 * station + 1)))->Z();
1229 extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
1231 SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
1232 extrapSegment->UpdateFromStationTrackParam
1233 (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
1234 trackParam1->GetInverseBendingMomentum());
1235 // same thing for corrected segment
1236 // better to use copy constructor, after checking that it works properly !!!!
1237 extrapCorrSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
1239 SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
1240 extrapCorrSegment->UpdateFromStationTrackParam
1241 (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
1242 trackParam1->GetInverseBendingMomentum());
1245 if (fPrintLevel >= 10) {
1246 cout << "FollowTracks: track candidate(0..): " << trackIndex
1247 << " Look for segment in station(0..): " << station << endl;
1249 // Loop over segments in station
1250 for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
1251 // Look for best compatible Segment in station
1252 // should consider all possibilities ????
1253 // multiple scattering ????
1254 // separation in 2 functions: Segment and HitForRec ????
1255 segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
1256 // correction of corrected segment (fBendingCoor and fNonBendingCoor)
1257 // according to real Z value of "segment" and slopes of "extrapSegment"
1259 SetBendingCoor(extrapSegment->GetBendingCoor() +
1260 extrapSegment->GetBendingSlope() *
1261 (segment->GetHitForRec1()->GetZ() -
1262 (&(pMUON->Chamber(2 * station)))->Z()));
1264 SetNonBendingCoor(extrapSegment->GetNonBendingCoor() +
1265 extrapSegment->GetNonBendingSlope() *
1266 (segment->GetHitForRec1()->GetZ() -
1267 (&(pMUON->Chamber(2 * station)))->Z()));
1269 NormalizedChi2WithSegment(extrapCorrSegment, maxSigma2Distance);
1270 if (chi2 < bestChi2) {
1271 // update best Chi2 and Segment if better found
1272 bestSegment = segment;
1277 // best segment found: add it to track candidate
1278 track->AddSegment(bestSegment);
1279 // set track parameters at these two TrakHit's
1280 track->SetTrackParamAtHit(track->GetNTrackHits() - 2, &(trackParam[0]));
1281 track->SetTrackParamAtHit(track->GetNTrackHits() - 1, &(trackParam[1]));
1282 if (fPrintLevel >= 10) {
1283 cout << "FollowTracks: track candidate(0..): " << trackIndex
1284 << " Added segment in station(0..): " << station << endl;
1285 track->RecursiveDump();
1289 // No best segment found:
1290 // Look for best compatible HitForRec in station:
1291 // should consider all possibilities ????
1292 // multiple scattering ???? do about like for extrapSegment !!!!
1293 extrapHit = new AliMUONHitForRec(); // empty hit
1294 extrapCorrHit = new AliMUONHitForRec(); // empty corrected hit
1297 if (fPrintLevel >= 10) {
1298 cout << "FollowTracks: track candidate(0..): " << trackIndex
1299 << " Segment not found, look for hit in station(0..): " << station
1302 // Loop over chambers of the station
1303 for (chInStation = 0; chInStation < 2; chInStation++) {
1304 // coordinates of extrapolated hit
1306 SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
1308 SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
1309 // resolutions from "extrapSegment"
1310 extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
1311 extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
1312 // same things for corrected hit
1313 // better to use copy constructor, after checking that it works properly !!!!
1315 SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
1317 SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
1318 extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
1319 extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
1320 // Loop over hits in the chamber
1321 ch = 2 * station + chInStation;
1322 for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
1323 iHit < fIndexOfFirstHitForRecPerChamber[ch] +
1324 fNHitsForRecPerChamber[ch];
1326 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
1327 // correction of corrected hit (fBendingCoor and fNonBendingCoor)
1328 // according to real Z value of "hit" and slopes of right "trackParam"
1330 SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor() +
1331 (&(trackParam[chInStation]))->GetBendingSlope() *
1333 (&(trackParam[chInStation]))->GetZ()));
1335 SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor() +
1336 (&(trackParam[chInStation]))->GetNonBendingSlope() *
1338 (&(trackParam[chInStation]))->GetZ()));
1339 // condition for hit not already in segment ????
1340 chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
1341 if (chi2 < bestChi2) {
1342 // update best Chi2 and HitForRec if better found
1345 chBestHit = chInStation;
1350 // best hit found: add it to track candidate
1351 track->AddHitForRec(bestHit);
1352 // set track parameters at this TrackHit
1353 track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
1354 &(trackParam[chBestHit]));
1355 if (fPrintLevel >= 10) {
1356 cout << "FollowTracks: track candidate(0..): " << trackIndex
1357 << " Added hit in station(0..): " << station << endl;
1358 track->RecursiveDump();
1362 // Remove current track candidate
1363 // and corresponding TrackHit's, ...
1365 delete extrapSegment;
1366 delete extrapCorrSegment;
1368 delete extrapCorrHit;
1369 break; // stop the search for this candidate:
1370 // exit from the loop over station
1373 delete extrapCorrHit;
1375 delete extrapSegment;
1376 delete extrapCorrSegment;
1377 // Sort track hits according to increasing Z
1378 track->GetTrackHitsPtr()->Sort();
1379 // Update track parameters at first track hit (smallest Z)
1380 trackParam1 = ((AliMUONTrackHit*)
1381 (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1382 bendingMomentum = 0.;
1383 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
1384 bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
1385 // Track removed if bendingMomentum not in window [min, max]
1386 if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
1388 break; // stop the search for this candidate:
1389 // exit from the loop over station
1392 // with multiple Coulomb scattering if all stations
1393 if (station == 0) track->SetFitMCS(1);
1394 // without multiple Coulomb scattering if not all stations
1395 else track->SetFitMCS(0);
1396 track->SetFitNParam(5); // with 5 parameters (momentum and position)
1397 track->SetFitStart(1); // from parameters at first hit
1399 Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
1400 if (numberOfDegFree > 0) {
1401 chi2Norm = track->GetFitFMin() / numberOfDegFree;
1405 // Track removed if normalized chi2 too high
1406 if (chi2Norm > fMaxChi2) {
1408 break; // stop the search for this candidate:
1409 // exit from the loop over station
1411 if (fPrintLevel >= 10) {
1412 cout << "FollowTracks: track candidate(0..): " << trackIndex
1413 << " after fit from station(0..): " << station << " to 4" << endl;
1414 track->RecursiveDump();
1416 // Track extrapolation to the vertex through the absorber (Branson)
1417 // after going through the first station
1419 trackParamVertex = *trackParam1;
1420 (&trackParamVertex)->ExtrapToVertex();
1421 track->SetTrackParamAtVertex(&trackParamVertex);
1422 if (fPrintLevel >= 1) {
1423 cout << "FollowTracks: track candidate(0..): " << trackIndex
1424 << " after extrapolation to vertex" << endl;
1425 track->RecursiveDump();
1428 } // for (station = 2;...
1429 // go really to next track
1432 // Compression of track array (necessary after Remove ????)
1433 fRecTracksPtr->Compress();
1437 //__________________________________________________________________________
1438 void AliMUONEventReconstructor::RemoveDoubleTracks(void)
1440 // To remove double tracks.
1441 // Tracks are considered identical
1442 // if they have at least half of their hits in common.
1443 // Among two identical tracks, one keeps the track with the larger number of hits
1444 // or, if these numbers are equal, the track with the minimum Chi2.
1445 AliMUONTrack *track1, *track2, *trackToRemove;
1446 Bool_t identicalTracks;
1447 Int_t hitsInCommon, nHits1, nHits2;
1448 identicalTracks = kTRUE;
1449 while (identicalTracks) {
1450 identicalTracks = kFALSE;
1451 // Loop over first track of the pair
1452 track1 = (AliMUONTrack*) fRecTracksPtr->First();
1453 while (track1 && (!identicalTracks)) {
1454 nHits1 = track1->GetNTrackHits();
1455 // Loop over second track of the pair
1456 track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1457 while (track2 && (!identicalTracks)) {
1458 nHits2 = track2->GetNTrackHits();
1459 // number of hits in common between two tracks
1460 hitsInCommon = track1->HitsInCommon(track2);
1461 // check for identical tracks
1462 if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
1463 identicalTracks = kTRUE;
1464 // decide which track to remove
1465 if (nHits1 > nHits2) trackToRemove = track2;
1466 else if (nHits1 < nHits2) trackToRemove = track1;
1467 else if ((track1->GetFitFMin()) < (track2->GetFitFMin()))
1468 trackToRemove = track2;
1469 else trackToRemove = track1;
1471 trackToRemove->Remove();
1473 track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1475 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1481 //__________________________________________________________________________
1482 void AliMUONEventReconstructor::EventDump(void)
1484 // Dump reconstructed event (track parameters at vertex and at first hit),
1485 // and the particle parameters
1487 AliMUONTrack *track;
1488 AliMUONTrackParam *trackParam, *trackParam1;
1490 Double_t bendingSlope, nonBendingSlope, pYZ;
1491 Double_t pX, pY, pZ, x, y, z, c;
1492 Int_t np, trackIndex, nTrackHits;
1494 if (fPrintLevel >= 1) cout << "****** enter EventDump ******" << endl;
1495 if (fPrintLevel >= 1) {
1496 cout << " Number of Reconstructed tracks :" << fNRecTracks << endl;
1498 fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1499 // Loop over reconstructed tracks
1500 for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1501 if (fTrackMethod != 1) continue; //AZ - skip the rest for now
1502 if (fPrintLevel >= 1)
1503 cout << " track number: " << trackIndex << endl;
1504 // function for each track for modularity ????
1505 track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1506 nTrackHits = track->GetNTrackHits();
1507 if (fPrintLevel >= 1)
1508 cout << " number of track hits: " << nTrackHits << endl;
1509 // track parameters at Vertex
1510 trackParam = track->GetTrackParamAtVertex();
1511 x = trackParam->GetNonBendingCoor();
1512 y = trackParam->GetBendingCoor();
1513 z = trackParam->GetZ();
1514 bendingSlope = trackParam->GetBendingSlope();
1515 nonBendingSlope = trackParam->GetNonBendingSlope();
1516 pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1517 pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1518 pX = pZ * nonBendingSlope;
1519 pY = pZ * bendingSlope;
1520 c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1521 if (fPrintLevel >= 1)
1522 printf(" track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1523 z, x, y, pX, pY, pZ, c);
1525 // track parameters at first hit
1526 trackParam1 = ((AliMUONTrackHit*)
1527 (track->GetTrackHitsPtr()->First()))->GetTrackParam();
1528 x = trackParam1->GetNonBendingCoor();
1529 y = trackParam1->GetBendingCoor();
1530 z = trackParam1->GetZ();
1531 bendingSlope = trackParam1->GetBendingSlope();
1532 nonBendingSlope = trackParam1->GetNonBendingSlope();
1533 pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1534 pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1535 pX = pZ * nonBendingSlope;
1536 pY = pZ * bendingSlope;
1537 c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1538 if (fPrintLevel >= 1)
1539 printf(" track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1540 z, x, y, pX, pY, pZ, c);
1542 // informations about generated particles
1543 np = gAlice->GetNtrack();
1544 printf(" **** number of generated particles: %d \n", np);
1546 for (Int_t iPart = 0; iPart < np; iPart++) {
1547 p = gAlice->Particle(iPart);
1548 printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1549 iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());
1554 void AliMUONEventReconstructor::FillEvent()
1556 // Create a new AliMUONRecoEvent, fill its track list, then add it as a
1557 // leaf in the Event branch of TreeRecoEvent tree
1558 cout << "Enter FillEvent() ...\n";
1561 fRecoEvent = new AliMUONRecoEvent();
1563 fRecoEvent->Clear();
1565 //save current directory
1566 TDirectory *current = gDirectory;
1567 if (!fTreeFile) fTreeFile = new TFile("tree_reco.root", "RECREATE");
1568 if (!fEventTree) fEventTree = new TTree("TreeRecoEvent", "MUON reconstructed events");
1569 //AZif (fRecoEvent->MakeDumpTracks(fRecTracksPtr)) {
1570 if (fRecoEvent->MakeDumpTracks(fMuons, fRecTracksPtr, this)) { //AZ
1571 if (fPrintLevel > 1) fRecoEvent->EventInfo();
1572 TBranch *branch = fEventTree->GetBranch("Event");
1573 if (!branch) branch = fEventTree->Branch("Event", "AliMUONRecoEvent", &fRecoEvent, 64000);
1574 branch->SetAutoDelete();
1579 // restore directory
1583 //__________________________________________________________________________
1584 void AliMUONEventReconstructor::MakeTrackCandidatesK(void)
1586 // To make initial tracks for Kalman filter from the list of segments
1588 AliMUONSegment *segment;
1589 AliMUONTrackK *trackK;
1591 if (fPrintLevel >= 1) cout << "enter MakeTrackCandidatesK" << endl;
1592 // Reset the TClonesArray of reconstructed tracks
1593 if (fRecTracksPtr) fRecTracksPtr->Delete();
1594 // Delete in order that the Track destructors are called,
1595 // hence the space for the TClonesArray of pointers to TrackHit's is freed
1598 AliMUONTrackK a(this, fHitsForRecPtr); // bad idea ???
1599 // Loop over stations(1...) 5 and 4
1600 for (istat=4; istat>=3; istat--) {
1601 // Loop over segments in the station
1602 for (iseg=0; iseg<fNSegments[istat]; iseg++) {
1603 // Transform segments to tracks and evaluate covariance matrix
1604 segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
1605 trackK = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrackK(segment);
1607 } // for (iseg=0;...)
1608 } // for (istat=4;...)
1612 //__________________________________________________________________________
1613 void AliMUONEventReconstructor::FollowTracksK(void)
1615 // Follow tracks using Kalman filter
1617 Int_t icand, ichamBeg, ichamEnd, chamBits;
1618 Double_t zDipole1, zDipole2;
1619 AliMUONTrackK *trackK;
1620 AliMUONHitForRec *hit;
1621 AliMUONRawCluster *clus;
1622 TClonesArray *rawclusters;
1624 clus = 0; rawclusters = 0;
1626 zDipole1 = GetSimpleBPosition() - GetSimpleBLength()/2;
1627 zDipole2 = zDipole1 + GetSimpleBLength();
1630 pMUON = (AliMUON*) gAlice->GetModule("MUON");
1631 for (Int_t i1=0; i1<fNHitsForRec; i1++) {
1632 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]);
1633 //if (hit->GetTHTrack() > 1 || hit->GetGeantSignal() == 0) continue;
1635 cout << " Hit #" << hit->GetChamberNumber() << " ";
1636 cout << hit->GetBendingCoor() << " ";
1637 cout << hit->GetNonBendingCoor() << " ";
1638 cout << hit->GetZ() << " ";
1639 cout << hit->GetGeantSignal() << " ";
1640 cout << hit->GetTHTrack() << endl;
1643 printf(" Hit # %d %10.4f %10.4f %10.4f",
1644 hit->GetChamberNumber(), hit->GetBendingCoor(),
1645 hit->GetNonBendingCoor(), hit->GetZ());
1646 if (fRecGeantHits) {
1648 printf(" %3d %3d \n", hit->GetGeantSignal(), hit->GetTHTrack());
1650 // from raw clusters
1651 rawclusters = pMUON->RawClustAddress(hit->GetChamberNumber());
1652 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->
1654 printf("%3d", clus->fTracks[1]-1);
1655 if (clus->fTracks[2] != 0) printf("%3d \n", clus->fTracks[2]-1);
1662 Int_t nSeeds = fNRecTracks; // starting number of seeds
1663 // Loop over track candidates
1664 while (icand < fNRecTracks-1) {
1666 trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1668 // Discard candidate which will produce the double track
1670 Ok = CheckCandidateK(icand,nSeeds);
1672 //trackK->SetRecover(-1); // mark candidate to be removed
1678 if (trackK->GetRecover() == 0) hit = (AliMUONHitForRec*)
1679 trackK->GetHitOnTrack()->Last(); // last hit
1680 else hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[1]; // 2'nd hit
1681 ichamBeg = hit->GetChamberNumber();
1683 // Check propagation direction
1684 if (trackK->GetTrackDir() > 0) {
1685 ichamEnd = 9; // forward propagation
1686 Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1688 ichamBeg = ichamEnd;
1689 ichamEnd = 6; // backward propagation
1690 // Change weight matrix and zero fChi2 for backpropagation
1691 trackK->StartBack();
1692 Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1693 ichamBeg = ichamEnd;
1697 if (trackK->GetBPFlag()) {
1699 ichamEnd = 6; // backward propagation
1700 // Change weight matrix and zero fChi2 for backpropagation
1701 trackK->StartBack();
1702 Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2);
1703 ichamBeg = ichamEnd;
1709 trackK->SetTrackDir(-1);
1710 trackK->SetBPFlag(kFALSE);
1711 Ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2);
1713 if (Ok) trackK->SetTrackQuality(0); // compute "track quality"
1714 else trackK->SetRecover(-1); // mark candidate to be removed
1716 // Majority 3 of 4 in first 2 stations
1718 for (Int_t i=0; i<trackK->GetNTrackHits(); i++) {
1719 hit = (AliMUONHitForRec*) (*trackK->GetHitOnTrack())[i];
1720 chamBits |= BIT(hit->GetChamberNumber()-1);
1722 //if (!((chamBits&3)==3 || (chamBits>>2&3)==3)) trackK->SetRecover(-1);
1723 //mark candidate to be removed
1726 for (Int_t i=0; i<fNRecTracks; i++) {
1727 trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1728 if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i);
1731 // Compress TClonesArray
1732 fRecTracksPtr->Compress();
1733 fNRecTracks = fRecTracksPtr->GetEntriesFast();
1737 //__________________________________________________________________________
1738 Bool_t AliMUONEventReconstructor::CheckCandidateK(Int_t icand, Int_t nSeeds)
1740 // Discards track candidate if it will produce the double track (having
1741 // the same seed segment hits as hits of a good track found before)
1742 AliMUONTrackK *track1, *track2;
1743 AliMUONHitForRec *hit1, *hit2, *hit;
1745 track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
1746 hit1 = (AliMUONHitForRec*) (*track1->GetHitOnTrack())[0]; // 1'st hit
1747 hit2 = (AliMUONHitForRec*) (*track1->GetHitOnTrack())[1]; // 2'nd hit
1749 for (Int_t i=0; i<icand; i++) {
1750 track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
1751 //if (track2->GetRecover() < 0) continue;
1752 if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
1754 if (track1->GetStartSegment() == track2->GetStartSegment()) {
1758 for (Int_t j=0; j<track2->GetNTrackHits(); j++) {
1759 hit = (AliMUONHitForRec*) (*track2->GetHitOnTrack())[j];
1760 if (hit == hit1 || hit == hit2) {
1762 if (nSame == 2) return kFALSE;
1764 } // for (Int_t j=0;
1766 } // for (Int_t i=0;
1770 //__________________________________________________________________________
1771 void AliMUONEventReconstructor::RemoveDoubleTracksK(void)
1773 // Removes double tracks (sharing more than half of their hits). Keeps
1774 // the track with higher quality
1775 AliMUONTrackK *track1, *track2, *trackToKill;
1777 // Sort tracks according to their quality
1778 fRecTracksPtr->Sort();
1780 // Loop over first track of the pair
1781 track1 = (AliMUONTrackK*) fRecTracksPtr->First();
1783 // Loop over second track of the pair
1784 track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1786 // Check whether or not to keep track2
1787 if (!track2->KeepTrack(track1)) {
1788 cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) <<
1789 " " << track2->GetTrackQuality() << endl;
1790 trackToKill = track2;
1791 track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1792 trackToKill->Kill();
1793 fRecTracksPtr->Compress();
1794 } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2);
1796 track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1);
1799 fNRecTracks = fRecTracksPtr->GetEntriesFast();
1800 cout << " Number of Kalman tracks: " << fNRecTracks << endl;
1803 //__________________________________________________________________________
1804 void AliMUONEventReconstructor::GoToVertex(void)
1806 // Propagates track to the vertex thru absorber
1807 // (using Branson correction for now)
1811 for (Int_t i=0; i<fNRecTracks; i++) {
1812 //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson();
1813 //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber
1814 ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2
1815 ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(); // with absorber