]>
Commit | Line | Data |
---|---|---|
8cc77df0 | 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 | ||
16 | /* $Id$ */ | |
17 | ||
18 | //////////////////////////////////// | |
19 | // | |
20 | // MUON track reconstructor using the kalman method | |
21 | // | |
22 | // This class contains as data: | |
23 | // * the parameters for the track reconstruction | |
24 | // | |
25 | // It contains as methods, among others: | |
26 | // * MakeTracks to build the tracks | |
27 | // | |
28 | //////////////////////////////////// | |
29 | ||
30 | #include <stdlib.h> | |
31 | #include <Riostream.h> | |
32 | #include <TDirectory.h> | |
33 | #include <TFile.h> | |
34 | #include <TMatrixD.h> | |
35 | ||
36 | #include "AliMUONVTrackReconstructor.h" | |
37 | #include "AliMUONTrackReconstructorK.h" | |
38 | #include "AliMUONData.h" | |
39 | #include "AliMUONConstants.h" | |
40 | #include "AliMUONHitForRec.h" | |
41 | #include "AliMUONRawCluster.h" | |
42 | #include "AliMUONSegment.h" | |
43 | #include "AliMUONTrackK.h" | |
44 | #include "AliLog.h" | |
45 | ||
46 | ClassImp(AliMUONTrackReconstructorK) // Class implementation in ROOT context | |
47 | ||
48 | //__________________________________________________________________________ | |
49 | AliMUONTrackReconstructorK::AliMUONTrackReconstructorK(AliMUONData* data, const Option_t* TrackMethod) | |
50 | : AliMUONVTrackReconstructor(data), | |
51 | fTrackMethod(2), //tracking method (2-Kalman 3-Combination-Kalman/Clustering) | |
52 | fMuons(0) | |
53 | { | |
54 | /// Constructor for class AliMUONTrackReconstructorK | |
55 | ||
56 | if (strstr(TrackMethod,"Kalman")) { | |
57 | cout << " *** Tracking with the Kalman filter *** " << endl; | |
58 | fTrackMethod = 2; | |
59 | } else if (strstr(TrackMethod,"Combi")) { | |
60 | cout << " *** Combined cluster / track finder ***" << endl; | |
61 | fTrackMethod = 3; | |
62 | } else AliFatal(Form("Tracking method %s not available",TrackMethod)); | |
63 | ||
64 | // Memory allocation for the TClonesArray of reconstructed tracks | |
65 | fRecTracksPtr = new TClonesArray("AliMUONTrackK", 10); | |
66 | } | |
67 | ||
68 | //__________________________________________________________________________ | |
69 | AliMUONTrackReconstructorK::~AliMUONTrackReconstructorK(void) | |
70 | { | |
71 | /// Destructor for class AliMUONTrackReconstructorK | |
72 | delete fRecTracksPtr; | |
73 | } | |
74 | ||
75 | //__________________________________________________________________________ | |
76 | void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters() | |
77 | { | |
78 | /// To add to the list of hits for reconstruction all the raw clusters | |
e6b25a6e | 79 | TTree *treeR; |
8cc77df0 | 80 | AliMUONHitForRec *hitForRec; |
81 | AliMUONRawCluster *clus; | |
82 | Int_t iclus, nclus, nTRentries; | |
83 | TClonesArray *rawclusters; | |
84 | AliDebug(1,"Enter AddHitsForRecFromRawClusters"); | |
85 | ||
e6b25a6e | 86 | treeR = fMUONData->TreeR(); |
87 | if (!treeR) { | |
88 | AliError("TreeR must be loaded"); | |
89 | exit(0); | |
90 | } | |
91 | ||
8cc77df0 | 92 | if (fTrackMethod != 3) { //AZ |
e6b25a6e | 93 | nTRentries = Int_t(treeR->GetEntries()); |
94 | ||
95 | if (!(fMUONData->IsRawClusterBranchesInTree())) { | |
96 | AliError(Form("RawCluster information is not avalaible, nTRentries = %d not equal to 1",nTRentries)); | |
8cc77df0 | 97 | exit(0); |
98 | } | |
e6b25a6e | 99 | |
100 | fMUONData->SetTreeAddress("RC"); //AZ | |
8cc77df0 | 101 | fMUONData->GetRawClusters(); // only one entry |
102 | } | |
103 | ||
104 | // Loop over tracking chambers | |
105 | for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) { | |
106 | // number of HitsForRec to 0 for the chamber | |
107 | fNHitsForRecPerChamber[ch] = 0; | |
108 | // index of first HitForRec for the chamber | |
109 | if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0; | |
110 | else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec; | |
111 | rawclusters =fMUONData->RawClusters(ch); | |
112 | nclus = (Int_t) (rawclusters->GetEntries()); | |
113 | // Loop over (cathode correlated) raw clusters | |
114 | for (iclus = 0; iclus < nclus; iclus++) { | |
115 | clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus); | |
116 | // new AliMUONHitForRec from raw cluster | |
117 | // and increment number of AliMUONHitForRec's (total and in chamber) | |
118 | hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus); | |
119 | fNHitsForRec++; | |
120 | (fNHitsForRecPerChamber[ch])++; | |
121 | // more information into HitForRec | |
122 | hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY()); | |
123 | hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX()); | |
124 | // original raw cluster | |
125 | hitForRec->SetChamberNumber(ch); | |
126 | hitForRec->SetHitNumber(iclus); | |
127 | // Z coordinate of the raw cluster (cm) | |
128 | hitForRec->SetZ(clus->GetZ(0)); | |
129 | ||
130 | StdoutToAliDebug(3, | |
131 | cout << "Chamber " << ch << | |
132 | " raw cluster " << iclus << " : " << endl; | |
133 | clus->Print("full"); | |
134 | cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl; | |
135 | hitForRec->Print("full"); | |
136 | ); | |
137 | } // end of cluster loop | |
138 | } // end of chamber loop | |
139 | SortHitsForRecWithIncreasingChamber(); | |
140 | ||
141 | AliDebug(1,"End of AddHitsForRecFromRawClusters"); | |
142 | if (AliLog::GetGlobalDebugLevel() > 0) { | |
143 | AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec)); | |
144 | for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) { | |
145 | AliDebug(1, Form("Chamber(0...): %d",ch)); | |
146 | AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch])); | |
147 | AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch])); | |
148 | for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch]; | |
149 | hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch]; | |
150 | hit++) { | |
151 | AliDebug(1, Form("HitForRec index(0...): %d",hit)); | |
152 | ((*fHitsForRecPtr)[hit])->Dump(); | |
153 | } | |
154 | } | |
155 | } | |
156 | ||
157 | return; | |
158 | } | |
159 | ||
160 | //__________________________________________________________________________ | |
161 | void AliMUONTrackReconstructorK::MakeSegments(void) | |
162 | { | |
163 | /// To make the list of segments in all stations, | |
164 | /// from the list of hits to be reconstructed | |
165 | AliDebug(1,"Enter MakeSegments"); | |
166 | //AZ ResetSegments(); | |
167 | // Loop over stations | |
168 | for (Int_t st = 3; st < AliMUONConstants::NTrackingCh()/2; st++) MakeSegmentsPerStation(st); | |
169 | ||
170 | StdoutToAliDebug(3, | |
171 | cout << "end of MakeSegments" << endl; | |
172 | for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) | |
173 | { | |
174 | cout << "station " << st | |
175 | << " has " << fNSegments[st] << " segments:" | |
176 | << endl; | |
177 | for (Int_t seg = 0; seg < fNSegments[st]; seg++) | |
178 | { | |
179 | ((*fSegmentsPtr[st])[seg])->Print(); | |
180 | } | |
181 | } | |
182 | ); | |
183 | return; | |
184 | } | |
185 | ||
186 | //__________________________________________________________________________ | |
187 | void AliMUONTrackReconstructorK::MakeTracks(void) | |
188 | { | |
189 | /// To make the tracks from the list of segments and points in all stations | |
190 | AliDebug(1,"Enter MakeTracks"); | |
191 | // The order may be important for the following Reset's | |
192 | //AZ ResetTracks(); | |
193 | MakeTrackCandidates(); | |
194 | if (fRecTracksPtr->GetEntriesFast() == 0) return; | |
195 | // Follow tracks in stations(1..) 3, 2 and 1 | |
196 | FollowTracks(); | |
197 | // Remove double tracks | |
198 | RemoveDoubleTracks(); | |
199 | // Propagate tracks to the vertex thru absorber | |
200 | GoToVertex(); | |
201 | // Fill AliMUONTrack data members | |
202 | FillMUONTrack(); | |
203 | } | |
204 | ||
205 | //__________________________________________________________________________ | |
206 | void AliMUONTrackReconstructorK::MakeTrackCandidates(void) | |
207 | { | |
208 | /// To make initial tracks for Kalman filter from the list of segments | |
209 | Int_t istat, iseg; | |
210 | AliMUONSegment *segment; | |
211 | AliMUONTrackK *trackK; | |
212 | ||
213 | AliDebug(1,"Enter MakeTrackCandidatesK"); | |
214 | ||
215 | AliMUONTrackK a(this, fHitsForRecPtr); | |
216 | // Loop over stations(1...) 5 and 4 | |
217 | for (istat=4; istat>=3; istat--) { | |
218 | // Loop over segments in the station | |
219 | for (iseg=0; iseg<fNSegments[istat]; iseg++) { | |
220 | // Transform segments to tracks and evaluate covariance matrix | |
221 | segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]); | |
222 | trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment); | |
223 | } // for (iseg=0;...) | |
224 | } // for (istat=4;...) | |
225 | return; | |
226 | } | |
227 | ||
228 | //__________________________________________________________________________ | |
229 | void AliMUONTrackReconstructorK::FollowTracks(void) | |
230 | { | |
231 | /// Follow tracks using Kalman filter | |
232 | Bool_t ok; | |
233 | Int_t icand, ichamBeg = 0, ichamEnd, chamBits; | |
234 | Double_t zDipole1, zDipole2; | |
235 | AliMUONTrackK *trackK; | |
236 | AliMUONHitForRec *hit; | |
237 | AliMUONRawCluster *clus; | |
238 | TClonesArray *rawclusters; | |
239 | clus = 0; rawclusters = 0; | |
240 | ||
241 | zDipole1 = fSimpleBPosition + fSimpleBLength/2; | |
242 | zDipole2 = zDipole1 - fSimpleBLength; | |
243 | ||
244 | // Print hits | |
245 | trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]); | |
246 | ||
247 | if (trackK->DebugLevel() > 0) { | |
248 | for (Int_t i1=0; i1<fNHitsForRec; i1++) { | |
249 | hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[i1]); | |
250 | printf(" Hit # %d %10.4f %10.4f %10.4f", | |
251 | hit->GetChamberNumber(), hit->GetBendingCoor(), | |
252 | hit->GetNonBendingCoor(), hit->GetZ()); | |
253 | ||
254 | // from raw clusters | |
255 | rawclusters = fMUONData->RawClusters(hit->GetChamberNumber()); | |
256 | clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit-> | |
257 | GetHitNumber()); | |
258 | printf(" %d", clus->GetTrack(1)); | |
259 | if (clus->GetTrack(2) != -1) printf(" %d \n", clus->GetTrack(2)); | |
260 | else printf("\n"); | |
261 | ||
262 | } | |
263 | } // if (trackK->DebugLevel() > 0) | |
264 | ||
265 | icand = -1; | |
266 | Int_t nSeeds; | |
267 | nSeeds = fNRecTracks; // starting number of seeds | |
268 | // Loop over track candidates | |
269 | while (icand < fNRecTracks-1) { | |
270 | icand ++; | |
271 | if (trackK->DebugLevel()>0) cout << " *** Kalman track candidate No. " << icand << endl; | |
272 | trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]); | |
273 | if (trackK->GetRecover() < 0) continue; // failed track | |
274 | ||
275 | // Discard candidate which will produce the double track | |
276 | /* | |
277 | if (icand > 0) { | |
278 | ok = CheckCandidate(icand,nSeeds); | |
279 | if (!ok) { | |
280 | trackK->SetRecover(-1); // mark candidate to be removed | |
281 | continue; | |
282 | } | |
283 | } | |
284 | */ | |
285 | ||
286 | ok = kTRUE; | |
287 | if (trackK->GetRecover() == 0) | |
288 | hit = (AliMUONHitForRec*) trackK->GetTrackHits()->Last(); // last hit | |
289 | else | |
290 | hit = trackK->GetHitLastOk(); // hit where track stopped | |
291 | ||
292 | if (hit) ichamBeg = hit->GetChamberNumber(); | |
293 | ichamEnd = 0; | |
294 | // Check propagation direction | |
295 | if (!hit) { ichamBeg = ichamEnd; AliFatal(" ??? "); } | |
296 | else if (trackK->GetTrackDir() < 0) { | |
297 | ichamEnd = 9; // forward propagation | |
298 | ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2); | |
299 | if (ok) { | |
300 | ichamBeg = ichamEnd; | |
301 | ichamEnd = 6; // backward propagation | |
302 | // Change weight matrix and zero fChi2 for backpropagation | |
303 | trackK->StartBack(); | |
304 | trackK->SetTrackDir(1); | |
305 | ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2); | |
306 | ichamBeg = ichamEnd; | |
307 | ichamEnd = 0; | |
308 | } | |
309 | } else { | |
310 | if (trackK->GetBPFlag()) { | |
311 | // backpropagation | |
312 | ichamEnd = 6; // backward propagation | |
313 | // Change weight matrix and zero fChi2 for backpropagation | |
314 | trackK->StartBack(); | |
315 | ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kTRUE,zDipole1,zDipole2); | |
316 | ichamBeg = ichamEnd; | |
317 | ichamEnd = 0; | |
318 | } | |
319 | } | |
320 | ||
321 | if (ok) { | |
322 | trackK->SetTrackDir(1); | |
323 | trackK->SetBPFlag(kFALSE); | |
324 | ok = trackK->KalmanFilter(ichamBeg,ichamEnd,kFALSE,zDipole1,zDipole2); | |
325 | } | |
326 | if (!ok) { trackK->SetRecover(-1); continue; } // mark candidate to be removed | |
327 | ||
328 | // Apply smoother | |
329 | if (trackK->GetRecover() >= 0) { | |
330 | ok = trackK->Smooth(); | |
331 | if (!ok) trackK->SetRecover(-1); // mark candidate to be removed | |
332 | } | |
333 | ||
334 | // Majority 3 of 4 in first 2 stations | |
335 | if (!ok) continue; | |
336 | chamBits = 0; | |
337 | Double_t chi2max = 0; | |
338 | for (Int_t i=0; i<trackK->GetNTrackHits(); i++) { | |
339 | hit = (AliMUONHitForRec*) (*trackK->GetTrackHits())[i]; | |
340 | chamBits |= BIT(hit->GetChamberNumber()); | |
341 | if (trackK->GetChi2PerPoint(i) > chi2max) chi2max = trackK->GetChi2PerPoint(i); | |
342 | } | |
343 | if (!((chamBits&3)==3 || (chamBits>>2&3)==3) && chi2max > 25) { | |
344 | //trackK->Recover(); | |
345 | trackK->SetRecover(-1); //mark candidate to be removed | |
346 | continue; | |
347 | } | |
348 | if (ok) trackK->SetTrackQuality(0); // compute "track quality" | |
349 | } // while | |
350 | ||
351 | for (Int_t i=0; i<fNRecTracks; i++) { | |
352 | trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[i]); | |
353 | if (trackK->GetRecover() < 0) fRecTracksPtr->RemoveAt(i); | |
354 | } | |
355 | ||
356 | // Compress TClonesArray | |
357 | fRecTracksPtr->Compress(); | |
358 | fNRecTracks = fRecTracksPtr->GetEntriesFast(); | |
359 | return; | |
360 | } | |
361 | ||
362 | //__________________________________________________________________________ | |
363 | Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) const | |
364 | { | |
365 | /// Discards track candidate if it will produce the double track (having | |
366 | /// the same seed segment hits as hits of a good track found before) | |
367 | AliMUONTrackK *track1, *track2; | |
368 | AliMUONHitForRec *hit1, *hit2, *hit; | |
369 | ||
370 | track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]); | |
371 | hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit | |
372 | hit2 = (AliMUONHitForRec*) (*track1->GetTrackHits())[1]; // 2'nd hit | |
373 | ||
374 | for (Int_t i=0; i<icand; i++) { | |
375 | track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]); | |
376 | //if (track2->GetRecover() < 0) continue; | |
377 | if (track2->GetRecover() < 0 && icand >= nSeeds) continue; | |
378 | ||
379 | if (track1->GetStartSegment() == track2->GetStartSegment()) { | |
380 | return kFALSE; | |
381 | } else { | |
382 | Int_t nSame = 0; | |
383 | for (Int_t j=0; j<track2->GetNTrackHits(); j++) { | |
384 | hit = (AliMUONHitForRec*) (*track2->GetTrackHits())[j]; | |
385 | if (hit == hit1 || hit == hit2) { | |
386 | nSame++; | |
387 | if (nSame == 2) return kFALSE; | |
388 | } | |
389 | } // for (Int_t j=0; | |
390 | } | |
391 | } // for (Int_t i=0; | |
392 | return kTRUE; | |
393 | } | |
394 | ||
395 | //__________________________________________________________________________ | |
396 | void AliMUONTrackReconstructorK::RemoveDoubleTracks(void) | |
397 | { | |
398 | /// Removes double tracks (sharing more than half of their hits). | |
399 | /// Keeps the track with higher quality | |
400 | AliMUONTrackK *track1, *track2, *trackToKill; | |
401 | ||
402 | // Sort tracks according to their quality | |
403 | fRecTracksPtr->Sort(); | |
404 | ||
405 | // Loop over first track of the pair | |
406 | track1 = (AliMUONTrackK*) fRecTracksPtr->First(); | |
407 | Int_t debug = track1->DebugLevel(); | |
408 | while (track1) { | |
409 | // Loop over second track of the pair | |
410 | track2 = (AliMUONTrackK*) fRecTracksPtr->After(track1); | |
411 | while (track2) { | |
412 | // Check whether or not to keep track2 | |
413 | if (!track2->KeepTrack(track1)) { | |
414 | if (debug >= 0) cout << " Killed track: " << 1/(*track2->GetTrackParameters())(4,0) << | |
415 | " " << track2->GetTrackQuality() << endl; | |
416 | trackToKill = track2; | |
417 | track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2); | |
418 | trackToKill->Kill(); | |
419 | fRecTracksPtr->Compress(); | |
420 | } else track2 = (AliMUONTrackK*) fRecTracksPtr->After(track2); | |
421 | } // track2 | |
422 | track1 = (AliMUONTrackK*) fRecTracksPtr->After(track1); | |
423 | } // track1 | |
424 | ||
425 | fNRecTracks = fRecTracksPtr->GetEntriesFast(); | |
426 | if (debug >= 0) cout << " Number of Kalman tracks: " << fNRecTracks << endl; | |
427 | } | |
428 | ||
429 | //__________________________________________________________________________ | |
430 | void AliMUONTrackReconstructorK::GoToVertex(void) | |
431 | { | |
432 | /// Propagates track to the vertex thru absorber | |
433 | /// (using Branson correction for now) | |
434 | Double_t zVertex; | |
435 | zVertex = 0; | |
436 | for (Int_t i=0; i<fNRecTracks; i++) { | |
437 | //((AliMUONTrackK*)(*fRecTracksPtr)[i])->Branson(); | |
438 | ((AliMUONTrackK*)(*fRecTracksPtr)[i])->SetTrackQuality(1); // compute Chi2 | |
439 | //((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToZ(zVertex); // w/out absorber | |
440 | ((AliMUONTrackK*)(*fRecTracksPtr)[i])->GoToVertex(1); // with absorber | |
441 | } | |
442 | } | |
443 | ||
444 | //__________________________________________________________________________ | |
445 | void AliMUONTrackReconstructorK::FillMUONTrack() | |
446 | { | |
447 | // Set track parameters at hits for Kalman track. Fill fTrackParamAtHit of AliMUONTrack's | |
448 | AliMUONTrackK *track; | |
449 | track = (AliMUONTrackK*) fRecTracksPtr->First(); | |
450 | while (track) { | |
451 | track->FillMUONTrack(); | |
452 | track = (AliMUONTrackK*) fRecTracksPtr->After(track); | |
453 | } | |
454 | return; | |
455 | } | |
456 | ||
457 | //__________________________________________________________________________ | |
458 | void AliMUONTrackReconstructorK::EventDump(void) | |
459 | { | |
460 | /// Dump reconstructed event (track parameters at vertex and at first hit), | |
461 | /// and the particle parameters | |
462 | AliDebug(1,"****** enter EventDump ******"); | |
463 | AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); | |
464 | ||
465 | AliWarning("AliMUONTrackReconstructorK::EventDump: not implemented"); | |
466 | ||
467 | // informations about generated particles NO !!!!!!!! | |
468 | ||
469 | // for (Int_t iPart = 0; iPart < np; iPart++) { | |
470 | // p = gAlice->Particle(iPart); | |
471 | // printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n", | |
472 | // iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode()); | |
473 | // } | |
474 | return; | |
475 | } | |
476 | ||
477 |