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
28 #include "AliMUONTrackReconstructor.h"
29 #include "AliMUONData.h"
30 #include "AliMUONConstants.h"
31 #include "AliMUONRawCluster.h"
32 #include "AliMUONHitForRec.h"
33 #include "AliMUONObjectPair.h"
34 #include "AliMUONTrack.h"
35 #include "AliMUONTrackParam.h"
36 #include "AliMUONTrackExtrap.h"
41 #include <Riostream.h>
44 // Functions to be minimized with Minuit
45 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
46 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
48 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param);
51 ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
54 //************* Defaults parameters for reconstruction
55 const Double_t AliMUONTrackReconstructor::fgkMaxNormChi2 = 100.0;
56 const Bool_t AliMUONTrackReconstructor::fgkTrackAllTracks = kFALSE;
58 //__________________________________________________________________________
59 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONData* data)
60 : AliMUONVTrackReconstructor(data)
62 /// Constructor for class AliMUONTrackReconstructor
64 // Memory allocation for the TClonesArray of reconstructed tracks
65 fRecTracksPtr = new TClonesArray("AliMUONTrack", 10);
68 //__________________________________________________________________________
69 AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void)
71 /// Destructor for class AliMUONTrackReconstructor
75 //__________________________________________________________________________
76 void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
78 /// To add to the list of hits for reconstruction all the raw clusters
80 AliMUONHitForRec *hitForRec;
81 AliMUONRawCluster *clus;
83 TClonesArray *rawclusters;
84 AliDebug(1,"Enter AddHitsForRecFromRawClusters");
86 treeR = fMUONData->TreeR();
88 AliError("TreeR must be loaded");
92 fMUONData->SetTreeAddress("RC");
93 fMUONData->GetRawClusters(); // only one entry
95 // Loop over tracking chambers
96 for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
97 rawclusters =fMUONData->RawClusters(ch);
98 nclus = (Int_t) (rawclusters->GetEntries());
99 // Loop over (cathode correlated) raw clusters
100 for (iclus = 0; iclus < nclus; iclus++) {
101 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus);
102 // new AliMUONHitForRec from raw cluster
103 // and increment number of AliMUONHitForRec's (total and in chamber)
104 hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
106 // more information into HitForRec
107 hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
108 hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
109 // original raw cluster
110 hitForRec->SetChamberNumber(ch);
111 hitForRec->SetHitNumber(iclus);
112 // Z coordinate of the raw cluster (cm)
113 hitForRec->SetZ(clus->GetZ(0));
115 cout << "Chamber " << ch <<
116 " raw cluster " << iclus << " : " << endl;
118 cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl;
119 hitForRec->Print("full");
121 } // end of cluster loop
122 } // end of chamber loop
123 SortHitsForRecWithIncreasingChamber();
125 AliDebug(1,"End of AddHitsForRecFromRawClusters");
126 if (AliLog::GetGlobalDebugLevel() > 0) {
127 AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec));
128 for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
129 AliDebug(1, Form("Chamber(0...): %d",ch));
130 AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch]));
131 AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch]));
132 for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch];
133 hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch];
135 AliDebug(1, Form("HitForRec index(0...): %d",hit));
136 ((*fHitsForRecPtr)[hit])->Dump();
144 //__________________________________________________________________________
145 void AliMUONTrackReconstructor::MakeTracks(void)
147 /// To make the tracks from the list of segments and points in all stations
148 AliDebug(1,"Enter MakeTracks");
149 // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5
150 MakeTrackCandidates();
151 // Follow tracks in stations(1..) 3, 2 and 1
153 // Remove double tracks
154 RemoveDoubleTracks();
155 // Propagate tracks to the vertex through absorber
156 ExtrapTracksToVertex();
157 // Fill out the AliMUONTrack's
161 //__________________________________________________________________________
162 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
164 /// To make track candidates:
165 /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
166 /// Good candidates are made of at least three hitForRec's.
167 /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
168 TClonesArray *segments;
169 AliMUONObjectPair *segment;
170 AliMUONHitForRec *hitForRec1, *hitForRec2;
172 AliMUONTrackParam *trackParamAtFirstHit;
173 Int_t iCandidate = 0;
175 AliDebug(1,"Enter MakeTrackCandidates");
177 // Loop over stations(1..) 5 and 4 and make track candidates
178 for (Int_t istat=4; istat>=3; istat--) {
179 // Make segments in the station
180 segments = MakeSegmentsInStation(istat);
181 // Loop over segments
182 for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
183 AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
184 // Transform segments to tracks and put them at the end of fRecTracksPtr
185 segment = (AliMUONObjectPair*) ((*segments)[iseg]);
186 hitForRec1 = (AliMUONHitForRec*) segment->First();
187 hitForRec2 = (AliMUONHitForRec*) segment->Second();
188 track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(hitForRec1, hitForRec2);
190 // Add MCS effects in parameter covariances
191 trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
192 AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
193 // Printout for debuging
194 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
195 cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
196 trackParamAtFirstHit->GetCovariances()->Print();
198 // Look for compatible hitForRec(s) in the other station
199 FollowTrackInStation(track,7-istat);
201 // delete the array of segments
204 fRecTracksPtr->Compress(); // this is essential before checking tracks
206 // Keep all different tracks or only the best ones as required
207 if (fgkTrackAllTracks) RemoveIdenticalTracks();
208 else RemoveDoubleTracks();
210 AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
214 //__________________________________________________________________________
215 void AliMUONTrackReconstructor::RemoveIdenticalTracks(void)
217 /// To remove identical tracks:
218 /// Tracks are considered identical if they have all their hits in common.
219 /// One keeps the track with the larger number of hits if need be
220 AliMUONTrack *track1, *track2, *trackToRemove;
221 Int_t hitsInCommon, nHits1, nHits2;
222 Bool_t removedTrack1;
223 // Loop over first track of the pair
224 track1 = (AliMUONTrack*) fRecTracksPtr->First();
226 removedTrack1 = kFALSE;
227 nHits1 = track1->GetNTrackHits();
228 // Loop over second track of the pair
229 track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
231 nHits2 = track2->GetNTrackHits();
232 // number of hits in common between two tracks
233 hitsInCommon = track1->HitsInCommon(track2);
234 // check for identical tracks
235 if ((hitsInCommon == nHits1) || (hitsInCommon == nHits2)) {
236 // decide which track to remove
237 if (nHits2 > nHits1) {
238 // remove track1 and continue the first loop with the track next to track1
239 trackToRemove = track1;
240 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
241 fRecTracksPtr->Remove(trackToRemove);
242 fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
244 removedTrack1 = kTRUE;
247 // remove track2 and continue the second loop with the track next to track2
248 trackToRemove = track2;
249 track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
250 fRecTracksPtr->Remove(trackToRemove);
251 fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
254 } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
256 if (removedTrack1) continue;
257 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
262 //__________________________________________________________________________
263 void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
265 /// To remove double tracks:
266 /// Tracks are considered identical if more than half of the hits of the track
267 /// which has the smaller number of hits are in common with the other track.
268 /// Among two identical tracks, one keeps the track with the larger number of hits
269 /// or, if these numbers are equal, the track with the minimum chi2.
270 AliMUONTrack *track1, *track2, *trackToRemove;
271 Int_t hitsInCommon, nHits1, nHits2;
272 Bool_t removedTrack1;
273 // Loop over first track of the pair
274 track1 = (AliMUONTrack*) fRecTracksPtr->First();
276 removedTrack1 = kFALSE;
277 nHits1 = track1->GetNTrackHits();
278 // Loop over second track of the pair
279 track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
281 nHits2 = track2->GetNTrackHits();
282 // number of hits in common between two tracks
283 hitsInCommon = track1->HitsInCommon(track2);
284 // check for identical tracks
285 if (((nHits1 < nHits2) && (2 * hitsInCommon > nHits1)) || (2 * hitsInCommon > nHits2)) {
286 // decide which track to remove
287 if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() <= track2->GetFitFMin()))) {
288 // remove track2 and continue the second loop with the track next to track2
289 trackToRemove = track2;
290 track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
291 fRecTracksPtr->Remove(trackToRemove);
292 fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
295 // else remove track1 and continue the first loop with the track next to track1
296 trackToRemove = track1;
297 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
298 fRecTracksPtr->Remove(trackToRemove);
299 fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
301 removedTrack1 = kTRUE;
304 } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
306 if (removedTrack1) continue;
307 track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
312 //__________________________________________________________________________
313 void AliMUONTrackReconstructor::FollowTracks(void)
315 /// Follow tracks in stations(1..) 3, 2 and 1
316 AliDebug(1,"Enter FollowTracks");
318 AliMUONTrack *track, *nextTrack;
319 AliMUONTrackParam *trackParamAtFirstHit;
320 Double_t numberOfDegFree, chi2Norm;
321 Int_t currentNRecTracks;
323 for (Int_t station = 2; station >= 0; station--) {
324 // Save the actual number of reconstructed track in case of
325 // tracks are added or suppressed during the tracking procedure
326 // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
327 currentNRecTracks = fNRecTracks;
328 for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) {
329 AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
330 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
332 // Do not take into account the multiple scattering to speed up the fit
333 // Calculate the track parameter covariance matrix
334 // If "station" is station(1..) 3 then use the vertex to better constrain the fit
336 SetVertexForFit(track);
337 track->SetFitWithVertex(kTRUE);
338 } else track->SetFitWithVertex(kFALSE);
339 Fit(track,kFALSE, kTRUE);
340 // Remove the track if the normalized chi2 is too high
341 numberOfDegFree = (2. * track->GetNTrackHits() - 5.);
342 if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree;
343 else chi2Norm = 1.e10;
344 if (chi2Norm > fgkMaxNormChi2) {
345 fRecTracksPtr->Remove(track);
349 // Add MCS effects in parameter covariances
350 trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
351 AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
352 // Printout for debuging
353 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
354 cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
355 trackParamAtFirstHit->GetCovariances()->Print();
357 // Look for compatible hitForRec in station(0..) "station"
358 FollowTrackInStation(track,station);
360 // Compress fRecTracksPtr for the next step
361 fRecTracksPtr->Compress();
362 // Keep only the best tracks if required
363 if (!fgkTrackAllTracks) RemoveDoubleTracks();
366 // Last fit of track candidates with all station
367 // Take into account the multiple scattering and remove bad tracks
368 Int_t trackIndex = -1;
369 track = (AliMUONTrack*) fRecTracksPtr->First();
372 nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
373 track->SetFitWithVertex(kFALSE); // just to be sure
374 Fit(track,kTRUE, kFALSE);
375 // Printout for debuging
376 if (AliLog::GetGlobalDebugLevel() >= 3) {
377 cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
378 track->RecursiveDump();
380 // Remove the track if the normalized chi2 is too high
381 numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
382 if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree;
383 else chi2Norm = 1.e10;
384 if (chi2Norm > fgkMaxNormChi2) {
385 fRecTracksPtr->Remove(track);
390 fRecTracksPtr->Compress();
394 //__________________________________________________________________________
395 void AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation)
397 /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s)
398 /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
399 /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
400 /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
401 /// Remove the obsolete "trackCandidate" at the end.
402 /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority.
403 AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
405 Int_t ch1 = 2*nextStation;
406 Int_t ch2 = 2*nextStation+1;
407 Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2);
408 Double_t chi2WithOneHitForRec = 1.e10;
409 Double_t chi2WithTwoHitForRec = 1.e10;
410 Double_t maxChi2WithOneHitForRec = 2.*fgkMaxNormChi2; // 2 because 2 quantities in chi2
411 Double_t maxChi2WithTwoHitForRec = 4.*fgkMaxNormChi2; // 4 because 4 quantities in chi2
412 Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec;
413 Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec;
414 AliMUONTrack *newTrack = 0x0;
415 AliMUONHitForRec *hitForRecCh1, *hitForRecCh2;
416 AliMUONHitForRec *bestHitForRec1 = 0x0, *bestHitForRec2 = 0x0;
417 Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]];
418 for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE;
420 //Extrapolate trackCandidate to chamber "ch2" to save computing time in the next steps
421 AliMUONTrackParam *extrapTrackParamPtr = trackCandidate->GetExtrapTrackParam();
422 *extrapTrackParamPtr = *((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First()));
423 AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, zCh2);
424 AliMUONTrackParam extrapTrackParamSave(*extrapTrackParamPtr);
426 // Printout for debuging
427 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
428 TMatrixD* paramCovForDebug = extrapTrackParamPtr->GetCovariances();
429 cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl;
430 paramCovForDebug->Print();
433 // look for candidates in chamber 2
434 // Printout for debuging
435 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
436 cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl;
438 for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) {
439 hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2);
440 // extrapolate track parameters and covariances only once for this hit
441 AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, hitForRecCh2->GetZ());
442 chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh2);
443 // if good chi2 then try to attach a hitForRec in the other chamber too
444 if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
445 // Printout for debuging
446 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
447 cout << "FollowTrackInStation: look for second hits in chamber(1..): " << ch1+1 << endl;
449 Bool_t foundSecondHit = kFALSE;
450 for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
451 hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
452 chi2WithTwoHitForRec = trackCandidate->TryTwoHitForRec(hitForRecCh2, hitForRecCh1); // order hits like that to save computing time
453 // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate"
454 if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) {
455 foundSecondHit = kTRUE;
456 if (fgkTrackAllTracks) {
457 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
458 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
460 AliMUONTrackParam trackParam1(extrapTrackParamSave);
461 AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh1->GetZ());
462 newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh1);
463 AliMUONTrackParam trackParam2(extrapTrackParamSave);
464 AliMUONTrackExtrap::ExtrapToZ(&trackParam2, hitForRecCh2->GetZ());
465 newTrack->AddTrackParamAtHit(&trackParam2,hitForRecCh2);
466 // Sort TrackParamAtHit according to increasing Z
467 newTrack->GetTrackParamAtHit()->Sort();
468 // Update the chi2 of the new track
469 if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithTwoHitForRec);
470 else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithTwoHitForRec);
471 // Tag hitForRecCh1 as used
472 hitForRecCh1Used[hit1] = kTRUE;
473 // Printout for debuging
474 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
475 cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1
476 << " (Chi2 = " << chi2WithTwoHitForRec << ")" << endl;
477 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
479 } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) {
480 // keep track of the best couple of hits
481 bestChi2WithTwoHitForRec = chi2WithTwoHitForRec;
482 bestHitForRec1 = hitForRecCh1;
483 bestHitForRec2 = hitForRecCh2;
487 // if no hitForRecCh1 found then consider to add hitForRecCh2 only
488 if (!foundSecondHit) {
489 if (fgkTrackAllTracks) {
490 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
491 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
493 AliMUONTrackParam trackParam1(extrapTrackParamSave);
494 AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh2->GetZ());
495 newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh2);
496 // Sort TrackParamAtHit according to increasing Z
497 newTrack->GetTrackParamAtHit()->Sort();
498 // Update the chi2 of the new track
499 if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
500 else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
501 // Printout for debuging
502 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
503 cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1
504 << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
505 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
507 } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
508 // keep track of the best single hitForRec except if a couple
509 // of hits has already been found (i.e. bestHitForRec2!=0x0)
510 bestChi2WithOneHitForRec = chi2WithOneHitForRec;
511 bestHitForRec1 = hitForRecCh2;
515 // reset the extrapolated track parameter for next step
516 trackCandidate->SetExtrapTrackParam(&extrapTrackParamSave);
519 // look for candidates in chamber 1 not already attached to a track
520 // if we want to keep all possible tracks or if no good couple of hitForRec has been found
521 if (fgkTrackAllTracks || !bestHitForRec2) {
522 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
523 cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl;
525 for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
526 hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
527 if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used
528 chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh1);
529 // if good chi2 then create a new track by adding the good hitForRec in "ch1" to the "trackCandidate"
530 // We do not try to attach a hitForRec in the other chamber too since it has already been done above
531 if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
532 if (fgkTrackAllTracks) {
533 // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
534 newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
536 AliMUONTrackParam trackParam1(extrapTrackParamSave);
537 AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh1->GetZ());
538 newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh1);
539 // Sort TrackParamAtHit according to increasing Z
540 newTrack->GetTrackParamAtHit()->Sort();
541 // Update the chi2 of the new track
542 if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
543 else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
544 // Printout for debuging
545 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
546 cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1
547 << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
548 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
550 } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
551 // keep track of the best single hitForRec except if a couple
552 // of hits has already been found (i.e. bestHitForRec1!=0x0)
553 bestChi2WithOneHitForRec = chi2WithOneHitForRec;
554 bestHitForRec1 = hitForRecCh1;
560 // fill out the best track if required else clean up the fRecTracksPtr array
561 if (!fgkTrackAllTracks && bestHitForRec1) {
562 AliMUONTrackParam trackParam1(extrapTrackParamSave);
563 AliMUONTrackExtrap::ExtrapToZ(&trackParam1, bestHitForRec1->GetZ());
564 trackCandidate->AddTrackParamAtHit(&trackParam1,bestHitForRec1);
565 if (bestHitForRec2) {
566 AliMUONTrackParam trackParam2(extrapTrackParamSave);
567 AliMUONTrackExtrap::ExtrapToZ(&trackParam2, bestHitForRec2->GetZ());
568 trackCandidate->AddTrackParamAtHit(&trackParam2,bestHitForRec2);
569 // Sort TrackParamAtHit according to increasing Z
570 trackCandidate->GetTrackParamAtHit()->Sort();
571 // Update the chi2 of the new track
572 if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithTwoHitForRec);
573 else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithTwoHitForRec);
574 // Printout for debuging
575 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
576 cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1
577 << " (Chi2 = " << bestChi2WithTwoHitForRec << ")" << endl;
578 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
581 // Sort TrackParamAtHit according to increasing Z
582 trackCandidate->GetTrackParamAtHit()->Sort();
583 // Update the chi2 of the new track
584 if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithOneHitForRec);
585 else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithOneHitForRec);
586 // Printout for debuging
587 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
588 cout << "FollowTrackInStation: added the best hit in station(1..): " << nextStation+1
589 << " (Chi2 = " << bestChi2WithOneHitForRec << ")" << endl;
590 if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
594 fRecTracksPtr->Remove(trackCandidate); // obsolete track
600 //__________________________________________________________________________
601 void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate)
603 /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate"
604 /// Compute the vertex resolution from natural vertex dispersion and
605 /// multiple scattering effets according to trackCandidate path in absorber
606 /// It is necessary to account for multiple scattering effects here instead of during the fit of
607 /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
608 AliDebug(1,"Enter SetVertexForFit");
610 Double_t nonBendingReso2 = fNonBendingVertexDispersion * fNonBendingVertexDispersion;
611 Double_t bendingReso2 = fBendingVertexDispersion * fBendingVertexDispersion;
612 // add multiple scattering effets
613 AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First())));
614 paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
615 AliMUONTrackExtrap::ExtrapToVertexUncorrected(¶mAtVertex,0.);
616 TMatrixD* paramCov = paramAtVertex.GetCovariances();
617 nonBendingReso2 += (*paramCov)(0,0);
618 bendingReso2 += (*paramCov)(2,2);
620 AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default
621 vertex.SetNonBendingReso2(nonBendingReso2);
622 vertex.SetBendingReso2(bendingReso2);
623 trackCandidate->SetVertex(&vertex);
626 //__________________________________________________________________________
627 void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov)
629 /// Fit the track "track" with or without multiple Coulomb scattering according to "includeMCS".
631 Double_t benC, errorParam, invBenP, nonBenC, x, y;
632 AliMUONTrackParam *trackParam;
633 Double_t arg[1], fedm, errdef, fitFMin;
635 Int_t status, covStatus;
637 // Clear MINUIT parameters
639 // Give the fitted track to MINUIT
640 gMinuit->SetObjectFit(track);
641 if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
642 // Define print level
644 gMinuit->mnexcm("SET PRI", arg, 1, status);
645 // Print covariance matrix
646 gMinuit->mnexcm("SHO COV", arg, 0, status);
649 gMinuit->mnexcm("SET PRI", arg, 1, status);
652 gMinuit->mnexcm("SET NOW", arg, 0, status);
655 //gMinuit->mnexcm("SET STR", arg, 1, status);
657 // Switch between available FCN according to "includeMCS"
658 if (includeMCS) gMinuit->SetFCN(TrackChi2MCS);
659 else gMinuit->SetFCN(TrackChi2);
661 // Set fitted parameters (!! The order is very important for the covariance matrix !!)
662 trackParam = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
663 // could be tried with no limits for the search (min=max=0) ????
664 // mandatory limits in non Bending to avoid NaN values of parameters
665 gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
666 gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status);
667 // mandatory limits in Bending to avoid NaN values of parameters
668 gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
669 gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status);
670 gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status);
673 gMinuit->mnexcm("MIGRAD", arg, 0, status);
675 // Calculate the covariance matrix more accurately if required
676 if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
678 // get results into "invBenP", "benC", "nonBenC" ("x", "y")
679 gMinuit->GetParameter(0, x, errorParam);
680 trackParam->SetNonBendingCoor(x);
681 gMinuit->GetParameter(1, nonBenC, errorParam);
682 trackParam->SetNonBendingSlope(nonBenC);
683 gMinuit->GetParameter(2, y, errorParam);
684 trackParam->SetBendingCoor(y);
685 gMinuit->GetParameter(3, benC, errorParam);
686 trackParam->SetBendingSlope(benC);
687 gMinuit->GetParameter(4, invBenP, errorParam);
688 trackParam->SetInverseBendingMomentum(invBenP);
690 // global result of the fit
691 gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus);
692 track->SetFitFMin(fitFMin);
694 // Get the covariance matrix if required
696 // Covariance matrix according to HESSE status
697 // If problem then keep only the diagonal terms (variances)
698 Double_t matrix[5][5];
699 gMinuit->mnemat(&matrix[0][0],5);
700 if (covStatus == 3) trackParam->SetCovariances(matrix);
701 else trackParam->SetVariances(matrix);
702 } else *(trackParam->GetCovariances()) = 0;
706 //__________________________________________________________________________
707 void TrackChi2(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
709 /// Return the "Chi2" to be minimized with Minuit for track fitting,
710 /// with "NParam" parameters
711 /// and their current values in array pointed to by "Param".
712 /// Assumes that the track hits are sorted according to increasing Z.
713 /// Track parameters at each TrackHit are updated accordingly.
714 /// Multiple Coulomb scattering is not taken into account
716 AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
717 // AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
718 AliMUONTrackParam param1;
719 AliMUONTrackParam* trackParamAtHit;
720 AliMUONHitForRec* hitForRec;
721 Chi2 = 0.0; // initialize Chi2
724 // copy of track parameters to be fitted
725 param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
726 param1.SetNonBendingCoor(Param[0]);
727 param1.SetNonBendingSlope(Param[1]);
728 param1.SetBendingCoor(Param[2]);
729 param1.SetBendingSlope(Param[3]);
730 param1.SetInverseBendingMomentum(Param[4]);
732 // Take the vertex into account in the fit if required
733 if (trackBeingFitted->GetFitWithVertex()) {
734 AliMUONTrackParam paramAtVertex(param1);
735 AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.);
736 AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
738 cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl;
741 dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
742 dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
743 Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
746 // Follow track through all planes of track hits
747 trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First());
748 while (trackParamAtHit) {
749 hitForRec = trackParamAtHit->GetHitForRecPtr();
750 // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
751 AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ());
752 // update track parameters of the current hit
753 trackParamAtHit->SetTrackParam(param1);
755 // done hit per hit, with hit resolution,
756 // and not with point and angle like in "reco_muon.F" !!!!
757 dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
758 dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor();
759 Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
760 trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit));
764 //__________________________________________________________________________
765 void TrackChi2MCS(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
767 /// Return the "Chi2" to be minimized with Minuit for track fitting,
768 /// with "NParam" parameters
769 /// and their current values in array pointed to by "Param".
770 /// Assumes that the track hits are sorted according to increasing Z.
771 /// Track parameters at each TrackHit are updated accordingly.
772 /// Multiple Coulomb scattering is taken into account with covariance matrix.
774 AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
775 // AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
776 AliMUONTrackParam param1;
777 AliMUONTrackParam* trackParamAtHit;
778 AliMUONHitForRec* hitForRec;
779 Chi2 = 0.0; // initialize Chi2
780 Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
782 AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3;
783 AliMUONHitForRec *hitForRec1, *hitForRec2;
784 Double_t hbc1, hbc2, pbc1, pbc2;
785 Double_t hnbc1, hnbc2, pnbc1, pnbc2;
786 Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
787 TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
788 TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
789 Double_t *msa2 = new Double_t[numberOfHit];
791 // copy of track parameters to be fitted
792 param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
793 param1.SetNonBendingCoor(Param[0]);
794 param1.SetNonBendingSlope(Param[1]);
795 param1.SetBendingCoor(Param[2]);
796 param1.SetBendingSlope(Param[3]);
797 param1.SetInverseBendingMomentum(Param[4]);
799 // Take the vertex into account in the fit if required
800 if (trackBeingFitted->GetFitWithVertex()) {
801 AliMUONTrackParam paramAtVertex(param1);
802 AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.);
803 AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
805 cout<<"Error in TrackChi2MCS: Want to use the vertex in tracking but it has not been created!!"<<endl;
808 Double_t dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
809 Double_t dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
810 Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
813 // Predicted coordinates and multiple scattering angles are first calculated
814 for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
815 trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
816 hitForRec = trackParamAtHit->GetHitForRecPtr();
817 // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit
818 AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ());
819 // update track parameters of the current hit
820 trackParamAtHit->SetTrackParam(param1);
821 // square of multiple scattering angle at current hit, with one chamber
822 msa2[hitNumber] = MultipleScatteringAngle2(¶m1);
823 // correction for eventual missing hits or multiple hits in a chamber,
824 // according to the number of chambers
825 // between the current hit and the previous one
826 chCurrent = hitForRec->GetChamberNumber();
827 if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev);
831 // Calculates the covariance matrix
832 for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {
833 trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
834 hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
835 z1 = hitForRec1->GetZ();
836 for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
837 trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
838 z2 = trackParamAtHit2->GetHitForRecPtr()->GetZ();
839 // initialization to 0 (diagonal plus upper triangular part)
840 (*covBending)(hitNumber2, hitNumber1) = 0.0;
841 // contribution from multiple scattering in bending plane:
842 // loop over upstream hits
843 for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {
844 trackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3));
845 z3 = trackParamAtHit3->GetHitForRecPtr()->GetZ();
846 (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]);
848 // equal contribution from multiple scattering in non bending plane
849 (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1);
850 if (hitNumber1 == hitNumber2) {
851 // Diagonal elements: add contribution from position measurements
853 (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + hitForRec1->GetBendingReso2();
854 // and in non bending plane
855 (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + hitForRec1->GetNonBendingReso2();
857 // Non diagonal elements: symmetrization
859 (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1);
860 // and non bending plane
861 (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1);
863 } // for (hitNumber2 = hitNumber1;...
864 } // for (hitNumber1 = 0;...
866 // Inversion of covariance matrices
868 gMinuit->mnvert(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
869 Int_t ifailNonBending;
870 gMinuit->mnvert(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
872 // It would be worth trying to calculate the inverse of the covariance matrix
873 // only once per fit, since it cannot change much in principle,
874 // and it would save a lot of computing time !!!!
877 if ((ifailBending == 0) && (ifailNonBending == 0)) {
878 // with Multiple Scattering if inversion correct
879 for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
880 trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
881 hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
882 hbc1 = hitForRec1->GetBendingCoor();
883 pbc1 = trackParamAtHit1->GetBendingCoor();
884 hnbc1 = hitForRec1->GetNonBendingCoor();
885 pnbc1 = trackParamAtHit1->GetNonBendingCoor();
886 for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) {
887 trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2));
888 hitForRec2 = trackParamAtHit2->GetHitForRecPtr();
889 hbc2 = hitForRec2->GetBendingCoor();
890 pbc2 = trackParamAtHit2->GetBendingCoor();
891 hnbc2 = hitForRec2->GetNonBendingCoor();
892 pnbc2 = trackParamAtHit2->GetNonBendingCoor();
893 Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) +
894 ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2));
898 // without Multiple Scattering if inversion impossible
899 for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) {
900 trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1));
901 hitForRec1 = trackParamAtHit1->GetHitForRecPtr();
902 hbc1 = hitForRec1->GetBendingCoor();
903 pbc1 = trackParamAtHit1->GetBendingCoor();
904 hnbc1 = hitForRec1->GetNonBendingCoor();
905 pnbc1 = trackParamAtHit1->GetNonBendingCoor();
906 Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / hitForRec1->GetBendingReso2()) +
907 ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / hitForRec1->GetNonBendingReso2());
912 delete covNonBending;
916 //__________________________________________________________________________
917 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
919 /// Returns square of multiple Coulomb scattering angle
920 /// from TrackParamAtHit pointed to by "param"
921 Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
922 Double_t varMultipleScatteringAngle;
923 slopeBending = param->GetBendingSlope();
924 slopeNonBending = param->GetNonBendingSlope();
925 // thickness in radiation length for the current track,
926 // taking local angle into account
927 radiationLength = AliMUONConstants::ChamberThicknessInX0() *
928 TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending);
929 inverseBendingMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
930 inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
931 (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending);
932 varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
933 // The velocity is assumed to be 1 !!!!
934 varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle;
935 return varMultipleScatteringAngle;
938 //__________________________________________________________________________
939 void AliMUONTrackReconstructor::ExtrapTracksToVertex()
941 /// propagate tracks to the vertex through the absorber (Branson)
943 AliMUONTrackParam trackParamVertex;
944 AliMUONHitForRec *vertex;
947 for (Int_t iRecTrack = 0; iRecTrack <fNRecTracks; iRecTrack++) {
948 track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
949 trackParamVertex = *((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First()));
950 vertex = track->GetVertex();
951 if (vertex) { // if track as a vertex defined, use it
952 vX = vertex->GetNonBendingCoor();
953 vY = vertex->GetBendingCoor();
955 } else { // else vertex = (0.,0.,0.)
960 AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, vX, vY, vZ);
961 track->SetTrackParamAtVertex(&trackParamVertex);
963 if (AliLog::GetGlobalDebugLevel() > 0) {
964 cout << "FollowTracks: track candidate(0..): " << iRecTrack
965 << " after extrapolation to vertex" << endl;
966 track->RecursiveDump();
972 //__________________________________________________________________________
973 void AliMUONTrackReconstructor::FillMUONTrack()
975 /// Fill fHitForRecAtHit of AliMUONTrack's
977 AliMUONTrackParam *trackParamAtHit;
978 track = (AliMUONTrack*) fRecTracksPtr->First();
980 trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
981 while (trackParamAtHit) {
982 track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr());
983 trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit));
985 track = (AliMUONTrack*) fRecTracksPtr->After(track);
990 //__________________________________________________________________________
991 void AliMUONTrackReconstructor::EventDump(void)
993 /// Dump reconstructed event (track parameters at vertex and at first hit),
994 /// and the particle parameters
996 AliMUONTrackParam *trackParam, *trackParam1;
997 Double_t bendingSlope, nonBendingSlope, pYZ;
998 Double_t pX, pY, pZ, x, y, z, c;
999 Int_t trackIndex, nTrackHits;
1001 AliDebug(1,"****** enter EventDump ******");
1002 AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks));
1004 fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole
1005 // Loop over reconstructed tracks
1006 for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) {
1007 AliDebug(1, Form("track number: %d", trackIndex));
1008 // function for each track for modularity ????
1009 track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]);
1010 nTrackHits = track->GetNTrackHits();
1011 AliDebug(1, Form("Number of track hits: %d ", nTrackHits));
1012 // track parameters at Vertex
1013 trackParam = track->GetTrackParamAtVertex();
1014 x = trackParam->GetNonBendingCoor();
1015 y = trackParam->GetBendingCoor();
1016 z = trackParam->GetZ();
1017 bendingSlope = trackParam->GetBendingSlope();
1018 nonBendingSlope = trackParam->GetNonBendingSlope();
1019 pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
1020 pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
1021 pX = pZ * nonBendingSlope;
1022 pY = pZ * bendingSlope;
1023 c = TMath::Sign(1.0, trackParam->GetInverseBendingMomentum());
1024 AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1025 z, x, y, pX, pY, pZ, c));
1027 // track parameters at first hit
1028 trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First();
1029 x = trackParam1->GetNonBendingCoor();
1030 y = trackParam1->GetBendingCoor();
1031 z = trackParam1->GetZ();
1032 bendingSlope = trackParam1->GetBendingSlope();
1033 nonBendingSlope = trackParam1->GetNonBendingSlope();
1034 pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum());
1035 pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
1036 pX = pZ * nonBendingSlope;
1037 pY = pZ * bendingSlope;
1038 c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum());
1039 AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n",
1040 z, x, y, pX, pY, pZ, c));
1042 // informations about generated particles NO !!!!!!!!
1044 // for (Int_t iPart = 0; iPart < np; iPart++) {
1045 // p = gAlice->Particle(iPart);
1046 // printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n",
1047 // iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode());