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 /// \class AliMUONTrackReconstructor
19 /// MUON track reconstructor using the original method
21 /// This class contains as data:
22 /// - the parameters for the track reconstruction
24 /// It contains as methods, among others:
25 /// - MakeTracks to build the tracks
29 #include <Riostream.h>
32 #include "AliMUONVTrackReconstructor.h"
33 #include "AliMUONTrackReconstructor.h"
34 #include "AliMUONData.h"
35 #include "AliMUONConstants.h"
36 #include "AliMUONRawCluster.h"
37 #include "AliMUONHitForRec.h"
38 #include "AliMUONSegment.h"
39 #include "AliMUONTrack.h"
40 #include "AliMUONTrackParam.h"
41 #include "AliMUONTrackExtrap.h"
43 #include <TVirtualFitter.h>
45 // Functions to be minimized with Minuit
46 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
47 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
49 void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
51 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param);
54 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
57 //************* Defaults parameters for reconstruction
58 const Double_t AliMUONTrackReconstructor::fgkDefaultMaxChi2 = 100.0;
60 TVirtualFitter* AliMUONTrackReconstructor::fgFitter = 0x0;
62 //__________________________________________________________________________
63 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONData* data)
64 : AliMUONVTrackReconstructor(data),
65 fMaxChi2(fgkDefaultMaxChi2)
67 /// Constructor for class AliMUONTrackReconstructor
69 // Memory allocation for the TClonesArray of reconstructed tracks
70 fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
73 //__________________________________________________________________________
74 AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void)
76 /// Destructor for class AliMUONTrackReconstructor
80 //__________________________________________________________________________
81 void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
83 /// To add to the list of hits for reconstruction all the raw clusters
84 TTree *TR = fMUONData->TreeR();
85 AliMUONHitForRec *hitForRec;
86 AliMUONRawCluster *clus;
87 Int_t iclus, nclus, nTRentries;
88 TClonesArray *rawclusters;
89 AliDebug(1,"Enter AddHitsForRecFromRawClusters");
91 fMUONData->SetTreeAddress("RC");
92 nTRentries = Int_t(TR->GetEntries());
93 if (nTRentries != 1) {
94 AliError(Form("nTRentries = %d not equal to 1 ",nTRentries));
97 fMUONData->GetRawClusters(); // only one entry
99 // Loop over tracking chambers
100 for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
101 // number of HitsForRec to 0 for the chamber
102 fNHitsForRecPerChamber[ch] = 0;
103 // index of first HitForRec for the chamber
104 if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
105 else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
106 rawclusters =fMUONData->RawClusters(ch);
107 nclus = (Int_t) (rawclusters->GetEntries());
108 // Loop over (cathode correlated) raw clusters
109 for (iclus = 0; iclus < nclus; iclus++) {
110 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
111 // new AliMUONHitForRec from raw cluster
112 // and increment number of AliMUONHitForRec's (total and in chamber)
113 hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
115 (fNHitsForRecPerChamber[ch])++;
116 // more information into HitForRec
117 hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
118 hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
119 // original raw cluster
120 hitForRec->SetChamberNumber(ch);
121 hitForRec->SetHitNumber(iclus);
122 // Z coordinate of the raw cluster (cm)
123 hitForRec->SetZ(clus->GetZ(0));
126 cout << "Chamber " << ch <<
127 " raw cluster " << iclus << " : " << endl;
129 cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
130 hitForRec->Print("full");
132 } // end of cluster loop
133 } // end of chamber loop
134 SortHitsForRecWithIncreasingChamber();
136 AliDebug(1,"End of AddHitsForRecFromRawClusters");
137 if (AliLog::GetGlobalDebugLevel() > 0) {
138 AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
139 for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
140 AliDebug(1, Form("Chamber(0...): %d",ch));
141 AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
142 AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
143 for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
144 hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
146 AliDebug(1, Form("HitForRec index(0...): %d",hit));
147 ((*fHitsForRecPtr)[hit])->Dump();
155 //__________________________________________________________________________
156 void AliMUONTrackReconstructor::MakeSegments(void)
158 /// To make the list of segments in all stations,
159 /// from the list of hits to be reconstructed
160 AliDebug(1,"Enter MakeSegments");
161 // Loop over stations
162 for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) MakeSegmentsPerStation(st);
165 cout << "end of MakeSegments" << endl;
166 for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
168 cout << "station " << st
169 << " has " << fNSegments[st] << " segments:"
171 for (Int_t seg = 0; seg < fNSegments[st]; seg++)
173 ((*fSegmentsPtr[st])[seg])->Print();
180 //__________________________________________________________________________
181 void AliMUONTrackReconstructor::MakeTracks(void)
183 /// To make the tracks from the list of segments and points in all stations
184 AliDebug(1,"Enter MakeTracks");
185 // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
186 MakeTrackCandidates();
187 // Follow tracks in stations(1..) 3, 2 and 1
189 // Remove double tracks
190 RemoveDoubleTracks();
191 UpdateHitForRecAtHit();
194 //__________________________________________________________________________
195 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
197 /// To make track candidates
198 /// with at least 3 aligned points in stations(1..) 4 and 5
199 /// (two Segment's or one Segment and one HitForRec)
200 Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
201 AliMUONSegment *begSegment;
202 AliDebug(1,"Enter MakeTrackCandidates");
203 // Loop over stations(1..) 5 and 4 for the beginning segment
204 for (begStation = 4; begStation > 2; begStation--) {
205 // Loop over segments in the beginning station
206 for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
207 // pointer to segment
208 begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
209 AliDebug(2,Form("Look for TrackCandidate's with Segment %d in Station(0..) %d", iBegSegment, begStation));
210 // Look for track candidates with two segments,
211 // "begSegment" and all compatible segments in other station.
212 // Only for beginning station(1..) 5
213 // because candidates with 2 segments have to looked for only once.
215 nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
216 // Look for track candidates with one segment and one point,
217 // "begSegment" and all compatible HitForRec's in other station.
218 // Only if "begSegment" does not belong already to a track candidate.
219 // Is that a too strong condition ????
220 if (!(begSegment->GetInTrack()))
221 nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
222 } // for (iBegSegment = 0;...
223 } // for (begStation = 4;...
227 //__________________________________________________________________________
228 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
230 /// To make track candidates with two segments in stations(1..) 4 and 5,
231 /// the first segment being pointed to by "BegSegment".
232 /// Returns the number of such track candidates.
233 Int_t endStation, iEndSegment, nbCan2Seg;
234 AliMUONSegment *endSegment;
235 AliMUONSegment *extrapSegment = NULL;
236 AliMUONTrack *recTrack;
238 AliDebug(1,"Enter MakeTrackCandidatesWithTwoSegments");
239 // Station for the end segment
240 endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
241 // multiple scattering factor corresponding to one chamber
243 GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
244 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
245 // linear extrapolation to end station
246 // number of candidates with 2 segments to 0
248 // Loop over segments in the end station
249 for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
250 // pointer to segment
251 endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
252 // test compatibility between current segment and "extrapSegment"
253 // 4 because 4 quantities in chi2
255 BegSegment->CreateSegmentFromLinearExtrapToStation(endSegment->GetZ(), mcsFactor);
257 NormalizedChi2WithSegment(extrapSegment,
258 fMaxSigma2Distance)) <= 4.0) {
259 // both segments compatible:
260 // make track candidate from "begSegment" and "endSegment"
261 AliDebug(2,Form("TrackCandidate with Segment %d in Station(0..) %d", iEndSegment, endStation));
262 // flag for both segments in one track:
263 // to be done in track constructor ????
264 BegSegment->SetInTrack(kTRUE);
265 endSegment->SetInTrack(kTRUE);
266 recTrack = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrack(BegSegment, endSegment);
267 // Set track parameters at vertex from last stations 4 & 5
268 CalcTrackParamAtVertex(recTrack);
270 if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
271 // increment number of track candidates with 2 segments
274 delete extrapSegment; // should not delete HitForRec's it points to !!!!
275 } // for (iEndSegment = 0;...
279 //__________________________________________________________________________
280 Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
282 /// To make track candidates with one segment and one point
283 /// in stations(1..) 4 and 5,
284 /// the segment being pointed to by "BegSegment".
285 Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
286 AliMUONHitForRec *extrapHitForRec= NULL;
287 AliMUONHitForRec *hit;
288 AliMUONTrack *recTrack;
290 AliDebug(1,"Enter MakeTrackCandidatesWithOneSegmentAndOnePoint");
291 // station for the end point
292 endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
293 // multiple scattering factor corresponding to one chamber
295 GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
296 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
297 // first and second chambers(0..) in the end station
298 ch1 = 2 * endStation;
300 // number of candidates to 0
302 // Loop over chambers of the end station
303 for (ch = ch2; ch >= ch1; ch--) {
304 // limits for the hit index in the loop
305 iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
306 iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
307 // Loop over HitForRec's in the chamber
308 for (iHit = iHitMin; iHit < iHitMax; iHit++) {
309 // pointer to HitForRec
310 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
311 // test compatibility between current HitForRec and "extrapHitForRec"
312 // 2 because 2 quantities in chi2
313 // linear extrapolation to chamber
315 BegSegment->CreateHitForRecFromLinearExtrapToChamber( hit->GetZ(), mcsFactor);
317 NormalizedChi2WithHitForRec(extrapHitForRec,
318 fMaxSigma2Distance)) <= 2.0) {
319 // both HitForRec's compatible:
320 // make track candidate from begSegment and current HitForRec
321 AliDebug(2, Form("TrackCandidate with HitForRec %d in Chamber(0..) %d", iHit, ch));
322 // flag for beginning segments in one track:
323 // to be done in track constructor ????
324 BegSegment->SetInTrack(kTRUE);
325 recTrack = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrack(BegSegment, hit);
326 // Set track parameters at vertex from last stations 4 & 5
327 CalcTrackParamAtVertex(recTrack);
328 // the right place to eliminate "double counting" ???? how ????
330 if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
331 // increment number of track candidates
334 delete extrapHitForRec;
335 } // for (iHit = iHitMin;...
336 } // for (ch = ch2;...
337 return nbCan1Seg1Hit;
340 //__________________________________________________________________________
341 void AliMUONTrackReconstructor::CalcTrackParamAtVertex(AliMUONTrack *Track) const
343 /// Set track parameters at vertex.
344 /// TrackHit's are assumed to be only in stations(1..) 4 and 5,
345 /// and sorted according to increasing Z.
346 /// Parameters are calculated from information in HitForRec's
347 /// of first and last TrackHit's.
348 AliMUONTrackParam *trackParamVertex = Track->GetTrackParamAtVertex(); // pointer to track parameters at vertex
349 // Pointer to HitForRec attached to first TrackParamAtHit
350 AliMUONHitForRec *firstHit = ((AliMUONTrackParam*) (Track->GetTrackParamAtHit()->First()))->GetHitForRecPtr();
351 // Pointer to HitForRec attached to last TrackParamAtHit
352 AliMUONHitForRec *lastHit = ((AliMUONTrackParam*) (Track->GetTrackParamAtHit()->Last()))->GetHitForRecPtr();
353 // Z difference between first and last hits
354 Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
355 // bending slope in stations(1..) 4 and 5
356 Double_t bendingSlope = (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
357 trackParamVertex->SetBendingSlope(bendingSlope);
359 Double_t impactParam = firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ();
360 // signed bending momentum
361 trackParamVertex->SetInverseBendingMomentum(1.0 / GetBendingMomentumFromImpactParam(impactParam));
362 // bending slope at vertex
363 trackParamVertex->SetBendingSlope(bendingSlope + impactParam / fSimpleBPosition);
365 trackParamVertex->SetNonBendingSlope((firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ);
366 // vertex coordinates at (0,0,0)
367 trackParamVertex->SetZ(0.0);
368 trackParamVertex->SetBendingCoor(0.0);
369 trackParamVertex->SetNonBendingCoor(0.0);
372 //__________________________________________________________________________
373 void AliMUONTrackReconstructor::FollowTracks(void)
375 /// Follow tracks in stations(1..) 3, 2 and 1
376 // too long: should be made more modular !!!!
377 AliMUONHitForRec *bestHit, *extrapHit, *hit;
378 AliMUONSegment *bestSegment, *extrapSegment, *segment;
379 AliMUONTrack *track, *nextTrack;
380 AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
381 // -1 to avoid compilation warnings
382 Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex;
383 Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
384 Double_t bendingMomentum, chi2Norm = 0.;
387 // local maxSigma2Distance, for easy increase in testing
388 maxSigma2Distance = fMaxSigma2Distance;
389 AliDebug(2,"Enter FollowTracks");
390 // Loop over track candidates
391 track = (AliMUONTrack*) fRecTracksPtr->First();
395 nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
396 AliDebug(2,Form("FollowTracks: track candidate(0..): %d", trackIndex));
397 // Fit track candidate from parameters at vertex
398 // -> with 3 parameters (X_vertex and Y_vertex are fixed)
399 // without multiple Coulomb scattering
401 if (AliLog::GetGlobalDebugLevel()> 2) {
402 cout << "FollowTracks: track candidate(0..): " << trackIndex
403 << " after fit in stations(0..) 3 and 4" << endl;
404 track->RecursiveDump();
406 // Loop over stations(1..) 3, 2 and 1
407 // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
408 // otherwise: majority coincidence 2 !!!!
409 for (station = 2; station >= 0; station--) {
410 // Track parameters at first track hit (smallest Z)
411 trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
412 // extrapolation to station
413 AliMUONTrackExtrap::ExtrapToStation(trackParam1, station, trackParam);
414 extrapSegment = new AliMUONSegment(); // empty segment
415 // multiple scattering factor corresponding to one chamber
416 // and momentum in bending plane (not total)
417 mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
418 mcsFactor = fChamberThicknessInX0 * mcsFactor * mcsFactor;
419 // Z difference from previous station
420 dZ1 = AliMUONConstants::DefaultChamberZ(2 * station) -
421 AliMUONConstants::DefaultChamberZ(2 * station + 2);
422 // Z difference between the two previous stations
423 dZ2 = AliMUONConstants::DefaultChamberZ(2 * station + 2) -
424 AliMUONConstants::DefaultChamberZ(2 * station + 4);
425 // Z difference between the two chambers in the previous station
426 dZ3 = AliMUONConstants::DefaultChamberZ(2 * station) -
427 AliMUONConstants::DefaultChamberZ(2 * station + 1);
428 extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
429 extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
430 extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
431 trackParam1->GetInverseBendingMomentum());
434 if (AliLog::GetGlobalDebugLevel() > 2) {
435 cout << "FollowTracks: track candidate(0..): " << trackIndex
436 << " Look for segment in station(0..): " << station << endl;
439 // Loop over segments in station
440 for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
441 // Look for best compatible Segment in station
442 // should consider all possibilities ????
443 // multiple scattering ????
444 // separation in 2 functions: Segment and HitForRec ????
445 segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
446 // correction of corrected segment (fBendingCoor and fNonBendingCoor)
447 // according to real Z value of "segment" and slopes of "extrapSegment"
448 AliMUONTrackExtrap::ExtrapToZ(&(trackParam[0]), segment->GetZ());
449 AliMUONTrackExtrap::ExtrapToZ(&(trackParam[1]), segment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
450 extrapSegment->SetBendingCoor((&(trackParam[0]))->GetBendingCoor());
451 extrapSegment->SetNonBendingCoor((&(trackParam[0]))->GetNonBendingCoor());
452 extrapSegment->SetBendingSlope((&(trackParam[0]))->GetBendingSlope());
453 extrapSegment->SetNonBendingSlope((&(trackParam[0]))->GetNonBendingSlope());
454 chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
455 if (chi2 < bestChi2) {
456 // update best Chi2 and Segment if better found
457 bestSegment = segment;
462 // best segment found: add it to track candidate
463 AliMUONTrackExtrap::ExtrapToZ(&(trackParam[0]), bestSegment->GetZ());
464 track->AddTrackParamAtHit(&(trackParam[0]),bestSegment->GetHitForRec1());
465 AliMUONTrackExtrap::ExtrapToZ(&(trackParam[1]), bestSegment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
466 track->AddTrackParamAtHit(&(trackParam[1]),bestSegment->GetHitForRec2());
467 AliDebug(3, Form("FollowTracks: track candidate(0..): %d Added segment in station(0..): %d", trackIndex, station));
468 if (AliLog::GetGlobalDebugLevel()>2) track->RecursiveDump();
470 // No best segment found:
471 // Look for best compatible HitForRec in station:
472 // should consider all possibilities ????
473 // multiple scattering ???? do about like for extrapSegment !!!!
474 extrapHit = new AliMUONHitForRec(); // empty hit
477 AliDebug(3, Form("FollowTracks: track candidate(0..): %d Segment not found, look for hit in station(0..): %d ",
478 trackIndex, station));
480 // Loop over chambers of the station
481 for (chInStation = 0; chInStation < 2; chInStation++) {
482 ch = 2 * station + chInStation;
483 for (iHit = fIndexOfFirstHitForRecPerChamber[ch]; iHit < fIndexOfFirstHitForRecPerChamber[ch]+fNHitsForRecPerChamber[ch]; iHit++) {
484 hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
485 // coordinates of extrapolated hit
486 AliMUONTrackExtrap::ExtrapToZ(&(trackParam[chInStation]), hit->GetZ());
487 extrapHit->SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
488 extrapHit->SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
489 // resolutions from "extrapSegment"
490 extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
491 extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
492 // Loop over hits in the chamber
493 // condition for hit not already in segment ????
494 chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
495 if (chi2 < bestChi2) {
496 // update best Chi2 and HitForRec if better found
499 chBestHit = chInStation;
504 // best hit found: add it to track candidate
505 AliMUONTrackExtrap::ExtrapToZ(&(trackParam[chBestHit]), bestHit->GetZ());
506 track->AddTrackParamAtHit(&(trackParam[chBestHit]),bestHit);
507 if (AliLog::GetGlobalDebugLevel() > 2) {
508 cout << "FollowTracks: track candidate(0..): " << trackIndex
509 << " Added hit in station(0..): " << station << endl;
510 track->RecursiveDump();
513 // Remove current track candidate
514 // and corresponding TrackHit's, ...
515 fRecTracksPtr->Remove(track);
517 delete extrapSegment;
519 break; // stop the search for this candidate:
520 // exit from the loop over station
524 delete extrapSegment;
525 // Sort TrackParamAtHit according to increasing Z
526 track->GetTrackParamAtHit()->Sort();
527 // Update track parameters at first track hit (smallest Z)
528 trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
529 bendingMomentum = 0.;
530 if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
531 bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
532 // Track removed if bendingMomentum not in window [min, max]
533 if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
534 fRecTracksPtr->Remove(track);
536 break; // stop the search for this candidate:
537 // exit from the loop over station
539 // Track fit from parameters at first hit
540 // -> with 5 parameters (momentum and position)
541 // with multiple Coulomb scattering if all stations
542 if (station == 0) Fit(track,1,1);
543 // without multiple Coulomb scattering if not all stations
545 Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
546 if (numberOfDegFree > 0) {
547 chi2Norm = track->GetFitFMin() / numberOfDegFree;
551 // Track removed if normalized chi2 too high
552 if (chi2Norm > fMaxChi2) {
553 fRecTracksPtr->Remove(track);
555 break; // stop the search for this candidate:
556 // exit from the loop over station
558 if (AliLog::GetGlobalDebugLevel() > 2) {
559 cout << "FollowTracks: track candidate(0..): " << trackIndex
560 << " after fit from station(0..): " << station << " to 4" << endl;
561 track->RecursiveDump();
563 // Track extrapolation to the vertex through the absorber (Branson)
564 // after going through the first station
566 trackParamVertex = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()));
567 AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, 0., 0., 0.);
568 track->SetTrackParamAtVertex(&trackParamVertex);
569 if (AliLog::GetGlobalDebugLevel() > 0) {
570 cout << "FollowTracks: track candidate(0..): " << trackIndex
571 << " after extrapolation to vertex" << endl;
572 track->RecursiveDump();
575 } // for (station = 2;...
576 // go really to next track
579 // Compression of track array (necessary after Remove)
580 fRecTracksPtr->Compress();
584 //__________________________________________________________________________
585 void AliMUONTrackReconstructor::Fit(AliMUONTrack *Track, Int_t FitStart, Int_t FitMCS)
587 /// Fit the track "Track",
588 /// with or without multiple Coulomb scattering according to "FitMCS",
589 /// starting, according to "FitStart",
590 /// with track parameters at vertex or at the first TrackHit.
592 if ((FitStart != 0) && (FitStart != 1)) {
593 cout << "ERROR in AliMUONTrackReconstructor::Fit(...)" << endl;
594 cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
597 if ((FitMCS != 0) && (FitMCS != 1)) {
598 cout << "ERROR in AliMUONTrackReconstructor::Fit(...)" << endl;
599 cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
603 Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
605 AliMUONTrackParam *trackParam;
606 // Check if Minuit is initialized...
607 fgFitter = TVirtualFitter::Fitter(Track,5);
608 fgFitter->Clear(); // necessary ???? probably yes
609 // how to reset the printout number at every fit ????
610 // is there any risk to leave it like that ????
611 // how to go faster ???? choice of Minuit parameters like EDM ????
612 // choice of function to be minimized according to fFitMCS
613 if (FitMCS == 0) fgFitter->SetFCN(TrackChi2);
614 else fgFitter->SetFCN(TrackChi2MCS);
615 // Switch off printout
617 fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
619 fgFitter->ExecuteCommand("SET NOW", arg, 0);
620 // Parameters according to "fFitStart"
621 // (should be a function to be used at every place where needed ????)
622 if (FitStart == 0) trackParam = Track->GetTrackParamAtVertex();
623 else trackParam = (AliMUONTrackParam*) (Track->GetTrackParamAtHit()->First());
624 // set first 3 Minuit parameters
625 // could be tried with no limits for the search (min=max=0) ????
626 fgFitter->SetParameter(0, "InvBenP",
627 trackParam->GetInverseBendingMomentum(),
629 fgFitter->SetParameter(1, "BenS",
630 trackParam->GetBendingSlope(),
632 fgFitter->SetParameter(2, "NonBenS",
633 trackParam->GetNonBendingSlope(),
636 // set last 2 Minuit parameters when we start from first track hit
637 // mandatory limits in Bending to avoid NaN values of parameters
638 fgFitter->SetParameter(3, "X",
639 trackParam->GetNonBendingCoor(),
640 0.03, -500.0, 500.0);
641 // mandatory limits in non Bending to avoid NaN values of parameters
642 fgFitter->SetParameter(4, "Y",
643 trackParam->GetBendingCoor(),
644 0.10, -500.0, 500.0);
646 // search without gradient calculation in the function
647 fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
649 fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
651 // fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
652 // get results into "invBenP", "benC", "nonBenC" ("x", "y")
653 fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
654 trackParam->SetInverseBendingMomentum(invBenP);
655 fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
656 trackParam->SetBendingSlope(benC);
657 fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
658 trackParam->SetNonBendingSlope(nonBenC);
660 fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
661 trackParam->SetNonBendingCoor(x);
662 fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
663 trackParam->SetBendingCoor(y);
665 // global result of the fit
666 Double_t fedm, errdef, fitFMin;
668 fgFitter->GetStats(fitFMin, fedm, errdef, npari, nparx);
669 Track->SetFitFMin(fitFMin);
672 //__________________________________________________________________________
673 void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
675 /// Return the "Chi2" to be minimized with Minuit for track fitting,
676 /// with "NParam" parameters
677 /// and their current values in array pointed to by "Param".
678 /// Assumes that the track hits are sorted according to increasing Z.
679 /// Track parameters at each TrackHit are updated accordingly.
680 /// Multiple Coulomb scattering is not taken into account
681 AliMUONTrack *trackBeingFitted;
682 AliMUONTrackParam param1;
683 AliMUONTrackParam* trackParamAtHit;
684 AliMUONHitForRec* hitForRec;
685 Chi2 = 0.0; // initialize Chi2
686 // copy of track parameters to be fitted
687 trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
688 // 3 parameters means fit track candidate from parameters at vertex (X_vertex and Y_vertex are fixed)
689 if (NParam == 3) param1 = *(trackBeingFitted->GetTrackParamAtVertex());
690 else param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
691 // Minuit parameters to be fitted into this copy
692 param1.SetInverseBendingMomentum(Param[0]);
693 param1.SetBendingSlope(Param[1]);
694 param1.SetNonBendingSlope(Param[2]);
696 param1.SetNonBendingCoor(Param[3]);
697 param1.SetBendingCoor(Param[4]);
699 // Follow track through all planes of track hits
700 trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First());
701 while (trackParamAtHit) {
702 hitForRec = trackParamAtHit->GetHitForRecPtr();
703 // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
704 AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ());
705 // update track parameters of the current hit
706 trackParamAtHit->SetTrackParam(param1);
708 // done hit per hit, with hit resolution,
709 // and not with point and angle like in "reco_muon.F" !!!!
710 // Needs to add multiple scattering contribution ????
711 Double_t dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
712 Double_t dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor();
713 Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
714 trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit));
718 //__________________________________________________________________________
719 void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
721 /// Return the "Chi2" to be minimized with Minuit for track fitting,
722 /// with "NParam" parameters
723 /// and their current values in array pointed to by "Param".
724 /// Assumes that the track hits are sorted according to increasing Z.
725 /// Track parameters at each TrackHit are updated accordingly.
726 /// Multiple Coulomb scattering is taken into account with covariance matrix.
727 AliMUONTrack *trackBeingFitted;
728 AliMUONTrackParam param1;
729 AliMUONTrackParam* trackParamAtHit;
730 AliMUONHitForRec* hitForRec;
731 Chi2 = 0.0; // initialize Chi2
732 // copy of track parameters to be fitted
733 trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
734 // 3 parameters means fit track candidate from parameters at vertex (X_vertex and Y_vertex are fixed)
735 if (NParam == 3) param1 = *(trackBeingFitted->GetTrackParamAtVertex());
736 else param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
737 // Minuit parameters to be fitted into this copy
738 param1.SetInverseBendingMomentum(Param[0]);
739 param1.SetBendingSlope(Param[1]);
740 param1.SetNonBendingSlope(Param[2]);
742 param1.SetNonBendingCoor(Param[3]);
743 param1.SetBendingCoor(Param[4]);
746 Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
748 AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3;
749 AliMUONHitForRec *hitForRec1, *hitForRec2;
750 Double_t hbc1, hbc2, pbc1, pbc2;
751 Double_t hnbc1, hnbc2, pnbc1, pnbc2;
752 Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
753 TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
754 TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
755 Double_t *msa2 = new Double_t[numberOfHit];
757 // Predicted coordinates and multiple scattering angles are first calculated
758 for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
759 trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
760 hitForRec = trackParamAtHit->GetHitForRecPtr();
761 // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
762 AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ());
763 // update track parameters of the current hit
764 trackParamAtHit->SetTrackParam(param1);
765 // square of multiple scattering angle at current hit, with one chamber
766 msa2[hitNumber] = MultipleScatteringAngle2(¶m1);
767 // correction for eventual missing hits or multiple hits in a chamber,
768 // according to the number of chambers
769 // between the current hit and the previous one
770 chCurrent = hitForRec->GetChamberNumber();
771 if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
775 // Calculates the covariance matrix
776 for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
777 trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
778 hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
779 z1 = hitForRec1->GetZ();
780 for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
781 trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
782 z2 = trackParamAtHit2->GetHitForRecPtr()->GetZ();
783 // initialization to 0 (diagonal plus upper triangular part)
784 (*covBending)(hitNumber2, hitNumber1) = 0.0;
785 // contribution from multiple scattering in bending plane:
786 // loop over upstream hits
787 for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {
788 trackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3));
789 z3 = trackParamAtHit3->GetHitForRecPtr()->GetZ();
790 (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]);
792 // equal contribution from multiple scattering in non bending plane
793 (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1);
794 if (hitNumber1 == hitNumber2) {
795 // Diagonal elements: add contribution from position measurements
797 (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + hitForRec1->GetBendingReso2();
798 // and in non bending plane
799 (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + hitForRec1->GetNonBendingReso2();
801 // Non diagonal elements: symmetrization
803 (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1);
804 // and non bending plane
805 (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1);
807 } // for (hitNumber2 = hitNumber1;...
808 } // for (hitNumber1 = 0;...
810 // Inversion of covariance matrices
811 // with "mnvertLocal", local "mnvert" function of Minuit.
812 // One cannot use directly "mnvert" since "TVirtualFitter" does not know it.
813 // One will have to replace this local function by the right inversion function
814 // from a specialized Root package for symmetric positive definite matrices,
815 // when available!!!!
817 mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
818 Int_t ifailNonBending;
819 mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
821 // It would be worth trying to calculate the inverse of the covariance matrix
822 // only once per fit, since it cannot change much in principle,
823 // and it would save a lot of computing time !!!!
826 if ((ifailBending == 0) && (ifailNonBending == 0)) {
827 // with Multiple Scattering if inversion correct
828 for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
829 trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
830 hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
831 hbc1 = hitForRec1->GetBendingCoor();
832 pbc1 = trackParamAtHit1->GetBendingCoor();
833 hnbc1 = hitForRec1->GetNonBendingCoor();
834 pnbc1 = trackParamAtHit1->GetNonBendingCoor();
835 for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
836 trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
837 hitForRec2 = trackParamAtHit2->GetHitForRecPtr();
838 hbc2 = hitForRec2->GetBendingCoor();
839 pbc2 = trackParamAtHit2->GetBendingCoor();
840 hnbc2 = hitForRec2->GetNonBendingCoor();
841 pnbc2 = trackParamAtHit2->GetNonBendingCoor();
842 Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) +
843 ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
847 // without Multiple Scattering if inversion impossible
848 for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
849 trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
850 hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
851 hbc1 = hitForRec1->GetBendingCoor();
852 pbc1 = trackParamAtHit1->GetBendingCoor();
853 hnbc1 = hitForRec1->GetNonBendingCoor();
854 pnbc1 = trackParamAtHit1->GetNonBendingCoor();
855 Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / hitForRec1->GetBendingReso2()) +
856 ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / hitForRec1->GetNonBendingReso2());
861 delete covNonBending;
865 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
867 /// Returns square of multiple Coulomb scattering angle
868 /// from TrackParamAtHit pointed to by "param"
869 Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
870 Double_t varMultipleScatteringAngle;
871 // Better implementation in AliMUONTrack class ????
872 slopeBending = param->GetBendingSlope();
873 slopeNonBending = param->GetNonBendingSlope();
874 // thickness in radiation length for the current track,
875 // taking local angle into account
876 radiationLength = AliMUONConstants::DefaultChamberThicknessInX0() *
877 TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending);
878 inverseBendingMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
879 inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
880 (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending);
881 varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
882 // The velocity is assumed to be 1 !!!!
883 varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle;
884 return varMultipleScatteringAngle;
887 //______________________________________________________________________________
888 void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
890 ///*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
891 ///*-* ==========================
892 ///*-* inverts a symmetric matrix. matrix is first scaled to
893 ///*-* have all ones on the diagonal (equivalent to change of units)
894 ///*-* but no pivoting is done since matrix is positive-definite.
895 ///*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
897 // taken from TMinuit package of Root (l>=n)
898 // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
899 // Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
900 Double_t * localVERTs = new Double_t[n];
901 Double_t * localVERTq = new Double_t[n];
902 Double_t * localVERTpp = new Double_t[n];
903 // fMaxint changed to localMaxint
904 Int_t localMaxint = n;
906 /* System generated locals */
909 /* Local variables */
911 Int_t i, j, k, kp1, km1;
913 /* Parameter adjustments */
919 if (n < 1) goto L100;
920 if (n > localMaxint) goto L100;
921 //*-*- scale matrix by sqrt of diag elements
922 for (i = 1; i <= n; ++i) {
924 if (si <= 0) goto L100;
925 localVERTs[i-1] = 1 / TMath::Sqrt(si);
927 for (i = 1; i <= n; ++i) {
928 for (j = 1; j <= n; ++j) {
929 a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
932 //*-*- . . . start main loop . . . .
933 for (i = 1; i <= n; ++i) {
935 //*-*- preparation for elimination step1
936 if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
938 localVERTpp[k-1] = 1;
942 if (km1 < 0) goto L100;
943 else if (km1 == 0) goto L50;
946 for (j = 1; j <= km1; ++j) {
947 localVERTpp[j-1] = a[j + k*l];
948 localVERTq[j-1] = a[j + k*l]*localVERTq[k-1];
952 if (k - n < 0) goto L51;
953 else if (k - n == 0) goto L60;
956 for (j = kp1; j <= n; ++j) {
957 localVERTpp[j-1] = a[k + j*l];
958 localVERTq[j-1] = -a[k + j*l]*localVERTq[k-1];
961 //*-*- elimination proper
963 for (j = 1; j <= n; ++j) {
964 for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
967 //*-*- elements of left diagonal and unscaling
968 for (j = 1; j <= n; ++j) {
969 for (k = 1; k <= j; ++k) {
970 a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
971 a[j + k*l] = a[k + j*l];
974 delete [] localVERTs;
975 delete [] localVERTq;
976 delete [] localVERTpp;
978 //*-*- failure return
980 delete [] localVERTs;
981 delete [] localVERTq;
982 delete [] localVERTpp;
986 //__________________________________________________________________________
987 void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
989 /// To remove double tracks.
990 /// Tracks are considered identical
991 /// if they have at least half of their hits in common.
992 /// Among two identical tracks, one keeps the track with the larger number of hits
993 /// or, if these numbers are equal, the track with the minimum Chi2.
994 AliMUONTrack *track1, *track2, *trackToRemove;
995 Int_t hitsInCommon, nHits1, nHits2;
996 Bool_t removedTrack1;
997 // Loop over first track of the pair
998 track1 = (AliMUONTrack*) fRecTracksPtr->First();
1000 removedTrack1 = kFALSE;
1001 nHits1 = track1->GetNTrackHits();
1002 // Loop over second track of the pair
1003 track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1005 nHits2 = track2->GetNTrackHits();
1006 // number of hits in common between two tracks
1007 hitsInCommon = track1->HitsInCommon(track2);
1008 // check for identical tracks
1009 if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
1010 // decide which track to remove
1011 if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() < track2->GetFitFMin()))) {
1012 // remove track2 and continue the second loop with the track next to track2
1013 trackToRemove = track2;
1014 track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1015 fRecTracksPtr->Remove(trackToRemove);
1017 fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
1019 // else remove track1 and continue the first loop with the track next to track1
1020 trackToRemove = track1;
1021 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1022 fRecTracksPtr->Remove(trackToRemove);
1024 fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
1025 removedTrack1 = kTRUE;
1028 } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
1030 if (removedTrack1) continue;
1031 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
1036 //__________________________________________________________________________
1037 void AliMUONTrackReconstructor::UpdateHitForRecAtHit()
1039 /// Set cluster parameters after track fitting. Fill fHitForRecAtHit of AliMUONTrack's
1040 AliMUONTrack *track;
1041 AliMUONTrackParam *trackParamAtHit;
1042 track = (AliMUONTrack*) fRecTracksPtr->First();
1044 trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
1045 while (trackParamAtHit) {
1046 track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
1047 trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
1049 track = (AliMUONTrack*) fRecTracksPtr->After(track);
1054 //__________________________________________________________________________
1055 void AliMUONTrackReconstructor::EventDump(void)
1057 /// Dump reconstructed event (track parameters at vertex and at first hit),
1058 /// and the particle parameters
1059 AliMUONTrack *track;
1060 AliMUONTrackParam *trackParam, *trackParam1;
1061 Double_t bendingSlope, nonBendingSlope, pYZ;
1062 Double_t pX, pY, pZ, x, y, z, c;
1063 Int_t trackIndex, nTrackHits;
1065 AliDebug(1,"****** enter EventDump ******");
1066 AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
1068 fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1069 // Loop over reconstructed tracks
1070 for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1071 AliDebug(1, Form("track number: %d", trackIndex));
1072 // function for each track for modularity ????
1073 track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1074 nTrackHits = track->GetNTrackHits();
1075 AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1076 // track parameters at Vertex
1077 trackParam = track->GetTrackParamAtVertex();
1078 x = trackParam->GetNonBendingCoor();
1079 y = trackParam->GetBendingCoor();
1080 z = trackParam->GetZ();
1081 bendingSlope = trackParam->GetBendingSlope();
1082 nonBendingSlope = trackParam->GetNonBendingSlope();
1083 pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1084 pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1085 pX = pZ * nonBendingSlope;
1086 pY = pZ * bendingSlope;
1087 c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1088 AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1089 z, x, y, pX, pY, pZ, c));
1091 // track parameters at first hit
1092 trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1093 x = trackParam1->GetNonBendingCoor();
1094 y = trackParam1->GetBendingCoor();
1095 z = trackParam1->GetZ();
1096 bendingSlope = trackParam1->GetBendingSlope();
1097 nonBendingSlope = trackParam1->GetNonBendingSlope();
1098 pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1099 pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1100 pX = pZ * nonBendingSlope;
1101 pY = pZ * bendingSlope;
1102 c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1103 AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1104 z, x, y, pX, pY, pZ, c));
1106 // informations about generated particles NO !!!!!!!!
1108 // for (Int_t iPart = 0; iPart < np; iPart++) {
1109 // p = gAlice->Particle(iPart);
1110 // printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1111 // iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());