]>
Commit | Line | Data |
---|---|---|
a9e2aefa | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
cc87ebcd | 16 | /* $Id$ */ |
17 | ||
7945aae7 | 18 | /// \class AliMUONTrackReconstructor |
19 | /// MUON track reconstructor using the original method | |
20 | /// | |
21 | /// This class contains as data: | |
22 | /// - the parameters for the track reconstruction | |
23 | /// | |
24 | /// It contains as methods, among others: | |
25 | /// - MakeTracks to build the tracks | |
26 | /// | |
a9e2aefa | 27 | |
29f1b13a | 28 | #include "AliMUONTrackReconstructor.h" |
7ec3b9cf | 29 | |
29fc2c86 | 30 | #include "AliMUONConstants.h" |
a9e2aefa | 31 | #include "AliMUONRawCluster.h" |
2457f726 | 32 | #include "AliMUONHitForRec.h" |
208f139e | 33 | #include "AliMUONObjectPair.h" |
a9e2aefa | 34 | #include "AliMUONTrack.h" |
37827b29 | 35 | #include "AliMUONTrackParam.h" |
36 | #include "AliMUONTrackExtrap.h" | |
8c343c7c | 37 | #include "AliLog.h" |
af743f54 | 38 | |
208f139e | 39 | #include <TMinuit.h> |
af743f54 | 40 | #include <Riostream.h> |
41 | #include <TMatrixD.h> | |
a9e2aefa | 42 | |
de2cd600 | 43 | // Functions to be minimized with Minuit |
44 | void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); | |
45 | void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); | |
46 | ||
de2cd600 | 47 | Double_t MultipleScatteringAngle2(AliMUONTrackParam *param); |
48 | ||
7945aae7 | 49 | /// \cond CLASSIMP |
29f1b13a | 50 | ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context |
7945aae7 | 51 | /// \endcond |
a9e2aefa | 52 | |
2457f726 | 53 | //************* Defaults parameters for reconstruction |
208f139e | 54 | const Double_t AliMUONTrackReconstructor::fgkMaxNormChi2 = 100.0; |
55 | const Bool_t AliMUONTrackReconstructor::fgkTrackAllTracks = kFALSE; | |
2457f726 | 56 | |
3bc8b580 | 57 | //__________________________________________________________________________ |
7ec3b9cf | 58 | AliMUONTrackReconstructor::AliMUONTrackReconstructor() |
59 | : AliMUONVTrackReconstructor() | |
a9e2aefa | 60 | { |
2457f726 | 61 | /// Constructor for class AliMUONTrackReconstructor |
de2cd600 | 62 | |
a9e2aefa | 63 | // Memory allocation for the TClonesArray of reconstructed tracks |
a9e2aefa | 64 | fRecTracksPtr = new TClonesArray("AliMUONTrack", 10); |
a9e2aefa | 65 | } |
9cbdf048 | 66 | |
7ec3b9cf | 67 | //__________________________________________________________________________ |
29f1b13a | 68 | AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void) |
a9e2aefa | 69 | { |
2457f726 | 70 | /// Destructor for class AliMUONTrackReconstructor |
de2cd600 | 71 | delete fRecTracksPtr; |
a9e2aefa | 72 | } |
a9e2aefa | 73 | |
7ec3b9cf | 74 | //__________________________________________________________________________ |
29f1b13a | 75 | void AliMUONTrackReconstructor::MakeTracks(void) |
a9e2aefa | 76 | { |
2457f726 | 77 | /// To make the tracks from the list of segments and points in all stations |
b840996b | 78 | AliDebug(1,"Enter MakeTracks"); |
2457f726 | 79 | // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5 |
80 | MakeTrackCandidates(); | |
81 | // Follow tracks in stations(1..) 3, 2 and 1 | |
82 | FollowTracks(); | |
83 | // Remove double tracks | |
84 | RemoveDoubleTracks(); | |
208f139e | 85 | // Fill out the AliMUONTrack's |
86 | FillMUONTrack(); | |
276c44b7 | 87 | } |
88 | ||
de2cd600 | 89 | //__________________________________________________________________________ |
90 | void AliMUONTrackReconstructor::MakeTrackCandidates(void) | |
91 | { | |
208f139e | 92 | /// To make track candidates: |
93 | /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4. | |
94 | /// Good candidates are made of at least three hitForRec's. | |
95 | /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks. | |
96 | TClonesArray *segments; | |
97 | AliMUONObjectPair *segment; | |
98 | AliMUONHitForRec *hitForRec1, *hitForRec2; | |
99 | AliMUONTrack *track; | |
100 | AliMUONTrackParam *trackParamAtFirstHit; | |
101 | Int_t iCandidate = 0; | |
102 | ||
de2cd600 | 103 | AliDebug(1,"Enter MakeTrackCandidates"); |
de2cd600 | 104 | |
208f139e | 105 | // Loop over stations(1..) 5 and 4 and make track candidates |
7ec3b9cf | 106 | for (Int_t istat=4; istat>=3; istat--) |
107 | { | |
208f139e | 108 | // Make segments in the station |
109 | segments = MakeSegmentsInStation(istat); | |
110 | // Loop over segments | |
7ec3b9cf | 111 | for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) |
112 | { | |
208f139e | 113 | AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate)); |
114 | // Transform segments to tracks and put them at the end of fRecTracksPtr | |
115 | segment = (AliMUONObjectPair*) ((*segments)[iseg]); | |
116 | hitForRec1 = (AliMUONHitForRec*) segment->First(); | |
117 | hitForRec2 = (AliMUONHitForRec*) segment->Second(); | |
118 | track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(hitForRec1, hitForRec2); | |
a9e2aefa | 119 | fNRecTracks++; |
208f139e | 120 | // Add MCS effects in parameter covariances |
121 | trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); | |
8cde4af5 | 122 | AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.); |
208f139e | 123 | // Printout for debuging |
7ec3b9cf | 124 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) |
125 | { | |
208f139e | 126 | cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl; |
127 | trackParamAtFirstHit->GetCovariances()->Print(); | |
128 | } | |
129 | // Look for compatible hitForRec(s) in the other station | |
130 | FollowTrackInStation(track,7-istat); | |
a9e2aefa | 131 | } |
208f139e | 132 | // delete the array of segments |
133 | delete segments; | |
134 | } | |
135 | fRecTracksPtr->Compress(); // this is essential before checking tracks | |
136 | ||
137 | // Keep all different tracks or only the best ones as required | |
138 | if (fgkTrackAllTracks) RemoveIdenticalTracks(); | |
139 | else RemoveDoubleTracks(); | |
140 | ||
141 | AliDebug(1,Form("Number of good candidates = %d",fNRecTracks)); | |
142 | ||
a9e2aefa | 143 | } |
144 | ||
145 | //__________________________________________________________________________ | |
208f139e | 146 | void AliMUONTrackReconstructor::RemoveIdenticalTracks(void) |
a9e2aefa | 147 | { |
208f139e | 148 | /// To remove identical tracks: |
149 | /// Tracks are considered identical if they have all their hits in common. | |
150 | /// One keeps the track with the larger number of hits if need be | |
151 | AliMUONTrack *track1, *track2, *trackToRemove; | |
152 | Int_t hitsInCommon, nHits1, nHits2; | |
153 | Bool_t removedTrack1; | |
154 | // Loop over first track of the pair | |
155 | track1 = (AliMUONTrack*) fRecTracksPtr->First(); | |
156 | while (track1) { | |
157 | removedTrack1 = kFALSE; | |
158 | nHits1 = track1->GetNTrackHits(); | |
159 | // Loop over second track of the pair | |
160 | track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); | |
161 | while (track2) { | |
162 | nHits2 = track2->GetNTrackHits(); | |
163 | // number of hits in common between two tracks | |
164 | hitsInCommon = track1->HitsInCommon(track2); | |
165 | // check for identical tracks | |
166 | if ((hitsInCommon == nHits1) || (hitsInCommon == nHits2)) { | |
167 | // decide which track to remove | |
168 | if (nHits2 > nHits1) { | |
169 | // remove track1 and continue the first loop with the track next to track1 | |
170 | trackToRemove = track1; | |
171 | track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); | |
172 | fRecTracksPtr->Remove(trackToRemove); | |
173 | fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards | |
174 | fNRecTracks--; | |
175 | removedTrack1 = kTRUE; | |
176 | break; | |
177 | } else { | |
178 | // remove track2 and continue the second loop with the track next to track2 | |
179 | trackToRemove = track2; | |
180 | track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); | |
181 | fRecTracksPtr->Remove(trackToRemove); | |
182 | fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards | |
183 | fNRecTracks--; | |
184 | } | |
185 | } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); | |
186 | } // track2 | |
187 | if (removedTrack1) continue; | |
188 | track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); | |
189 | } // track1 | |
190 | return; | |
a9e2aefa | 191 | } |
192 | ||
193 | //__________________________________________________________________________ | |
208f139e | 194 | void AliMUONTrackReconstructor::RemoveDoubleTracks(void) |
a9e2aefa | 195 | { |
208f139e | 196 | /// To remove double tracks: |
197 | /// Tracks are considered identical if more than half of the hits of the track | |
198 | /// which has the smaller number of hits are in common with the other track. | |
199 | /// Among two identical tracks, one keeps the track with the larger number of hits | |
200 | /// or, if these numbers are equal, the track with the minimum chi2. | |
201 | AliMUONTrack *track1, *track2, *trackToRemove; | |
202 | Int_t hitsInCommon, nHits1, nHits2; | |
203 | Bool_t removedTrack1; | |
204 | // Loop over first track of the pair | |
205 | track1 = (AliMUONTrack*) fRecTracksPtr->First(); | |
206 | while (track1) { | |
207 | removedTrack1 = kFALSE; | |
208 | nHits1 = track1->GetNTrackHits(); | |
209 | // Loop over second track of the pair | |
210 | track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); | |
211 | while (track2) { | |
212 | nHits2 = track2->GetNTrackHits(); | |
213 | // number of hits in common between two tracks | |
214 | hitsInCommon = track1->HitsInCommon(track2); | |
215 | // check for identical tracks | |
216 | if (((nHits1 < nHits2) && (2 * hitsInCommon > nHits1)) || (2 * hitsInCommon > nHits2)) { | |
217 | // decide which track to remove | |
218 | if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() <= track2->GetFitFMin()))) { | |
219 | // remove track2 and continue the second loop with the track next to track2 | |
220 | trackToRemove = track2; | |
221 | track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); | |
222 | fRecTracksPtr->Remove(trackToRemove); | |
223 | fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards | |
224 | fNRecTracks--; | |
225 | } else { | |
226 | // else remove track1 and continue the first loop with the track next to track1 | |
227 | trackToRemove = track1; | |
228 | track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); | |
229 | fRecTracksPtr->Remove(trackToRemove); | |
230 | fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards | |
231 | fNRecTracks--; | |
232 | removedTrack1 = kTRUE; | |
233 | break; | |
234 | } | |
235 | } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); | |
236 | } // track2 | |
237 | if (removedTrack1) continue; | |
238 | track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); | |
239 | } // track1 | |
240 | return; | |
a9e2aefa | 241 | } |
242 | ||
243 | //__________________________________________________________________________ | |
29f1b13a | 244 | void AliMUONTrackReconstructor::FollowTracks(void) |
a9e2aefa | 245 | { |
2457f726 | 246 | /// Follow tracks in stations(1..) 3, 2 and 1 |
208f139e | 247 | AliDebug(1,"Enter FollowTracks"); |
248 | ||
04b5ea16 | 249 | AliMUONTrack *track, *nextTrack; |
208f139e | 250 | AliMUONTrackParam *trackParamAtFirstHit; |
251 | Double_t numberOfDegFree, chi2Norm; | |
af743f54 | 252 | Int_t currentNRecTracks; |
208f139e | 253 | |
254 | for (Int_t station = 2; station >= 0; station--) { | |
255 | // Save the actual number of reconstructed track in case of | |
256 | // tracks are added or suppressed during the tracking procedure | |
257 | // !! Do not compress fRecTracksPtr until the end of the loop over tracks !! | |
af743f54 | 258 | currentNRecTracks = fNRecTracks; |
259 | for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) { | |
208f139e | 260 | AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1)); |
261 | track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack); | |
262 | // Fit the track: | |
263 | // Do not take into account the multiple scattering to speed up the fit | |
264 | // Calculate the track parameter covariance matrix | |
265 | // If "station" is station(1..) 3 then use the vertex to better constrain the fit | |
266 | if (station==2) { | |
267 | SetVertexForFit(track); | |
268 | track->SetFitWithVertex(kTRUE); | |
269 | } else track->SetFitWithVertex(kFALSE); | |
270 | Fit(track,kFALSE, kTRUE); | |
271 | // Remove the track if the normalized chi2 is too high | |
272 | numberOfDegFree = (2. * track->GetNTrackHits() - 5.); | |
273 | if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree; | |
274 | else chi2Norm = 1.e10; | |
275 | if (chi2Norm > fgkMaxNormChi2) { | |
276 | fRecTracksPtr->Remove(track); | |
277 | fNRecTracks--; | |
278 | continue; | |
279 | } | |
280 | // Add MCS effects in parameter covariances | |
281 | trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); | |
8cde4af5 | 282 | AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.); |
208f139e | 283 | // Printout for debuging |
284 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
285 | cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl; | |
286 | trackParamAtFirstHit->GetCovariances()->Print(); | |
287 | } | |
288 | // Look for compatible hitForRec in station(0..) "station" | |
289 | FollowTrackInStation(track,station); | |
290 | } | |
291 | // Compress fRecTracksPtr for the next step | |
292 | fRecTracksPtr->Compress(); | |
293 | // Keep only the best tracks if required | |
294 | if (!fgkTrackAllTracks) RemoveDoubleTracks(); | |
295 | } | |
296 | ||
297 | // Last fit of track candidates with all station | |
298 | // Take into account the multiple scattering and remove bad tracks | |
299 | Int_t trackIndex = -1; | |
04b5ea16 | 300 | track = (AliMUONTrack*) fRecTracksPtr->First(); |
04b5ea16 | 301 | while (track) { |
04b5ea16 | 302 | trackIndex++; |
303 | nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track | |
208f139e | 304 | track->SetFitWithVertex(kFALSE); // just to be sure |
b58638a9 | 305 | Fit(track,kTRUE, kTRUE); |
208f139e | 306 | // Printout for debuging |
307 | if (AliLog::GetGlobalDebugLevel() >= 3) { | |
308 | cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl; | |
a9e2aefa | 309 | track->RecursiveDump(); |
208f139e | 310 | } |
311 | // Remove the track if the normalized chi2 is too high | |
312 | numberOfDegFree = (2.0 * track->GetNTrackHits() - 5); | |
313 | if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree; | |
314 | else chi2Norm = 1.e10; | |
315 | if (chi2Norm > fgkMaxNormChi2) { | |
316 | fRecTracksPtr->Remove(track); | |
317 | fNRecTracks--; | |
a9e2aefa | 318 | } |
208f139e | 319 | track = nextTrack; |
320 | } | |
321 | fRecTracksPtr->Compress(); | |
322 | ||
323 | } | |
1a38e749 | 324 | |
208f139e | 325 | //__________________________________________________________________________ |
326 | void AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation) | |
327 | { | |
328 | /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s) | |
329 | /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks: | |
330 | /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of | |
331 | /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure. | |
332 | /// Remove the obsolete "trackCandidate" at the end. | |
333 | /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority. | |
334 | AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1)); | |
335 | ||
336 | Int_t ch1 = 2*nextStation; | |
337 | Int_t ch2 = 2*nextStation+1; | |
338 | Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2); | |
339 | Double_t chi2WithOneHitForRec = 1.e10; | |
340 | Double_t chi2WithTwoHitForRec = 1.e10; | |
341 | Double_t maxChi2WithOneHitForRec = 2.*fgkMaxNormChi2; // 2 because 2 quantities in chi2 | |
342 | Double_t maxChi2WithTwoHitForRec = 4.*fgkMaxNormChi2; // 4 because 4 quantities in chi2 | |
343 | Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec; | |
344 | Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec; | |
345 | AliMUONTrack *newTrack = 0x0; | |
346 | AliMUONHitForRec *hitForRecCh1, *hitForRecCh2; | |
347 | AliMUONHitForRec *bestHitForRec1 = 0x0, *bestHitForRec2 = 0x0; | |
847cbaef | 348 | Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]]; |
349 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE; | |
208f139e | 350 | // |
351 | //Extrapolate trackCandidate to chamber "ch2" to save computing time in the next steps | |
352 | AliMUONTrackParam *extrapTrackParamPtr = trackCandidate->GetExtrapTrackParam(); | |
353 | *extrapTrackParamPtr = *((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First())); | |
354 | AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, zCh2); | |
355 | AliMUONTrackParam extrapTrackParamSave(*extrapTrackParamPtr); | |
356 | // | |
357 | // Printout for debuging | |
358 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
359 | TMatrixD* paramCovForDebug = extrapTrackParamPtr->GetCovariances(); | |
360 | cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl; | |
361 | paramCovForDebug->Print(); | |
362 | } | |
363 | // | |
364 | // look for candidates in chamber 2 | |
365 | // Printout for debuging | |
366 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
367 | cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl; | |
368 | } | |
369 | for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) { | |
370 | hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2); | |
371 | // extrapolate track parameters and covariances only once for this hit | |
372 | AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, hitForRecCh2->GetZ()); | |
373 | chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh2); | |
374 | // if good chi2 then try to attach a hitForRec in the other chamber too | |
375 | if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) { | |
376 | // Printout for debuging | |
377 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
378 | cout << "FollowTrackInStation: look for second hits in chamber(1..): " << ch1+1 << endl; | |
a9e2aefa | 379 | } |
208f139e | 380 | Bool_t foundSecondHit = kFALSE; |
381 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { | |
382 | hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); | |
383 | chi2WithTwoHitForRec = trackCandidate->TryTwoHitForRec(hitForRecCh2, hitForRecCh1); // order hits like that to save computing time | |
384 | // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate" | |
385 | if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) { | |
386 | foundSecondHit = kTRUE; | |
387 | if (fgkTrackAllTracks) { | |
388 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
389 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate); | |
390 | fNRecTracks++; | |
af743f54 | 391 | AliMUONTrackParam trackParam1(extrapTrackParamSave); |
392 | AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh1->GetZ()); | |
393 | newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh1); | |
394 | AliMUONTrackParam trackParam2(extrapTrackParamSave); | |
395 | AliMUONTrackExtrap::ExtrapToZ(&trackParam2, hitForRecCh2->GetZ()); | |
396 | newTrack->AddTrackParamAtHit(&trackParam2,hitForRecCh2); | |
208f139e | 397 | // Sort TrackParamAtHit according to increasing Z |
398 | newTrack->GetTrackParamAtHit()->Sort(); | |
399 | // Update the chi2 of the new track | |
400 | if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithTwoHitForRec); | |
401 | else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithTwoHitForRec); | |
847cbaef | 402 | // Tag hitForRecCh1 as used |
403 | hitForRecCh1Used[hit1] = kTRUE; | |
208f139e | 404 | // Printout for debuging |
405 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
406 | cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1 | |
407 | << " (Chi2 = " << chi2WithTwoHitForRec << ")" << endl; | |
408 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
a9e2aefa | 409 | } |
208f139e | 410 | } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) { |
411 | // keep track of the best couple of hits | |
412 | bestChi2WithTwoHitForRec = chi2WithTwoHitForRec; | |
413 | bestHitForRec1 = hitForRecCh1; | |
414 | bestHitForRec2 = hitForRecCh2; | |
415 | } | |
a9e2aefa | 416 | } |
d0bfce8d | 417 | } |
208f139e | 418 | // if no hitForRecCh1 found then consider to add hitForRecCh2 only |
419 | if (!foundSecondHit) { | |
420 | if (fgkTrackAllTracks) { | |
421 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
422 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate); | |
423 | fNRecTracks++; | |
af743f54 | 424 | AliMUONTrackParam trackParam1(extrapTrackParamSave); |
425 | AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh2->GetZ()); | |
426 | newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh2); | |
208f139e | 427 | // Sort TrackParamAtHit according to increasing Z |
428 | newTrack->GetTrackParamAtHit()->Sort(); | |
429 | // Update the chi2 of the new track | |
430 | if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec); | |
431 | else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec); | |
432 | // Printout for debuging | |
433 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
434 | cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1 | |
435 | << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl; | |
436 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
437 | } | |
438 | } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) { | |
439 | // keep track of the best single hitForRec except if a couple | |
440 | // of hits has already been found (i.e. bestHitForRec2!=0x0) | |
441 | bestChi2WithOneHitForRec = chi2WithOneHitForRec; | |
442 | bestHitForRec1 = hitForRecCh2; | |
443 | } | |
d0bfce8d | 444 | } |
208f139e | 445 | } |
446 | // reset the extrapolated track parameter for next step | |
447 | trackCandidate->SetExtrapTrackParam(&extrapTrackParamSave); | |
448 | } | |
449 | // | |
450 | // look for candidates in chamber 1 not already attached to a track | |
451 | // if we want to keep all possible tracks or if no good couple of hitForRec has been found | |
452 | if (fgkTrackAllTracks || !bestHitForRec2) { | |
453 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
454 | cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl; | |
455 | } | |
456 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { | |
457 | hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); | |
847cbaef | 458 | if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used |
208f139e | 459 | chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh1); |
460 | // if good chi2 then create a new track by adding the good hitForRec in "ch1" to the "trackCandidate" | |
461 | // We do not try to attach a hitForRec in the other chamber too since it has already been done above | |
462 | if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) { | |
463 | if (fgkTrackAllTracks) { | |
464 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
465 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate); | |
466 | fNRecTracks++; | |
af743f54 | 467 | AliMUONTrackParam trackParam1(extrapTrackParamSave); |
468 | AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh1->GetZ()); | |
469 | newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh1); | |
208f139e | 470 | // Sort TrackParamAtHit according to increasing Z |
471 | newTrack->GetTrackParamAtHit()->Sort(); | |
472 | // Update the chi2 of the new track | |
473 | if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec); | |
474 | else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec); | |
475 | // Printout for debuging | |
476 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
477 | cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1 | |
478 | << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl; | |
479 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
480 | } | |
481 | } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) { | |
482 | // keep track of the best single hitForRec except if a couple | |
483 | // of hits has already been found (i.e. bestHitForRec1!=0x0) | |
484 | bestChi2WithOneHitForRec = chi2WithOneHitForRec; | |
485 | bestHitForRec1 = hitForRecCh1; | |
486 | } | |
d0bfce8d | 487 | } |
208f139e | 488 | } |
489 | } | |
490 | // | |
491 | // fill out the best track if required else clean up the fRecTracksPtr array | |
492 | if (!fgkTrackAllTracks && bestHitForRec1) { | |
af743f54 | 493 | AliMUONTrackParam trackParam1(extrapTrackParamSave); |
494 | AliMUONTrackExtrap::ExtrapToZ(&trackParam1, bestHitForRec1->GetZ()); | |
495 | trackCandidate->AddTrackParamAtHit(&trackParam1,bestHitForRec1); | |
208f139e | 496 | if (bestHitForRec2) { |
af743f54 | 497 | AliMUONTrackParam trackParam2(extrapTrackParamSave); |
498 | AliMUONTrackExtrap::ExtrapToZ(&trackParam2, bestHitForRec2->GetZ()); | |
499 | trackCandidate->AddTrackParamAtHit(&trackParam2,bestHitForRec2); | |
208f139e | 500 | // Sort TrackParamAtHit according to increasing Z |
501 | trackCandidate->GetTrackParamAtHit()->Sort(); | |
502 | // Update the chi2 of the new track | |
503 | if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithTwoHitForRec); | |
504 | else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithTwoHitForRec); | |
505 | // Printout for debuging | |
506 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
507 | cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1 | |
508 | << " (Chi2 = " << bestChi2WithTwoHitForRec << ")" << endl; | |
509 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
a9e2aefa | 510 | } |
208f139e | 511 | } else { |
512 | // Sort TrackParamAtHit according to increasing Z | |
513 | trackCandidate->GetTrackParamAtHit()->Sort(); | |
514 | // Update the chi2 of the new track | |
515 | if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithOneHitForRec); | |
516 | else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithOneHitForRec); | |
517 | // Printout for debuging | |
518 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
519 | cout << "FollowTrackInStation: added the best hit in station(1..): " << nextStation+1 | |
520 | << " (Chi2 = " << bestChi2WithOneHitForRec << ")" << endl; | |
521 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); | |
04b5ea16 | 522 | } |
208f139e | 523 | } |
524 | } else { | |
525 | fRecTracksPtr->Remove(trackCandidate); // obsolete track | |
526 | fNRecTracks--; | |
527 | } | |
528 | ||
04b5ea16 | 529 | } |
530 | ||
de2cd600 | 531 | //__________________________________________________________________________ |
208f139e | 532 | void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate) |
de2cd600 | 533 | { |
208f139e | 534 | /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate" |
535 | /// Compute the vertex resolution from natural vertex dispersion and | |
536 | /// multiple scattering effets according to trackCandidate path in absorber | |
537 | /// It is necessary to account for multiple scattering effects here instead of during the fit of | |
538 | /// the "trackCandidate" to do not influence the result by changing track resolution at vertex | |
539 | AliDebug(1,"Enter SetVertexForFit"); | |
de2cd600 | 540 | |
208f139e | 541 | Double_t nonBendingReso2 = fNonBendingVertexDispersion * fNonBendingVertexDispersion; |
542 | Double_t bendingReso2 = fBendingVertexDispersion * fBendingVertexDispersion; | |
543 | // add multiple scattering effets | |
544 | AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First()))); | |
545 | paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering | |
546 | AliMUONTrackExtrap::ExtrapToVertexUncorrected(¶mAtVertex,0.); | |
547 | TMatrixD* paramCov = paramAtVertex.GetCovariances(); | |
548 | nonBendingReso2 += (*paramCov)(0,0); | |
549 | bendingReso2 += (*paramCov)(2,2); | |
550 | // Set the vertex | |
551 | AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default | |
552 | vertex.SetNonBendingReso2(nonBendingReso2); | |
553 | vertex.SetBendingReso2(bendingReso2); | |
554 | trackCandidate->SetVertex(&vertex); | |
555 | } | |
556 | ||
557 | //__________________________________________________________________________ | |
558 | void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov) | |
559 | { | |
560 | /// Fit the track "track" with or without multiple Coulomb scattering according to "includeMCS". | |
de2cd600 | 561 | |
208f139e | 562 | Double_t benC, errorParam, invBenP, nonBenC, x, y; |
de2cd600 | 563 | AliMUONTrackParam *trackParam; |
208f139e | 564 | Double_t arg[1], fedm, errdef, fitFMin; |
565 | Int_t npari, nparx; | |
566 | Int_t status, covStatus; | |
567 | ||
8cde4af5 | 568 | // Instantiate gMinuit if not already done |
569 | if (!gMinuit) gMinuit = new TMinuit(6); | |
208f139e | 570 | // Clear MINUIT parameters |
571 | gMinuit->mncler(); | |
572 | // Give the fitted track to MINUIT | |
573 | gMinuit->SetObjectFit(track); | |
574 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
575 | // Define print level | |
576 | arg[0] = 1; | |
577 | gMinuit->mnexcm("SET PRI", arg, 1, status); | |
578 | // Print covariance matrix | |
579 | gMinuit->mnexcm("SHO COV", arg, 0, status); | |
580 | } else { | |
581 | arg[0] = -1; | |
582 | gMinuit->mnexcm("SET PRI", arg, 1, status); | |
583 | } | |
de2cd600 | 584 | // No warnings |
208f139e | 585 | gMinuit->mnexcm("SET NOW", arg, 0, status); |
586 | // Define strategy | |
587 | //arg[0] = 2; | |
588 | //gMinuit->mnexcm("SET STR", arg, 1, status); | |
589 | ||
590 | // Switch between available FCN according to "includeMCS" | |
591 | if (includeMCS) gMinuit->SetFCN(TrackChi2MCS); | |
592 | else gMinuit->SetFCN(TrackChi2); | |
593 | ||
594 | // Set fitted parameters (!! The order is very important for the covariance matrix !!) | |
595 | trackParam = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); | |
de2cd600 | 596 | // could be tried with no limits for the search (min=max=0) ???? |
208f139e | 597 | // mandatory limits in non Bending to avoid NaN values of parameters |
598 | gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status); | |
599 | gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status); | |
600 | // mandatory limits in Bending to avoid NaN values of parameters | |
601 | gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status); | |
602 | gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status); | |
603 | gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status); | |
604 | ||
de2cd600 | 605 | // minimization |
208f139e | 606 | gMinuit->mnexcm("MIGRAD", arg, 0, status); |
607 | ||
608 | // Calculate the covariance matrix more accurately if required | |
609 | if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status); | |
8cde4af5 | 610 | |
de2cd600 | 611 | // get results into "invBenP", "benC", "nonBenC" ("x", "y") |
208f139e | 612 | gMinuit->GetParameter(0, x, errorParam); |
613 | trackParam->SetNonBendingCoor(x); | |
614 | gMinuit->GetParameter(1, nonBenC, errorParam); | |
de2cd600 | 615 | trackParam->SetNonBendingSlope(nonBenC); |
208f139e | 616 | gMinuit->GetParameter(2, y, errorParam); |
617 | trackParam->SetBendingCoor(y); | |
618 | gMinuit->GetParameter(3, benC, errorParam); | |
619 | trackParam->SetBendingSlope(benC); | |
620 | gMinuit->GetParameter(4, invBenP, errorParam); | |
621 | trackParam->SetInverseBendingMomentum(invBenP); | |
622 | ||
de2cd600 | 623 | // global result of the fit |
208f139e | 624 | gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus); |
625 | track->SetFitFMin(fitFMin); | |
626 | ||
627 | // Get the covariance matrix if required | |
628 | if (calcCov) { | |
629 | // Covariance matrix according to HESSE status | |
630 | // If problem then keep only the diagonal terms (variances) | |
631 | Double_t matrix[5][5]; | |
632 | gMinuit->mnemat(&matrix[0][0],5); | |
633 | if (covStatus == 3) trackParam->SetCovariances(matrix); | |
634 | else trackParam->SetVariances(matrix); | |
635 | } else *(trackParam->GetCovariances()) = 0; | |
636 | ||
de2cd600 | 637 | } |
638 | ||
639 | //__________________________________________________________________________ | |
208f139e | 640 | void TrackChi2(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) |
de2cd600 | 641 | { |
2457f726 | 642 | /// Return the "Chi2" to be minimized with Minuit for track fitting, |
643 | /// with "NParam" parameters | |
644 | /// and their current values in array pointed to by "Param". | |
645 | /// Assumes that the track hits are sorted according to increasing Z. | |
646 | /// Track parameters at each TrackHit are updated accordingly. | |
647 | /// Multiple Coulomb scattering is not taken into account | |
208f139e | 648 | |
649 | AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit(); | |
650 | // AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit(); | |
de2cd600 | 651 | AliMUONTrackParam param1; |
2457f726 | 652 | AliMUONTrackParam* trackParamAtHit; |
653 | AliMUONHitForRec* hitForRec; | |
de2cd600 | 654 | Chi2 = 0.0; // initialize Chi2 |
208f139e | 655 | Double_t dX, dY; |
656 | ||
de2cd600 | 657 | // copy of track parameters to be fitted |
208f139e | 658 | param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First())); |
659 | param1.SetNonBendingCoor(Param[0]); | |
660 | param1.SetNonBendingSlope(Param[1]); | |
661 | param1.SetBendingCoor(Param[2]); | |
662 | param1.SetBendingSlope(Param[3]); | |
663 | param1.SetInverseBendingMomentum(Param[4]); | |
664 | ||
665 | // Take the vertex into account in the fit if required | |
666 | if (trackBeingFitted->GetFitWithVertex()) { | |
667 | AliMUONTrackParam paramAtVertex(param1); | |
668 | AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.); | |
669 | AliMUONHitForRec *vertex = trackBeingFitted->GetVertex(); | |
670 | if (!vertex) { | |
671 | cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl; | |
672 | exit(-1); | |
673 | } | |
674 | dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor(); | |
675 | dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor(); | |
676 | Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2(); | |
de2cd600 | 677 | } |
208f139e | 678 | |
de2cd600 | 679 | // Follow track through all planes of track hits |
2457f726 | 680 | trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()); |
681 | while (trackParamAtHit) { | |
682 | hitForRec = trackParamAtHit->GetHitForRecPtr(); | |
683 | // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit | |
37827b29 | 684 | AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ()); |
de2cd600 | 685 | // update track parameters of the current hit |
2457f726 | 686 | trackParamAtHit->SetTrackParam(param1); |
de2cd600 | 687 | // Increment Chi2 |
688 | // done hit per hit, with hit resolution, | |
689 | // and not with point and angle like in "reco_muon.F" !!!! | |
208f139e | 690 | dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor(); |
691 | dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor(); | |
2457f726 | 692 | Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2(); |
693 | trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit)); | |
de2cd600 | 694 | } |
695 | } | |
696 | ||
697 | //__________________________________________________________________________ | |
208f139e | 698 | void TrackChi2MCS(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) |
de2cd600 | 699 | { |
2457f726 | 700 | /// Return the "Chi2" to be minimized with Minuit for track fitting, |
701 | /// with "NParam" parameters | |
702 | /// and their current values in array pointed to by "Param". | |
703 | /// Assumes that the track hits are sorted according to increasing Z. | |
704 | /// Track parameters at each TrackHit are updated accordingly. | |
705 | /// Multiple Coulomb scattering is taken into account with covariance matrix. | |
208f139e | 706 | |
707 | AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit(); | |
708 | // AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit(); | |
de2cd600 | 709 | AliMUONTrackParam param1; |
2457f726 | 710 | AliMUONTrackParam* trackParamAtHit; |
711 | AliMUONHitForRec* hitForRec; | |
de2cd600 | 712 | Chi2 = 0.0; // initialize Chi2 |
de2cd600 | 713 | Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3; |
714 | Double_t z1, z2, z3; | |
2457f726 | 715 | AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3; |
716 | AliMUONHitForRec *hitForRec1, *hitForRec2; | |
de2cd600 | 717 | Double_t hbc1, hbc2, pbc1, pbc2; |
718 | Double_t hnbc1, hnbc2, pnbc1, pnbc2; | |
719 | Int_t numberOfHit = trackBeingFitted->GetNTrackHits(); | |
720 | TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit); | |
721 | TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit); | |
722 | Double_t *msa2 = new Double_t[numberOfHit]; | |
208f139e | 723 | |
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]); | |
731 | ||
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(); | |
737 | if (!vertex) { | |
738 | cout<<"Error in TrackChi2MCS: Want to use the vertex in tracking but it has not been created!!"<<endl; | |
739 | exit(-1); | |
740 | } | |
741 | Double_t dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor(); | |
742 | Double_t dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor(); | |
743 | Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2(); | |
744 | } | |
745 | ||
746 | // Predicted coordinates and multiple scattering angles are first calculated | |
de2cd600 | 747 | for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) { |
2457f726 | 748 | trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber)); |
749 | hitForRec = trackParamAtHit->GetHitForRecPtr(); | |
750 | // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit | |
37827b29 | 751 | AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ()); |
de2cd600 | 752 | // update track parameters of the current hit |
2457f726 | 753 | trackParamAtHit->SetTrackParam(param1); |
de2cd600 | 754 | // square of multiple scattering angle at current hit, with one chamber |
755 | msa2[hitNumber] = MultipleScatteringAngle2(¶m1); | |
756 | // correction for eventual missing hits or multiple hits in a chamber, | |
757 | // according to the number of chambers | |
758 | // between the current hit and the previous one | |
2457f726 | 759 | chCurrent = hitForRec->GetChamberNumber(); |
de2cd600 | 760 | if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev); |
761 | chPrev = chCurrent; | |
762 | } | |
763 | ||
764 | // Calculates the covariance matrix | |
765 | for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { | |
2457f726 | 766 | trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1)); |
767 | hitForRec1 = trackParamAtHit1->GetHitForRecPtr(); | |
768 | z1 = hitForRec1->GetZ(); | |
de2cd600 | 769 | for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) { |
2457f726 | 770 | trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2)); |
771 | z2 = trackParamAtHit2->GetHitForRecPtr()->GetZ(); | |
de2cd600 | 772 | // initialization to 0 (diagonal plus upper triangular part) |
773 | (*covBending)(hitNumber2, hitNumber1) = 0.0; | |
774 | // contribution from multiple scattering in bending plane: | |
775 | // loop over upstream hits | |
776 | for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) { | |
2457f726 | 777 | trackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3)); |
778 | z3 = trackParamAtHit3->GetHitForRecPtr()->GetZ(); | |
de2cd600 | 779 | (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); |
780 | } | |
781 | // equal contribution from multiple scattering in non bending plane | |
782 | (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1); | |
783 | if (hitNumber1 == hitNumber2) { | |
784 | // Diagonal elements: add contribution from position measurements | |
785 | // in bending plane | |
2457f726 | 786 | (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + hitForRec1->GetBendingReso2(); |
de2cd600 | 787 | // and in non bending plane |
2457f726 | 788 | (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + hitForRec1->GetNonBendingReso2(); |
de2cd600 | 789 | } else { |
790 | // Non diagonal elements: symmetrization | |
791 | // for bending plane | |
792 | (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1); | |
793 | // and non bending plane | |
794 | (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1); | |
795 | } | |
796 | } // for (hitNumber2 = hitNumber1;... | |
797 | } // for (hitNumber1 = 0;... | |
798 | ||
799 | // Inversion of covariance matrices | |
de2cd600 | 800 | Int_t ifailBending; |
208f139e | 801 | gMinuit->mnvert(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending); |
de2cd600 | 802 | Int_t ifailNonBending; |
208f139e | 803 | gMinuit->mnvert(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending); |
de2cd600 | 804 | |
805 | // It would be worth trying to calculate the inverse of the covariance matrix | |
806 | // only once per fit, since it cannot change much in principle, | |
807 | // and it would save a lot of computing time !!!! | |
808 | ||
809 | // Calculates Chi2 | |
810 | if ((ifailBending == 0) && (ifailNonBending == 0)) { | |
811 | // with Multiple Scattering if inversion correct | |
812 | for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { | |
2457f726 | 813 | trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1)); |
814 | hitForRec1 = trackParamAtHit1->GetHitForRecPtr(); | |
815 | hbc1 = hitForRec1->GetBendingCoor(); | |
816 | pbc1 = trackParamAtHit1->GetBendingCoor(); | |
817 | hnbc1 = hitForRec1->GetNonBendingCoor(); | |
818 | pnbc1 = trackParamAtHit1->GetNonBendingCoor(); | |
de2cd600 | 819 | for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) { |
2457f726 | 820 | trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2)); |
821 | hitForRec2 = trackParamAtHit2->GetHitForRecPtr(); | |
822 | hbc2 = hitForRec2->GetBendingCoor(); | |
823 | pbc2 = trackParamAtHit2->GetBendingCoor(); | |
824 | hnbc2 = hitForRec2->GetNonBendingCoor(); | |
825 | pnbc2 = trackParamAtHit2->GetNonBendingCoor(); | |
de2cd600 | 826 | Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) + |
827 | ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2)); | |
828 | } | |
829 | } | |
830 | } else { | |
831 | // without Multiple Scattering if inversion impossible | |
832 | for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { | |
2457f726 | 833 | trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1)); |
834 | hitForRec1 = trackParamAtHit1->GetHitForRecPtr(); | |
835 | hbc1 = hitForRec1->GetBendingCoor(); | |
836 | pbc1 = trackParamAtHit1->GetBendingCoor(); | |
837 | hnbc1 = hitForRec1->GetNonBendingCoor(); | |
838 | pnbc1 = trackParamAtHit1->GetNonBendingCoor(); | |
839 | Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / hitForRec1->GetBendingReso2()) + | |
840 | ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / hitForRec1->GetNonBendingReso2()); | |
de2cd600 | 841 | } |
842 | } | |
843 | ||
844 | delete covBending; | |
845 | delete covNonBending; | |
846 | delete [] msa2; | |
847 | } | |
848 | ||
208f139e | 849 | //__________________________________________________________________________ |
de2cd600 | 850 | Double_t MultipleScatteringAngle2(AliMUONTrackParam *param) |
851 | { | |
2457f726 | 852 | /// Returns square of multiple Coulomb scattering angle |
853 | /// from TrackParamAtHit pointed to by "param" | |
de2cd600 | 854 | Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2; |
855 | Double_t varMultipleScatteringAngle; | |
de2cd600 | 856 | slopeBending = param->GetBendingSlope(); |
857 | slopeNonBending = param->GetNonBendingSlope(); | |
858 | // thickness in radiation length for the current track, | |
859 | // taking local angle into account | |
208f139e | 860 | radiationLength = AliMUONConstants::ChamberThicknessInX0() * |
de2cd600 | 861 | TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending); |
862 | inverseBendingMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum(); | |
863 | inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) / | |
864 | (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); | |
865 | varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength)); | |
866 | // The velocity is assumed to be 1 !!!! | |
867 | varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * varMultipleScatteringAngle * varMultipleScatteringAngle; | |
868 | return varMultipleScatteringAngle; | |
d837040f | 869 | } |
870 | ||
b8dc484b | 871 | //__________________________________________________________________________ |
71a2d3aa | 872 | void AliMUONTrackReconstructor::FillMUONTrack(void) |
b8dc484b | 873 | { |
b58638a9 | 874 | /// Fill AliMUONTrack's fHitForRecAtHit array |
875 | /// Recompute track parameters and covariances at each attached cluster from those at the first one | |
b8dc484b | 876 | AliMUONTrack *track; |
b58638a9 | 877 | AliMUONTrackParam trackParam; |
de2cd600 | 878 | AliMUONTrackParam *trackParamAtHit; |
b58638a9 | 879 | AliMUONHitForRec *hitForRecAtHit; |
880 | ||
b8dc484b | 881 | track = (AliMUONTrack*) fRecTracksPtr->First(); |
882 | while (track) { | |
de2cd600 | 883 | trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); |
b58638a9 | 884 | trackParam = *trackParamAtHit; |
de2cd600 | 885 | while (trackParamAtHit) { |
b58638a9 | 886 | hitForRecAtHit = trackParamAtHit->GetHitForRecPtr(); |
887 | // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit | |
888 | AliMUONTrackExtrap::ExtrapToZCov(&trackParam, hitForRecAtHit->GetZ()); | |
889 | // update track parameters of the current hit | |
890 | trackParamAtHit->SetTrackParam(trackParam); | |
891 | // update covariance matrix of track parameters of the current hit | |
892 | trackParamAtHit->SetCovariances(trackParam.GetCovariances()); | |
893 | // update array of track hit | |
894 | track->AddHitForRecAtHit(hitForRecAtHit); | |
895 | // prepare next step, add MCS effects in parameter covariances | |
896 | AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); | |
de2cd600 | 897 | trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit)); |
898 | } | |
b8dc484b | 899 | track = (AliMUONTrack*) fRecTracksPtr->After(track); |
de2cd600 | 900 | } |
b8dc484b | 901 | return; |
902 | } | |
903 | ||
8429a5e4 | 904 | //__________________________________________________________________________ |
29f1b13a | 905 | void AliMUONTrackReconstructor::EventDump(void) |
04b5ea16 | 906 | { |
2457f726 | 907 | /// Dump reconstructed event (track parameters at vertex and at first hit), |
908 | /// and the particle parameters | |
04b5ea16 | 909 | AliMUONTrack *track; |
8cde4af5 | 910 | AliMUONTrackParam trackParam, *trackParam1; |
04b5ea16 | 911 | Double_t bendingSlope, nonBendingSlope, pYZ; |
912 | Double_t pX, pY, pZ, x, y, z, c; | |
5e671e06 | 913 | Int_t trackIndex, nTrackHits; |
04b5ea16 | 914 | |
b840996b | 915 | AliDebug(1,"****** enter EventDump ******"); |
916 | AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); | |
917 | ||
04b5ea16 | 918 | fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole |
919 | // Loop over reconstructed tracks | |
920 | for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) { | |
b840996b | 921 | AliDebug(1, Form("track number: %d", trackIndex)); |
04b5ea16 | 922 | // function for each track for modularity ???? |
923 | track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]); | |
924 | nTrackHits = track->GetNTrackHits(); | |
b840996b | 925 | AliDebug(1, Form("Number of track hits: %d ", nTrackHits)); |
04b5ea16 | 926 | // track parameters at Vertex |
8cde4af5 | 927 | trackParam = (*((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())); |
928 | AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.); | |
929 | x = trackParam.GetNonBendingCoor(); | |
930 | y = trackParam.GetBendingCoor(); | |
931 | z = trackParam.GetZ(); | |
932 | bendingSlope = trackParam.GetBendingSlope(); | |
933 | nonBendingSlope = trackParam.GetNonBendingSlope(); | |
934 | pYZ = 1/TMath::Abs(trackParam.GetInverseBendingMomentum()); | |
04b5ea16 | 935 | pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope); |
936 | pX = pZ * nonBendingSlope; | |
937 | pY = pZ * bendingSlope; | |
8cde4af5 | 938 | c = TMath::Sign(1.0, trackParam.GetInverseBendingMomentum()); |
b840996b | 939 | AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n", |
940 | z, x, y, pX, pY, pZ, c)); | |
04b5ea16 | 941 | |
942 | // track parameters at first hit | |
de2cd600 | 943 | trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First(); |
04b5ea16 | 944 | x = trackParam1->GetNonBendingCoor(); |
945 | y = trackParam1->GetBendingCoor(); | |
946 | z = trackParam1->GetZ(); | |
947 | bendingSlope = trackParam1->GetBendingSlope(); | |
948 | nonBendingSlope = trackParam1->GetNonBendingSlope(); | |
949 | pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum()); | |
950 | pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope); | |
951 | pX = pZ * nonBendingSlope; | |
952 | pY = pZ * bendingSlope; | |
953 | c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum()); | |
b840996b | 954 | AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n", |
955 | z, x, y, pX, pY, pZ, c)); | |
04b5ea16 | 956 | } |
8f373020 | 957 | // informations about generated particles NO !!!!!!!! |
04b5ea16 | 958 | |
88cb7938 | 959 | // for (Int_t iPart = 0; iPart < np; iPart++) { |
960 | // p = gAlice->Particle(iPart); | |
961 | // printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n", | |
962 | // iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode()); | |
963 | // } | |
a9e2aefa | 964 | return; |
965 | } | |
966 | ||
276c44b7 | 967 |