]>
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" |
2457f726 | 31 | #include "AliMUONHitForRec.h" |
a9e2aefa | 32 | #include "AliMUONTrack.h" |
37827b29 | 33 | #include "AliMUONTrackParam.h" |
34 | #include "AliMUONTrackExtrap.h" | |
8c343c7c | 35 | #include "AliLog.h" |
af743f54 | 36 | |
208f139e | 37 | #include <TMinuit.h> |
af743f54 | 38 | #include <Riostream.h> |
ea94c18b | 39 | #include <TMath.h> |
af743f54 | 40 | #include <TMatrixD.h> |
a9e2aefa | 41 | |
de2cd600 | 42 | // Functions to be minimized with Minuit |
ea94c18b | 43 | void TrackChi2(Int_t &nParam, Double_t *gradient, Double_t &chi2, Double_t *param, Int_t flag); |
de2cd600 | 44 | |
7945aae7 | 45 | /// \cond CLASSIMP |
29f1b13a | 46 | ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context |
7945aae7 | 47 | /// \endcond |
a9e2aefa | 48 | |
ea94c18b | 49 | //************* Parameters for reconstruction |
50 | const Double_t AliMUONTrackReconstructor::fgkBendingVertexDispersion = 10.; | |
51 | const Double_t AliMUONTrackReconstructor::fgkNonBendingVertexDispersion = 10.; | |
52 | ||
2457f726 | 53 | |
3bc8b580 | 54 | //__________________________________________________________________________ |
7ec3b9cf | 55 | AliMUONTrackReconstructor::AliMUONTrackReconstructor() |
56 | : AliMUONVTrackReconstructor() | |
a9e2aefa | 57 | { |
2457f726 | 58 | /// Constructor for class AliMUONTrackReconstructor |
ea94c18b | 59 | AliInfo("*** Original tracking ***"); |
a9e2aefa | 60 | } |
9cbdf048 | 61 | |
7ec3b9cf | 62 | //__________________________________________________________________________ |
ea94c18b | 63 | AliMUONTrackReconstructor::~AliMUONTrackReconstructor() |
a9e2aefa | 64 | { |
ea94c18b | 65 | /// Destructor |
66 | } | |
276c44b7 | 67 | |
de2cd600 | 68 | //__________________________________________________________________________ |
ea94c18b | 69 | void AliMUONTrackReconstructor::MakeTrackCandidates() |
de2cd600 | 70 | { |
208f139e | 71 | /// To make track candidates: |
72 | /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4. | |
73 | /// Good candidates are made of at least three hitForRec's. | |
74 | /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks. | |
75 | TClonesArray *segments; | |
208f139e | 76 | AliMUONTrack *track; |
208f139e | 77 | Int_t iCandidate = 0; |
78 | ||
de2cd600 | 79 | AliDebug(1,"Enter MakeTrackCandidates"); |
de2cd600 | 80 | |
208f139e | 81 | // Loop over stations(1..) 5 and 4 and make track candidates |
ea94c18b | 82 | for (Int_t istat=4; istat>=3; istat--) { |
83 | ||
208f139e | 84 | // Make segments in the station |
85 | segments = MakeSegmentsInStation(istat); | |
ea94c18b | 86 | |
208f139e | 87 | // Loop over segments |
7ec3b9cf | 88 | for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) |
89 | { | |
208f139e | 90 | AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate)); |
ea94c18b | 91 | |
208f139e | 92 | // Transform segments to tracks and put them at the end of fRecTracksPtr |
ea94c18b | 93 | track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg])); |
a9e2aefa | 94 | fNRecTracks++; |
ea94c18b | 95 | |
208f139e | 96 | // Printout for debuging |
7ec3b9cf | 97 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) |
98 | { | |
208f139e | 99 | cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl; |
ea94c18b | 100 | ((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())->GetCovariances().Print(); |
208f139e | 101 | } |
ea94c18b | 102 | |
208f139e | 103 | // Look for compatible hitForRec(s) in the other station |
ea94c18b | 104 | if (!FollowTrackInStation(*track,7-istat)) { |
105 | fRecTracksPtr->Remove(track); | |
106 | fNRecTracks--; | |
107 | } | |
108 | ||
a9e2aefa | 109 | } |
ea94c18b | 110 | |
208f139e | 111 | // delete the array of segments |
112 | delete segments; | |
113 | } | |
ea94c18b | 114 | |
208f139e | 115 | fRecTracksPtr->Compress(); // this is essential before checking tracks |
116 | ||
117 | // Keep all different tracks or only the best ones as required | |
118 | if (fgkTrackAllTracks) RemoveIdenticalTracks(); | |
119 | else RemoveDoubleTracks(); | |
120 | ||
121 | AliDebug(1,Form("Number of good candidates = %d",fNRecTracks)); | |
122 | ||
a9e2aefa | 123 | } |
124 | ||
125 | //__________________________________________________________________________ | |
ea94c18b | 126 | void AliMUONTrackReconstructor::FollowTracks() |
a9e2aefa | 127 | { |
2457f726 | 128 | /// Follow tracks in stations(1..) 3, 2 and 1 |
208f139e | 129 | AliDebug(1,"Enter FollowTracks"); |
130 | ||
04b5ea16 | 131 | AliMUONTrack *track, *nextTrack; |
af743f54 | 132 | Int_t currentNRecTracks; |
ea94c18b | 133 | Bool_t hitFound; |
208f139e | 134 | |
135 | for (Int_t station = 2; station >= 0; station--) { | |
ea94c18b | 136 | |
208f139e | 137 | // Save the actual number of reconstructed track in case of |
138 | // tracks are added or suppressed during the tracking procedure | |
139 | // !! Do not compress fRecTracksPtr until the end of the loop over tracks !! | |
af743f54 | 140 | currentNRecTracks = fNRecTracks; |
ea94c18b | 141 | |
af743f54 | 142 | for (Int_t iRecTrack = 0; iRecTrack <currentNRecTracks; iRecTrack++) { |
208f139e | 143 | AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1)); |
ea94c18b | 144 | |
208f139e | 145 | track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack); |
ea94c18b | 146 | |
208f139e | 147 | // Fit the track: |
148 | // Do not take into account the multiple scattering to speed up the fit | |
149 | // Calculate the track parameter covariance matrix | |
150 | // If "station" is station(1..) 3 then use the vertex to better constrain the fit | |
151 | if (station==2) { | |
ea94c18b | 152 | SetVertexForFit(*track); |
208f139e | 153 | track->SetFitWithVertex(kTRUE); |
154 | } else track->SetFitWithVertex(kFALSE); | |
ea94c18b | 155 | Fit(*track, kFALSE, kTRUE); |
156 | ||
208f139e | 157 | // Remove the track if the normalized chi2 is too high |
ea94c18b | 158 | if (track->GetNormalizedChi2() > fgkSigmaToCutForTracking * fgkSigmaToCutForTracking) { |
208f139e | 159 | fRecTracksPtr->Remove(track); |
160 | fNRecTracks--; | |
161 | continue; | |
162 | } | |
ea94c18b | 163 | |
208f139e | 164 | // Printout for debuging |
165 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
166 | cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl; | |
ea94c18b | 167 | ((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())->GetCovariances().Print(); |
208f139e | 168 | } |
ea94c18b | 169 | |
208f139e | 170 | // Look for compatible hitForRec in station(0..) "station" |
ea94c18b | 171 | hitFound = FollowTrackInStation(*track,station); |
172 | ||
173 | // Try to recover track if required | |
174 | if (!hitFound && fgkRecoverTracks) hitFound = RecoverTrack(*track,station); | |
175 | ||
176 | // remove track if no hit found | |
177 | if (!hitFound) { | |
178 | fRecTracksPtr->Remove(track); | |
179 | fNRecTracks--; | |
180 | } | |
181 | ||
208f139e | 182 | } |
ea94c18b | 183 | |
208f139e | 184 | // Compress fRecTracksPtr for the next step |
185 | fRecTracksPtr->Compress(); | |
ea94c18b | 186 | |
208f139e | 187 | // Keep only the best tracks if required |
188 | if (!fgkTrackAllTracks) RemoveDoubleTracks(); | |
ea94c18b | 189 | |
208f139e | 190 | } |
191 | ||
192 | // Last fit of track candidates with all station | |
193 | // Take into account the multiple scattering and remove bad tracks | |
194 | Int_t trackIndex = -1; | |
04b5ea16 | 195 | track = (AliMUONTrack*) fRecTracksPtr->First(); |
04b5ea16 | 196 | while (track) { |
ea94c18b | 197 | |
04b5ea16 | 198 | trackIndex++; |
199 | nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track | |
ea94c18b | 200 | |
208f139e | 201 | track->SetFitWithVertex(kFALSE); // just to be sure |
ea94c18b | 202 | Fit(*track, kTRUE, kTRUE); |
203 | ||
208f139e | 204 | // Printout for debuging |
205 | if (AliLog::GetGlobalDebugLevel() >= 3) { | |
206 | cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl; | |
a9e2aefa | 207 | track->RecursiveDump(); |
208f139e | 208 | } |
ea94c18b | 209 | |
208f139e | 210 | // Remove the track if the normalized chi2 is too high |
ea94c18b | 211 | if (track->GetNormalizedChi2() > fgkSigmaToCutForTracking * fgkSigmaToCutForTracking) { |
208f139e | 212 | fRecTracksPtr->Remove(track); |
213 | fNRecTracks--; | |
a9e2aefa | 214 | } |
ea94c18b | 215 | |
208f139e | 216 | track = nextTrack; |
ea94c18b | 217 | |
208f139e | 218 | } |
ea94c18b | 219 | |
208f139e | 220 | fRecTracksPtr->Compress(); |
221 | ||
222 | } | |
1a38e749 | 223 | |
208f139e | 224 | //__________________________________________________________________________ |
ea94c18b | 225 | Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, Int_t nextStation) |
208f139e | 226 | { |
227 | /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s) | |
228 | /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks: | |
229 | /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of | |
230 | /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure. | |
231 | /// Remove the obsolete "trackCandidate" at the end. | |
232 | /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority. | |
233 | AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1)); | |
234 | ||
ea94c18b | 235 | // Order the chamber according to the propagation direction (tracking starts with chamber 2): |
236 | // - nextStation == station(1...) 5 => forward propagation | |
237 | // - nextStation < station(1...) 5 => backward propagation | |
238 | Int_t ch1, ch2; | |
239 | if (nextStation==4) { | |
240 | ch1 = 2*nextStation+1; | |
241 | ch2 = 2*nextStation; | |
242 | } else { | |
243 | ch1 = 2*nextStation; | |
244 | ch2 = 2*nextStation+1; | |
245 | } | |
246 | ||
208f139e | 247 | Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2); |
248 | Double_t chi2WithOneHitForRec = 1.e10; | |
249 | Double_t chi2WithTwoHitForRec = 1.e10; | |
ea94c18b | 250 | Double_t maxChi2WithOneHitForRec = 2. * fgkSigmaToCutForTracking * fgkSigmaToCutForTracking; // 2 because 2 quantities in chi2 |
251 | Double_t maxChi2WithTwoHitForRec = 4. * fgkSigmaToCutForTracking * fgkSigmaToCutForTracking; // 4 because 4 quantities in chi2 | |
208f139e | 252 | Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec; |
253 | Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec; | |
ea94c18b | 254 | Bool_t foundOneHit = kFALSE; |
255 | Bool_t foundTwoHits = kFALSE; | |
208f139e | 256 | AliMUONTrack *newTrack = 0x0; |
257 | AliMUONHitForRec *hitForRecCh1, *hitForRecCh2; | |
ea94c18b | 258 | AliMUONTrackParam extrapTrackParamAtHit1; |
259 | AliMUONTrackParam extrapTrackParamAtHit2; | |
260 | AliMUONTrackParam bestTrackParamAtHit1; | |
261 | AliMUONTrackParam bestTrackParamAtHit2; | |
847cbaef | 262 | Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]]; |
263 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE; | |
ea94c18b | 264 | |
265 | // Get track parameters and extrapolate them to chamber "ch2" to save computing time in the next steps | |
266 | AliMUONTrackParam extrapTrackParamAtCh2(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->First()); | |
267 | ||
268 | // Add MCS effect | |
269 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
270 | ||
271 | // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria) | |
272 | if (ch1 < ch2 && extrapTrackParamAtCh2.GetHitForRecPtr()->GetChamberNumber() > ch2 + 1) { | |
273 | // extrapolation to the missing chamber | |
274 | AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh2, AliMUONConstants::DefaultChamberZ(ch2 + 1)); | |
275 | // add MCS effect | |
276 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
277 | } | |
278 | ||
279 | //Extrapolate trackCandidate to chamber "ch2" | |
280 | AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh2, zCh2); | |
281 | ||
208f139e | 282 | // Printout for debuging |
283 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { | |
208f139e | 284 | cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl; |
ea94c18b | 285 | extrapTrackParamAtCh2.GetCovariances().Print(); |
208f139e | 286 | } |
ea94c18b | 287 | |
208f139e | 288 | // Printout for debuging |
289 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
290 | cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl; | |
291 | } | |
ea94c18b | 292 | |
293 | // look for candidates in chamber 2 | |
208f139e | 294 | for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) { |
ea94c18b | 295 | |
208f139e | 296 | hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2); |
ea94c18b | 297 | |
298 | // try to add the current hit | |
299 | chi2WithOneHitForRec = TryOneHitForRec(extrapTrackParamAtCh2, hitForRecCh2, extrapTrackParamAtHit2); | |
300 | ||
208f139e | 301 | // if good chi2 then try to attach a hitForRec in the other chamber too |
302 | if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) { | |
ea94c18b | 303 | Bool_t foundSecondHit = kFALSE; |
304 | ||
208f139e | 305 | // Printout for debuging |
306 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
ea94c18b | 307 | cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch2+1 |
308 | << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl; | |
309 | cout << " look for second hits in chamber(1..): " << ch1+1 << " ..." << endl; | |
a9e2aefa | 310 | } |
ea94c18b | 311 | |
312 | // add MCS effect for next step | |
313 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtHit2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
314 | ||
208f139e | 315 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { |
ea94c18b | 316 | |
317 | hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); | |
318 | ||
319 | // try to add the current hit in addition to the one found on the previous chamber | |
320 | chi2WithTwoHitForRec = TryTwoHitForRec(extrapTrackParamAtHit2, hitForRecCh1, extrapTrackParamAtHit1); | |
321 | ||
322 | // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate" | |
208f139e | 323 | if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) { |
324 | foundSecondHit = kTRUE; | |
ea94c18b | 325 | foundTwoHits = kTRUE; |
326 | ||
327 | // Printout for debuging | |
328 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
329 | cout << "FollowTrackInStation: found second hit in chamber(1..): " << ch1+1 | |
330 | << " (Global Chi2 = " << chi2WithTwoHitForRec << ")" << endl; | |
331 | } | |
332 | ||
333 | if (fgkTrackAllTracks) { | |
208f139e | 334 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's |
ea94c18b | 335 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); |
336 | UpdateTrack(*newTrack,extrapTrackParamAtHit1,extrapTrackParamAtHit2); | |
208f139e | 337 | fNRecTracks++; |
ea94c18b | 338 | |
847cbaef | 339 | // Tag hitForRecCh1 as used |
340 | hitForRecCh1Used[hit1] = kTRUE; | |
ea94c18b | 341 | |
208f139e | 342 | // Printout for debuging |
343 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
ea94c18b | 344 | cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1 << endl; |
208f139e | 345 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); |
a9e2aefa | 346 | } |
ea94c18b | 347 | |
208f139e | 348 | } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) { |
349 | // keep track of the best couple of hits | |
350 | bestChi2WithTwoHitForRec = chi2WithTwoHitForRec; | |
ea94c18b | 351 | bestTrackParamAtHit1 = extrapTrackParamAtHit1; |
352 | bestTrackParamAtHit2 = extrapTrackParamAtHit2; | |
208f139e | 353 | } |
ea94c18b | 354 | |
a9e2aefa | 355 | } |
ea94c18b | 356 | |
d0bfce8d | 357 | } |
ea94c18b | 358 | |
208f139e | 359 | // if no hitForRecCh1 found then consider to add hitForRecCh2 only |
360 | if (!foundSecondHit) { | |
ea94c18b | 361 | foundOneHit = kTRUE; |
362 | ||
363 | if (fgkTrackAllTracks) { | |
208f139e | 364 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's |
ea94c18b | 365 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); |
366 | UpdateTrack(*newTrack,extrapTrackParamAtHit2); | |
208f139e | 367 | fNRecTracks++; |
ea94c18b | 368 | |
208f139e | 369 | // Printout for debuging |
370 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
ea94c18b | 371 | cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1 << endl; |
208f139e | 372 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); |
373 | } | |
ea94c18b | 374 | |
375 | } else if (!foundTwoHits && chi2WithOneHitForRec < bestChi2WithOneHitForRec) { | |
208f139e | 376 | // keep track of the best single hitForRec except if a couple |
377 | // of hits has already been found (i.e. bestHitForRec2!=0x0) | |
378 | bestChi2WithOneHitForRec = chi2WithOneHitForRec; | |
ea94c18b | 379 | bestTrackParamAtHit1 = extrapTrackParamAtHit2; |
208f139e | 380 | } |
ea94c18b | 381 | |
d0bfce8d | 382 | } |
ea94c18b | 383 | |
208f139e | 384 | } |
ea94c18b | 385 | |
208f139e | 386 | } |
ea94c18b | 387 | |
208f139e | 388 | // look for candidates in chamber 1 not already attached to a track |
389 | // if we want to keep all possible tracks or if no good couple of hitForRec has been found | |
ea94c18b | 390 | if (fgkTrackAllTracks || !foundTwoHits) { |
391 | ||
392 | // Printout for debuging | |
208f139e | 393 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { |
394 | cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl; | |
395 | } | |
ea94c18b | 396 | |
397 | // add MCS effect for next step | |
398 | AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh2,AliMUONConstants::ChamberThicknessInX0(),1.); | |
399 | ||
208f139e | 400 | for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { |
ea94c18b | 401 | |
208f139e | 402 | hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); |
ea94c18b | 403 | |
847cbaef | 404 | if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used |
ea94c18b | 405 | |
406 | // try to add the current hit | |
407 | chi2WithOneHitForRec = TryOneHitForRec(extrapTrackParamAtCh2, hitForRecCh1, extrapTrackParamAtHit1); | |
408 | ||
409 | // if good chi2 then consider to add hitForRecCh1 | |
208f139e | 410 | // We do not try to attach a hitForRec in the other chamber too since it has already been done above |
411 | if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) { | |
ea94c18b | 412 | foundOneHit = kTRUE; |
413 | ||
414 | // Printout for debuging | |
415 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
416 | cout << "FollowTrackInStation: found one hit in chamber(1..): " << ch1+1 | |
417 | << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl; | |
418 | } | |
419 | ||
208f139e | 420 | if (fgkTrackAllTracks) { |
421 | // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's | |
ea94c18b | 422 | newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); |
423 | UpdateTrack(*newTrack,extrapTrackParamAtHit1); | |
208f139e | 424 | fNRecTracks++; |
ea94c18b | 425 | |
426 | // Printout for debuging | |
208f139e | 427 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { |
ea94c18b | 428 | cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1 << endl; |
208f139e | 429 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); |
430 | } | |
ea94c18b | 431 | |
432 | } else if (chi2WithOneHitForRec < bestChi2WithOneHitForRec) { | |
208f139e | 433 | // keep track of the best single hitForRec except if a couple |
434 | // of hits has already been found (i.e. bestHitForRec1!=0x0) | |
435 | bestChi2WithOneHitForRec = chi2WithOneHitForRec; | |
ea94c18b | 436 | bestTrackParamAtHit1 = extrapTrackParamAtHit1; |
208f139e | 437 | } |
ea94c18b | 438 | |
d0bfce8d | 439 | } |
ea94c18b | 440 | |
208f139e | 441 | } |
ea94c18b | 442 | |
208f139e | 443 | } |
ea94c18b | 444 | |
208f139e | 445 | // fill out the best track if required else clean up the fRecTracksPtr array |
ea94c18b | 446 | if (!fgkTrackAllTracks) { |
447 | if (foundTwoHits) { | |
448 | UpdateTrack(trackCandidate,bestTrackParamAtHit1,bestTrackParamAtHit2); | |
449 | ||
208f139e | 450 | // Printout for debuging |
451 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
ea94c18b | 452 | cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1 << endl; |
208f139e | 453 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); |
a9e2aefa | 454 | } |
ea94c18b | 455 | |
456 | } else if (foundOneHit) { | |
457 | UpdateTrack(trackCandidate,bestTrackParamAtHit1); | |
458 | ||
208f139e | 459 | // Printout for debuging |
460 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
ea94c18b | 461 | cout << "FollowTrackInStation: added the best hit in chamber(1..): " << bestTrackParamAtHit1.GetHitForRecPtr()->GetChamberNumber()+1 << endl; |
208f139e | 462 | if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); |
04b5ea16 | 463 | } |
ea94c18b | 464 | |
465 | } else return kFALSE; | |
466 | ||
467 | } else if (foundOneHit || foundTwoHits) { | |
468 | ||
469 | // remove obsolete track | |
470 | fRecTracksPtr->Remove(&trackCandidate); | |
471 | fNRecTracks--; | |
472 | ||
473 | } else return kFALSE; | |
474 | ||
475 | return kTRUE; | |
476 | ||
477 | } | |
478 | ||
479 | //__________________________________________________________________________ | |
480 | Double_t AliMUONTrackReconstructor::TryTwoHitForRec(const AliMUONTrackParam &trackParamAtHit1, AliMUONHitForRec* hitForRec2, AliMUONTrackParam &trackParamAtHit2) | |
481 | { | |
482 | /// Test the compatibility between the track and the 2 hitForRec together (using trackParam's covariance matrix): | |
483 | /// return the corresponding Chi2 accounting for covariances between the 2 hitForRec | |
484 | /// return trackParamAtHit1 & 2 | |
485 | ||
486 | if (!trackParamAtHit1.CovariancesExist()) AliWarning(" track parameter covariance matrix does not exist"); | |
487 | AliMUONHitForRec* hitForRec1 = trackParamAtHit1.GetHitForRecPtr(); | |
488 | if (!hitForRec1) AliFatal(" no hitForRec attached to trackParamAtHit1"); | |
489 | ||
490 | // extrapolate track parameters at the z position of the second hit | |
491 | trackParamAtHit2.SetParameters(trackParamAtHit1.GetParameters()); | |
492 | trackParamAtHit2.SetZ(trackParamAtHit1.GetZ()); | |
493 | AliMUONTrackExtrap::ExtrapToZ(&trackParamAtHit2, hitForRec2->GetZ()); | |
494 | ||
495 | // set pointer to hit2 into trackParamAtHit2 | |
496 | trackParamAtHit2.SetHitForRecPtr(hitForRec2); | |
497 | ||
498 | // Set differences between track and the 2 hitForRec in the bending and non bending directions | |
499 | TMatrixD dPos(4,1); | |
500 | dPos(0,0) = hitForRec1->GetNonBendingCoor() - trackParamAtHit1.GetNonBendingCoor(); | |
501 | dPos(1,0) = hitForRec1->GetBendingCoor() - trackParamAtHit1.GetBendingCoor(); | |
502 | dPos(2,0) = hitForRec2->GetNonBendingCoor() - trackParamAtHit2.GetNonBendingCoor(); | |
503 | dPos(3,0) = hitForRec2->GetBendingCoor() - trackParamAtHit2.GetBendingCoor(); | |
504 | ||
505 | // quick tests of hitForRec compatibility within a wide road of x*y = 1*1 cm2 to save computing time | |
506 | if (TMath::Abs(dPos(0,0)) > fgkMaxTrackingDistanceNonBending || | |
507 | TMath::Abs(dPos(1,0)) > fgkMaxTrackingDistanceBending || | |
508 | TMath::Abs(dPos(2,0)) > fgkMaxTrackingDistanceNonBending || | |
509 | TMath::Abs(dPos(3,0)) > fgkMaxTrackingDistanceBending) return 1.e10; | |
510 | ||
511 | // Calculate the error matrix from the track parameter covariances at first hitForRec | |
512 | TMatrixD error(4,4); | |
513 | error.Zero(); | |
514 | if (trackParamAtHit1.CovariancesExist()) { | |
515 | // Save track parameters at first hitForRec | |
516 | AliMUONTrackParam trackParamAtHit1Save(trackParamAtHit1); | |
517 | TMatrixD paramAtHit1Save(trackParamAtHit1Save.GetParameters()); | |
518 | Double_t z1 = trackParamAtHit1Save.GetZ(); | |
519 | ||
520 | // Save track coordinates at second hitForRec | |
521 | Double_t nonBendingCoor2 = trackParamAtHit2.GetNonBendingCoor(); | |
522 | Double_t bendingCoor2 = trackParamAtHit2.GetBendingCoor(); | |
523 | ||
524 | // add MCS effect at first hitForRec | |
525 | AliMUONTrackExtrap::AddMCSEffect(&trackParamAtHit1Save,AliMUONConstants::ChamberThicknessInX0(),1.); | |
526 | ||
527 | // Get the pointer to the parameter covariance matrix at first hitForRec | |
528 | const TMatrixD& kParamCov = trackParamAtHit1Save.GetCovariances(); | |
529 | ||
530 | // Calculate the jacobian related to the transformation between track parameters | |
531 | // at first hitForRec and track coordinates at the 2 hitForRec z-position | |
532 | TMatrixD jacob(4,5); | |
533 | jacob.Zero(); | |
534 | // first derivative at the first hitForRec: | |
535 | jacob(0,0) = 1.; // dx1/dx | |
536 | jacob(1,2) = 1.; // dy1/dy | |
537 | // first derivative at the second hitForRec: | |
538 | TMatrixD dParam(5,1); | |
539 | for (Int_t i=0; i<5; i++) { | |
540 | // Skip jacobian calculation for parameters with no associated error | |
541 | if (kParamCov(i,i) == 0.) continue; | |
542 | // Small variation of parameter i only | |
543 | for (Int_t j=0; j<5; j++) { | |
544 | if (j==i) { | |
545 | dParam(j,0) = TMath::Sqrt(kParamCov(i,i)); | |
546 | if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramAtHit1Save(4,0)); // variation always in the same direction | |
547 | } else dParam(j,0) = 0.; | |
548 | } | |
549 | ||
550 | // Set new track parameters at first hitForRec | |
551 | trackParamAtHit1Save.SetParameters(paramAtHit1Save); | |
552 | trackParamAtHit1Save.AddParameters(dParam); | |
553 | trackParamAtHit1Save.SetZ(z1); | |
554 | ||
555 | // Extrapolate new track parameters to the z position of the second hitForRec | |
556 | AliMUONTrackExtrap::ExtrapToZ(&trackParamAtHit1Save,hitForRec2->GetZ()); | |
557 | ||
558 | // Calculate the jacobian | |
559 | jacob(2,i) = (trackParamAtHit1Save.GetNonBendingCoor() - nonBendingCoor2) / dParam(i,0); // dx2/dParami | |
560 | jacob(3,i) = (trackParamAtHit1Save.GetBendingCoor() - bendingCoor2 ) / dParam(i,0); // dy2/dParami | |
208f139e | 561 | } |
ea94c18b | 562 | |
563 | // Calculate the error matrix | |
564 | TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov); | |
565 | error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob); | |
566 | } | |
567 | ||
568 | // Add hitForRec resolution to the error matrix | |
569 | error(0,0) += hitForRec1->GetNonBendingReso2(); | |
570 | error(1,1) += hitForRec1->GetBendingReso2(); | |
571 | error(2,2) += hitForRec2->GetNonBendingReso2(); | |
572 | error(3,3) += hitForRec2->GetBendingReso2(); | |
573 | ||
574 | // invert the error matrix for Chi2 calculation | |
575 | if (error.Determinant() != 0) { | |
576 | error.Invert(); | |
208f139e | 577 | } else { |
ea94c18b | 578 | AliWarning(" Determinant error=0"); |
579 | return 1.e10; | |
208f139e | 580 | } |
581 | ||
ea94c18b | 582 | // Compute the Chi2 value |
583 | TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error); | |
584 | TMatrixD result(tmp2,TMatrixD::kMult,dPos); | |
585 | ||
586 | return result(0,0); | |
587 | ||
588 | } | |
589 | ||
590 | //__________________________________________________________________________ | |
591 | void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit) | |
592 | { | |
593 | /// Add 1 hit to the track candidate | |
594 | /// Update chi2 of the track | |
595 | ||
596 | // Compute local chi2 | |
597 | AliMUONHitForRec* hit = trackParamAtHit.GetHitForRecPtr(); | |
598 | Double_t deltaX = trackParamAtHit.GetNonBendingCoor() - hit->GetNonBendingCoor(); | |
599 | Double_t deltaY = trackParamAtHit.GetBendingCoor() - hit->GetBendingCoor(); | |
600 | Double_t localChi2 = deltaX*deltaX / hit->GetNonBendingReso2() + | |
601 | deltaY*deltaY / hit->GetBendingReso2(); | |
602 | ||
603 | // Update the chi2 of the new track | |
604 | track.SetFitFMin(track.GetFitFMin() + localChi2); | |
605 | ||
606 | // Update TrackParamAtHit | |
607 | track.AddTrackParamAtHit(&trackParamAtHit,trackParamAtHit.GetHitForRecPtr()); | |
608 | track.GetTrackParamAtHit()->Sort(); | |
609 | ||
04b5ea16 | 610 | } |
611 | ||
de2cd600 | 612 | //__________________________________________________________________________ |
ea94c18b | 613 | void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtHit1, AliMUONTrackParam &trackParamAtHit2) |
614 | { | |
615 | /// Add 2 hits to the track candidate | |
616 | /// Update track and local chi2 | |
617 | ||
618 | // Update local chi2 at first hit | |
619 | AliMUONHitForRec* hit1 = trackParamAtHit1.GetHitForRecPtr(); | |
620 | Double_t deltaX = trackParamAtHit1.GetNonBendingCoor() - hit1->GetNonBendingCoor(); | |
621 | Double_t deltaY = trackParamAtHit1.GetBendingCoor() - hit1->GetBendingCoor(); | |
622 | Double_t localChi2AtHit1 = deltaX*deltaX / hit1->GetNonBendingReso2() + | |
623 | deltaY*deltaY / hit1->GetBendingReso2(); | |
624 | trackParamAtHit1.SetLocalChi2(localChi2AtHit1); | |
625 | ||
626 | // Flag first hit as being removable | |
627 | trackParamAtHit1.SetRemovable(kTRUE); | |
628 | ||
629 | // Update local chi2 at second hit | |
630 | AliMUONHitForRec* hit2 = trackParamAtHit2.GetHitForRecPtr(); | |
631 | deltaX = trackParamAtHit2.GetNonBendingCoor() - hit2->GetNonBendingCoor(); | |
632 | deltaY = trackParamAtHit2.GetBendingCoor() - hit2->GetBendingCoor(); | |
633 | Double_t localChi2AtHit2 = deltaX*deltaX / hit2->GetNonBendingReso2() + | |
634 | deltaY*deltaY / hit2->GetBendingReso2(); | |
635 | trackParamAtHit2.SetLocalChi2(localChi2AtHit2); | |
636 | ||
637 | // Flag first hit as being removable | |
638 | trackParamAtHit2.SetRemovable(kTRUE); | |
639 | ||
640 | // Update the chi2 of the new track | |
641 | track.SetFitFMin(track.GetFitFMin() + localChi2AtHit1 + localChi2AtHit2); | |
642 | ||
643 | // Update TrackParamAtHit | |
644 | track.AddTrackParamAtHit(&trackParamAtHit1,trackParamAtHit1.GetHitForRecPtr()); | |
645 | track.AddTrackParamAtHit(&trackParamAtHit2,trackParamAtHit2.GetHitForRecPtr()); | |
646 | track.GetTrackParamAtHit()->Sort(); | |
647 | ||
648 | } | |
649 | ||
650 | //__________________________________________________________________________ | |
651 | Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, Int_t nextStation) | |
652 | { | |
653 | /// Try to recover the track candidate in the next station | |
654 | /// by removing the worst of the two hits attached in the current station | |
655 | /// Return kTRUE if recovering succeeds | |
656 | AliDebug(1,"Enter RecoverTrack"); | |
657 | ||
658 | // Do not try to recover track until we have attached hit(s) on station(1..) 3 | |
659 | if (nextStation > 1) return kFALSE; | |
660 | ||
661 | Int_t worstHitNumber = -1; | |
662 | Double_t localChi2, worstLocalChi2 = 0.; | |
663 | ||
664 | // Look for the hit to remove | |
665 | for (Int_t hitNumber = 0; hitNumber < 2; hitNumber++) { | |
666 | AliMUONTrackParam *trackParamAtHit = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(hitNumber); | |
667 | ||
668 | // check if current hit is removable | |
669 | if (!trackParamAtHit->IsRemovable()) return kFALSE; | |
670 | ||
671 | // Pick up hit with the worst chi2 | |
672 | localChi2 = trackParamAtHit->GetLocalChi2(); | |
673 | if (localChi2 > worstLocalChi2) { | |
674 | worstLocalChi2 = localChi2; | |
675 | worstHitNumber = hitNumber; | |
676 | } | |
677 | } | |
678 | ||
679 | // Reset best hit as being NOT removable | |
680 | ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt((worstHitNumber+1)%2))->SetRemovable(kFALSE); | |
681 | ||
682 | // Remove the worst hit | |
683 | trackCandidate.RemoveTrackParamAtHit((AliMUONTrackParam*)trackCandidate.GetTrackParamAtHit()->UncheckedAt(worstHitNumber)); | |
684 | ||
685 | // Re-fit the track: | |
686 | // Do not take into account the multiple scattering to speed up the fit | |
687 | // Calculate the track parameter covariance matrix | |
688 | trackCandidate.SetFitWithVertex(kFALSE); // To be sure | |
689 | Fit(trackCandidate, kFALSE, kTRUE); | |
690 | ||
691 | // Look for new hit(s) in next station | |
692 | return FollowTrackInStation(trackCandidate,nextStation); | |
693 | ||
694 | } | |
695 | ||
696 | //__________________________________________________________________________ | |
697 | void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack &trackCandidate) | |
de2cd600 | 698 | { |
208f139e | 699 | /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate" |
700 | /// Compute the vertex resolution from natural vertex dispersion and | |
701 | /// multiple scattering effets according to trackCandidate path in absorber | |
702 | /// It is necessary to account for multiple scattering effects here instead of during the fit of | |
703 | /// the "trackCandidate" to do not influence the result by changing track resolution at vertex | |
704 | AliDebug(1,"Enter SetVertexForFit"); | |
de2cd600 | 705 | |
ea94c18b | 706 | Double_t nonBendingReso2 = fgkNonBendingVertexDispersion * fgkNonBendingVertexDispersion; |
707 | Double_t bendingReso2 = fgkBendingVertexDispersion * fgkBendingVertexDispersion; | |
208f139e | 708 | // add multiple scattering effets |
ea94c18b | 709 | AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtHit()->First()))); |
208f139e | 710 | paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering |
711 | AliMUONTrackExtrap::ExtrapToVertexUncorrected(¶mAtVertex,0.); | |
ea94c18b | 712 | const TMatrixD& kParamCov = paramAtVertex.GetCovariances(); |
713 | nonBendingReso2 += kParamCov(0,0); | |
714 | bendingReso2 += kParamCov(2,2); | |
208f139e | 715 | // Set the vertex |
716 | AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default | |
717 | vertex.SetNonBendingReso2(nonBendingReso2); | |
718 | vertex.SetBendingReso2(bendingReso2); | |
ea94c18b | 719 | trackCandidate.SetVertex(&vertex); |
208f139e | 720 | } |
721 | ||
722 | //__________________________________________________________________________ | |
ea94c18b | 723 | void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t calcCov) |
208f139e | 724 | { |
ea94c18b | 725 | /// Fit the track "track" w/wo multiple Coulomb scattering according to "includeMCS". |
de2cd600 | 726 | |
208f139e | 727 | Double_t benC, errorParam, invBenP, nonBenC, x, y; |
de2cd600 | 728 | AliMUONTrackParam *trackParam; |
208f139e | 729 | Double_t arg[1], fedm, errdef, fitFMin; |
730 | Int_t npari, nparx; | |
731 | Int_t status, covStatus; | |
732 | ||
8cde4af5 | 733 | // Instantiate gMinuit if not already done |
734 | if (!gMinuit) gMinuit = new TMinuit(6); | |
208f139e | 735 | // Clear MINUIT parameters |
736 | gMinuit->mncler(); | |
737 | // Give the fitted track to MINUIT | |
ea94c18b | 738 | gMinuit->SetObjectFit(&track); |
208f139e | 739 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { |
740 | // Define print level | |
741 | arg[0] = 1; | |
742 | gMinuit->mnexcm("SET PRI", arg, 1, status); | |
743 | // Print covariance matrix | |
744 | gMinuit->mnexcm("SHO COV", arg, 0, status); | |
745 | } else { | |
746 | arg[0] = -1; | |
747 | gMinuit->mnexcm("SET PRI", arg, 1, status); | |
748 | } | |
de2cd600 | 749 | // No warnings |
208f139e | 750 | gMinuit->mnexcm("SET NOW", arg, 0, status); |
751 | // Define strategy | |
752 | //arg[0] = 2; | |
753 | //gMinuit->mnexcm("SET STR", arg, 1, status); | |
754 | ||
ea94c18b | 755 | // set flag w/wo multiple scattering according to "includeMCS" |
756 | track.SetFitWithMCS(includeMCS); | |
757 | if (includeMCS) { | |
758 | // compute hit weights only once | |
759 | if (!track.ComputeHitWeights()) { | |
760 | AliWarning("cannot take into account the multiple scattering effects"); | |
761 | track.SetFitWithMCS(kFALSE); | |
762 | } | |
763 | } | |
764 | ||
765 | // Set fitting function | |
766 | gMinuit->SetFCN(TrackChi2); | |
208f139e | 767 | |
768 | // Set fitted parameters (!! The order is very important for the covariance matrix !!) | |
ea94c18b | 769 | trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtHit()->First()); |
de2cd600 | 770 | // could be tried with no limits for the search (min=max=0) ???? |
208f139e | 771 | // mandatory limits in non Bending to avoid NaN values of parameters |
772 | gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status); | |
773 | gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status); | |
774 | // mandatory limits in Bending to avoid NaN values of parameters | |
775 | gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status); | |
776 | gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status); | |
777 | gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status); | |
778 | ||
de2cd600 | 779 | // minimization |
208f139e | 780 | gMinuit->mnexcm("MIGRAD", arg, 0, status); |
781 | ||
782 | // Calculate the covariance matrix more accurately if required | |
783 | if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status); | |
8cde4af5 | 784 | |
de2cd600 | 785 | // get results into "invBenP", "benC", "nonBenC" ("x", "y") |
208f139e | 786 | gMinuit->GetParameter(0, x, errorParam); |
787 | trackParam->SetNonBendingCoor(x); | |
788 | gMinuit->GetParameter(1, nonBenC, errorParam); | |
de2cd600 | 789 | trackParam->SetNonBendingSlope(nonBenC); |
208f139e | 790 | gMinuit->GetParameter(2, y, errorParam); |
791 | trackParam->SetBendingCoor(y); | |
792 | gMinuit->GetParameter(3, benC, errorParam); | |
793 | trackParam->SetBendingSlope(benC); | |
794 | gMinuit->GetParameter(4, invBenP, errorParam); | |
795 | trackParam->SetInverseBendingMomentum(invBenP); | |
796 | ||
de2cd600 | 797 | // global result of the fit |
208f139e | 798 | gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus); |
ea94c18b | 799 | track.SetFitFMin(fitFMin); |
208f139e | 800 | |
801 | // Get the covariance matrix if required | |
802 | if (calcCov) { | |
803 | // Covariance matrix according to HESSE status | |
804 | // If problem then keep only the diagonal terms (variances) | |
805 | Double_t matrix[5][5]; | |
806 | gMinuit->mnemat(&matrix[0][0],5); | |
807 | if (covStatus == 3) trackParam->SetCovariances(matrix); | |
808 | else trackParam->SetVariances(matrix); | |
ea94c18b | 809 | } else trackParam->DeleteCovariances(); |
208f139e | 810 | |
de2cd600 | 811 | } |
812 | ||
813 | //__________________________________________________________________________ | |
ea94c18b | 814 | void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/) |
de2cd600 | 815 | { |
ea94c18b | 816 | /// Return the "Chi2" to be minimized with Minuit for track fitting. |
2457f726 | 817 | /// Assumes that the track hits are sorted according to increasing Z. |
818 | /// Track parameters at each TrackHit are updated accordingly. | |
ea94c18b | 819 | /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()". |
820 | /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()". | |
208f139e | 821 | |
822 | AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit(); | |
ea94c18b | 823 | AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtHit()->First(); |
208f139e | 824 | Double_t dX, dY; |
ea94c18b | 825 | chi2 = 0.; // initialize chi2 |
208f139e | 826 | |
ea94c18b | 827 | // update track parameters |
828 | trackParamAtHit->SetNonBendingCoor(param[0]); | |
829 | trackParamAtHit->SetNonBendingSlope(param[1]); | |
830 | trackParamAtHit->SetBendingCoor(param[2]); | |
831 | trackParamAtHit->SetBendingSlope(param[3]); | |
832 | trackParamAtHit->SetInverseBendingMomentum(param[4]); | |
833 | trackBeingFitted->UpdateTrackParamAtHit(); | |
208f139e | 834 | |
835 | // Take the vertex into account in the fit if required | |
836 | if (trackBeingFitted->GetFitWithVertex()) { | |
ea94c18b | 837 | AliMUONTrackParam paramAtVertex(*trackParamAtHit); |
208f139e | 838 | AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.); |
839 | AliMUONHitForRec *vertex = trackBeingFitted->GetVertex(); | |
840 | if (!vertex) { | |
841 | cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl; | |
842 | exit(-1); | |
843 | } | |
844 | dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor(); | |
845 | dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor(); | |
ea94c18b | 846 | chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2(); |
de2cd600 | 847 | } |
208f139e | 848 | |
ea94c18b | 849 | // compute chi2 w/wo multiple scattering |
850 | if (trackBeingFitted->GetFitWithMCS()) chi2 += trackBeingFitted->ComputeGlobalChi2(kTRUE); | |
851 | else chi2 += trackBeingFitted->ComputeGlobalChi2(kFALSE); | |
852 | ||
de2cd600 | 853 | } |
854 | ||
855 | //__________________________________________________________________________ | |
ea94c18b | 856 | void AliMUONTrackReconstructor::ImproveTracks() |
de2cd600 | 857 | { |
ea94c18b | 858 | /// Improve tracks by removing clusters with local chi2 highter than the defined cut |
859 | /// Recompute track parameters and covariances at the remaining clusters | |
860 | AliDebug(1,"Enter ImproveTracks"); | |
208f139e | 861 | |
ea94c18b | 862 | Double_t localChi2, worstLocalChi2; |
863 | Int_t worstChamber; | |
864 | AliMUONTrackParam *trackParamAtHit, *worstTrackParamAtHit; | |
208f139e | 865 | |
ea94c18b | 866 | // Remove double track to improve only "good" tracks |
867 | RemoveDoubleTracks(); | |
868 | ||
869 | AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First(); | |
870 | while (track) { | |
871 | ||
872 | while (!track->IsImproved()) { | |
873 | ||
874 | // Update track parameters and covariances | |
875 | track->UpdateCovTrackParamAtHit(); | |
876 | ||
877 | // Compute local chi2 of each hits | |
878 | track->ComputeLocalChi2(kTRUE); | |
879 | ||
880 | // Look for the hit to remove | |
881 | worstTrackParamAtHit = 0; | |
882 | worstLocalChi2 = 0.; | |
883 | trackParamAtHit = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First(); | |
884 | while (trackParamAtHit) { | |
885 | ||
886 | // Pick up hit with the worst chi2 | |
887 | localChi2 = trackParamAtHit->GetLocalChi2(); | |
888 | if (localChi2 > worstLocalChi2) { | |
889 | worstLocalChi2 = localChi2; | |
890 | worstTrackParamAtHit = trackParamAtHit; | |
891 | } | |
892 | ||
893 | trackParamAtHit = (AliMUONTrackParam*) track->GetTrackParamAtHit()->After(trackParamAtHit); | |
de2cd600 | 894 | } |
ea94c18b | 895 | |
896 | // Check if bad hit found | |
897 | if (!worstTrackParamAtHit) { | |
898 | track->SetImproved(kTRUE); | |
899 | break; | |
de2cd600 | 900 | } |
ea94c18b | 901 | |
902 | // check whether the worst hit is removable or not | |
903 | if (!worstTrackParamAtHit->IsRemovable()) { | |
904 | track->SetImproved(kTRUE); | |
905 | break; | |
de2cd600 | 906 | } |
ea94c18b | 907 | |
908 | // Check whether the worst chi2 is under requirement or not | |
909 | if (worstLocalChi2 < 2. * fgkSigmaToCutForImprovement * fgkSigmaToCutForImprovement) { // 2 because 2 quantities in chi2 | |
910 | track->SetImproved(kTRUE); | |
911 | break; | |
912 | } | |
913 | ||
914 | // Reset the second hit in the same station as the bad one as being NOT removable | |
915 | worstChamber = worstTrackParamAtHit->GetHitForRecPtr()->GetChamberNumber(); | |
916 | if (worstChamber%2 == 0) ((AliMUONTrackParam*)track->GetTrackParamAtHit()->After(worstTrackParamAtHit))->SetRemovable(kFALSE); | |
917 | else ((AliMUONTrackParam*)track->GetTrackParamAtHit()->Before(worstTrackParamAtHit))->SetRemovable(kFALSE); | |
918 | ||
919 | // Remove the worst hit | |
920 | track->RemoveTrackParamAtHit(worstTrackParamAtHit); | |
921 | ||
922 | // Re-fit the track: | |
923 | // Take into account the multiple scattering | |
924 | // Calculate the track parameter covariance matrix | |
925 | track->SetFitWithVertex(kFALSE); // To be sure | |
926 | Fit(*track, kTRUE, kTRUE); | |
927 | ||
928 | // Printout for debuging | |
929 | if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { | |
930 | cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(track)+1 << " improved " << endl; | |
931 | } | |
932 | ||
de2cd600 | 933 | } |
ea94c18b | 934 | |
935 | track = (AliMUONTrack*) fRecTracksPtr->After(track); | |
de2cd600 | 936 | } |
937 | ||
de2cd600 | 938 | } |
939 | ||
208f139e | 940 | //__________________________________________________________________________ |
ea94c18b | 941 | void AliMUONTrackReconstructor::Finalize() |
b8dc484b | 942 | { |
b58638a9 | 943 | /// Fill AliMUONTrack's fHitForRecAtHit array |
944 | /// Recompute track parameters and covariances at each attached cluster from those at the first one | |
b8dc484b | 945 | AliMUONTrack *track; |
de2cd600 | 946 | AliMUONTrackParam *trackParamAtHit; |
b58638a9 | 947 | |
b8dc484b | 948 | track = (AliMUONTrack*) fRecTracksPtr->First(); |
949 | while (track) { | |
ea94c18b | 950 | |
951 | // update track parameters if not already done | |
952 | if (!track->IsImproved()) track->UpdateCovTrackParamAtHit(); | |
953 | ||
de2cd600 | 954 | trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); |
955 | while (trackParamAtHit) { | |
ea94c18b | 956 | |
b58638a9 | 957 | // update array of track hit |
ea94c18b | 958 | track->AddHitForRecAtHit(trackParamAtHit->GetHitForRecPtr()); |
959 | ||
960 | trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit)); | |
de2cd600 | 961 | } |
ea94c18b | 962 | |
b8dc484b | 963 | track = (AliMUONTrack*) fRecTracksPtr->After(track); |
ea94c18b | 964 | |
de2cd600 | 965 | } |
ea94c18b | 966 | |
a9e2aefa | 967 | } |
968 | ||
276c44b7 | 969 |