Line data Source code
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 : #include "TClonesArray.h"
19 : #include "TMath.h"
20 :
21 : #include "AliCodeTimer.h"
22 : #include "AliLog.h"
23 :
24 : #include "AliMFTTrackReconstructor.h"
25 : #include "AliMFTTrackParam.h"
26 : #include "AliMFTTrack.h"
27 : #include "AliMFTTrackExtrap.h"
28 : #include "AliMFTCATrack.h"
29 : #include "AliMFTCACell.h"
30 : #include "AliMFTConstants.h"
31 :
32 : /// \cond CLASSIMP
33 12 : ClassImp(AliMFTTrackReconstructor); // Class implementation in ROOT context
34 : /// \endcond
35 :
36 :
37 : //=============================================================================================
38 :
39 0 : AliMFTTrackReconstructor::AliMFTTrackReconstructor():TObject()
40 0 : {
41 : /// Default constructor
42 :
43 : // set the magnetic field for track extrapolations
44 0 : AliMFTTrackExtrap::SetField();
45 :
46 0 : }
47 :
48 :
49 : //=============================================================================================
50 :
51 :
52 0 : AliMFTTrackReconstructor::~AliMFTTrackReconstructor() {
53 :
54 0 : }
55 : //__________________________________________________________________________
56 : Bool_t AliMFTTrackReconstructor::TraceTrack(AliMFTTrack *currentTrack ){
57 :
58 : Bool_t extrapStatus = kTRUE;
59 : Double_t addChi2TrackAtCluster;
60 :
61 0 : AliMFTCATrack * currentCATrack = currentTrack->GetCATrack();
62 0 : Int_t nCells = currentCATrack->GetNcells();
63 :
64 0 : if (nCells < 2) return kFALSE; // Skip tracks wih less than 2 cells
65 :
66 : // Initiate the seed starting from the last cells (the more downstream one)
67 0 : AliMFTCACell * caCell = currentCATrack->GetCell(0);
68 0 : caCell->PrintCell("");
69 :
70 : // Evaluate the track sign and Pt
71 0 : currentCATrack->EvalSignedPt();
72 :
73 : Double_t *caHit1, *caHit2;
74 0 : caHit1 = caCell->GetHit1();
75 0 : caHit2 = caCell->GetHit2();
76 :
77 0 : Double_t dX = caHit1[0] - caHit2[0];
78 0 : Double_t dY = caHit1[1] - caHit2[1];
79 0 : Double_t dZ = caHit1[2] - caHit2[2];
80 0 : Double_t dr = TMath::Sqrt(dX*dX+dY*dY);
81 0 : Double_t slopeX_Z = dX / dZ;
82 0 : Double_t slopeY_Z = dY / dZ;
83 : Double_t slopeY_R = dY / dr;
84 0 : Double_t slope2 = slopeX_Z*slopeX_Z + slopeY_Z*slopeY_Z;
85 :
86 : // Inverse momentum
87 0 : Double_t inversePt = 1./currentCATrack->GetPt(); // Signed Pt estimation
88 :
89 : // Set track parameters at second cluster
90 0 : AliMFTTrackParam* trackParamAtLastCluster = (AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
91 :
92 0 : trackParamAtLastCluster->SetX(caHit2[0]);
93 0 : trackParamAtLastCluster->SetY(caHit2[1]);
94 0 : trackParamAtLastCluster->SetZ(caHit2[2]);
95 0 : trackParamAtLastCluster->SetSlopeX(slopeX_Z);
96 0 : trackParamAtLastCluster->SetSlopeY(slopeY_Z);
97 0 : trackParamAtLastCluster->SetInverseTransverseMomentum(inversePt);
98 :
99 : Double_t inverseMomentum;
100 :
101 : // Compute and set track parameters covariances at first cluster
102 0 : TMatrixD paramCov(5,5);
103 0 : paramCov.Zero();
104 0 : Double_t errX2 = AliMFTConstants::kXPixelPitch * AliMFTConstants::kXPixelPitch / 12.;
105 0 : Double_t errY2 = AliMFTConstants::kYPixelPitch * AliMFTConstants::kYPixelPitch / 12.;
106 :
107 :
108 0 : paramCov(0,0) = errX2;
109 0 : paramCov(1,1) = errY2;
110 0 : paramCov(2,2) = ( 1001. * errX2 )/ dZ / dZ;// Weight 1 for the last cluster and 1000 for the first one. the first cluster error will be taken into account at the first step of the Kalman filter
111 0 : paramCov(2,0) = -errX2 / dZ;
112 0 : paramCov(0,2) = paramCov(2,0);
113 :
114 0 : paramCov(3,3) = ( 1001. * errY2 )/ dZ / dZ;// Weight 1 for the last cluster and 1000 for the first one. the first cluster error will be taken into account at the first step of the Kalman filter
115 0 : paramCov(3,1) = -errY2 / dZ;
116 0 : paramCov(1,3) = paramCov(3,1);
117 :
118 :
119 :
120 : //take 100% error on inverse momentum
121 :
122 0 : paramCov(4,4) = inversePt*inversePt
123 0 : *(0.1*0.1+ (slopeX_Z*slopeX_Z *errX2 *1001. + slopeY_Z*slopeY_Z *errY2 *1001.)/dZ/dZ / slope2 / slope2);
124 0 : paramCov(4,0) = inversePt / dZ / slope2 * errX2 * slopeX_Z;
125 0 : paramCov(4,1) = inversePt / dZ / slope2 * errY2 * slopeY_Z;
126 0 : paramCov(0,4) = paramCov(4,0);
127 0 : paramCov(1,4) = paramCov(4,1);
128 0 : paramCov(4,2) = - inversePt / slope2 * slopeX_Z * paramCov(2,2);
129 0 : paramCov(4,3) = - inversePt / slope2 * slopeY_Z * paramCov(3,3);
130 0 : paramCov(2,4) = paramCov(4,2);
131 0 : paramCov(3,4) = paramCov(4,3);
132 :
133 0 : trackParamAtLastCluster->SetCovariances(paramCov);
134 0 : AliInfo("Starting Covariance Matrix");
135 0 : paramCov.Print();
136 : // Reset the track chi2
137 0 : trackParamAtLastCluster->SetTrackChi2(0.);
138 0 : trackParamAtLastCluster->Print("FULL");
139 :
140 :
141 : // Follow the track going upstream
142 : AliMFTTrackParam * startingTrackParam = NULL;
143 : startingTrackParam = trackParamAtLastCluster;
144 0 : AliMFTTrackParam * trackParamAtCluster = (AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->After(startingTrackParam);
145 0 : currentTrack->SetMCLabel(currentCATrack->GetMCindex());
146 :
147 0 : for (Int_t iCell = 0 ; iCell < nCells; iCell++) {
148 0 : caCell = currentCATrack->GetCell(iCell);
149 0 : caHit1 = caCell->GetHit1();
150 0 : caHit2 = caCell->GetHit2();
151 :
152 0 : caCell->PrintCell("MC");
153 :
154 : // reset track parameters and their covariances
155 0 : trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
156 0 : trackParamAtCluster->SetZ(startingTrackParam->GetZ());
157 0 : trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
158 :
159 : // add MCS effect
160 0 : extrapStatus = AddMCSEffect(caCell, trackParamAtCluster);
161 :
162 : // extrapolation to the plane of the cluster attached to the current trackParamAtCluster (update the propagator)
163 0 : if (!AliMFTTrackExtrap::ExtrapToZCov(trackParamAtCluster, caHit1[2],
164 0 : kFALSE)) extrapStatus = kFALSE;
165 :
166 : // Compute new track parameters using kalman filter
167 0 : addChi2TrackAtCluster = RunKalmanFilter(*trackParamAtCluster);
168 :
169 : // Update the track chi2
170 0 : currentTrack->SetChi2(currentTrack->GetChi2() + addChi2TrackAtCluster);
171 :
172 0 : trackParamAtCluster->SetTrackChi2(currentTrack->GetChi2());
173 0 : trackParamAtCluster->SetLocalChi2(addChi2TrackAtCluster);
174 0 : AliInfo("Param at cluster");
175 0 : trackParamAtCluster->Print("FULL");
176 :
177 :
178 : // prepare next step
179 : startingTrackParam = trackParamAtCluster;
180 0 : trackParamAtCluster = (AliMFTTrackParam*) (currentTrack->GetTrackParamAtCluster()->After(startingTrackParam));
181 :
182 0 : if( (iCell == nCells-1) && trackParamAtCluster )
183 0 : AliWarning("Reaching last cell but still some AliMFTTrackParameter objects in the array ...");
184 :
185 : }
186 :
187 0 : currentTrack->SetChi2(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->GetTrackChi2()/(nCells+1));
188 0 : currentTrack->SetPhi(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->GetPhi());
189 0 : currentTrack->SetTheta(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->GetTheta());
190 0 : currentTrack->SetP(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->P());
191 0 : AliInfo(Form("Input Pt %f ",currentCATrack->GetPt()));
192 :
193 0 : currentTrack->SetPt(currentCATrack->GetPt());
194 0 : currentTrack->Print();
195 0 : return extrapStatus;
196 0 : }
197 : //__________________________________________________________________________
198 : Bool_t AliMFTTrackReconstructor::AddMCSEffect(AliMFTCACell *currentCell, AliMFTTrackParam *trackParam)
199 : {
200 : Bool_t returnStatus = kTRUE;
201 0 : Int_t *planeIDs = currentCell->GetLayers();
202 0 : Int_t startingPlaneID = planeIDs[1];
203 0 : Int_t startingDiskID = startingPlaneID/2;
204 0 : AliInfo(Form("Layer ID = %d ; Length = %d ",startingPlaneID,currentCell->GetLength()));
205 :
206 0 : if(currentCell->GetLength()==1){
207 : // No MCS effect between front face and back face of 2 different disks (it is the air between two disks, neglect it for now)
208 0 : if(startingPlaneID%2==1){
209 0 : AliInfo(Form("Add MCS of Disk %d ",startingDiskID));
210 :
211 0 : AliMFTTrackExtrap::AddMCSEffect(trackParam,-1.4,AliMFTConstants::DiskThicknessInX0(startingDiskID));
212 0 : }
213 0 : return kTRUE;
214 : }
215 :
216 :
217 : // Applying MCS taking into account missing planes if any
218 : Int_t currentPlaneID = startingPlaneID;
219 0 : AliInfo(Form("Layer IDs %d - %d ",planeIDs[0], planeIDs[1]));
220 0 : while (currentPlaneID != planeIDs[0]) {
221 : // extrapolation to the missing chamber (update the propagator)
222 : // add MCS effect if second hit is in the back plane face of a disk
223 0 : if(currentPlaneID%2==1){
224 0 : AliInfo(Form("Add MCS of Disk %d ",currentPlaneID/2));
225 0 : AliMFTTrackExtrap::AddMCSEffect(trackParam,-1.4,AliMFTConstants::DiskThicknessInX0(currentPlaneID/2));
226 :
227 0 : }
228 0 : AliInfo(Form("Extrapolate to Plane ID %d ",currentPlaneID-1));
229 0 : if (!AliMFTTrackExtrap::ExtrapToZCov(trackParam, AliMFTConstants::DefaultPlaneZ(currentPlaneID-1), kFALSE)) returnStatus = kFALSE;
230 0 : AliInfo(Form("After extrapolation to Z = %f",trackParam->GetZ()));
231 0 : trackParam->Print("FULL");
232 :
233 : // // add MCS effect
234 : // AliInfo(Form("Add MCS of Disk %d ",currentPlaneID/2));
235 : // AliMFTTrackExtrap::AddMCSEffect(trackParam,-1.4,AliMFTConstants::DiskThicknessInX0(startingDiskID));
236 :
237 0 : currentPlaneID--;
238 : }
239 :
240 0 : return returnStatus;
241 :
242 0 : }
243 : //__________________________________________________________________________
244 : void AliMFTTrackReconstructor::EventReconstruct(TClonesArray *fMFTTracks )
245 : {
246 : /// To reconstruct one event
247 0 : AliDebug(1,"");
248 0 : AliCodeTimerAuto("",0);
249 :
250 : // Stop tracking if no track candidate found
251 :
252 0 : Int_t nTracks = fMFTTracks->GetEntriesFast();
253 0 : if (nTracks == 0) return;
254 :
255 :
256 0 : for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
257 0 : AliInfo(Form("Treating Track %d",iTrack));
258 0 : TraceTrack((AliMFTTrack *) fMFTTracks->UncheckedAt(iTrack));
259 0 : ((AliMFTTrack *) fMFTTracks->UncheckedAt(iTrack))->Print("");
260 : }
261 :
262 :
263 :
264 0 : AliInfo("I am HERE ....");
265 : // Reset array of tracks
266 : // ResetTracks();
267 : //
268 : // // Look for candidates from clusters in stations(1..) 4 and 5 (abort in case of failure)
269 : // if (!MakeTrackCandidates(clusterStore)) return;
270 : //
271 : // // Look for extra candidates from clusters in stations(1..) 4 and 5 (abort in case of failure)
272 : // if (GetRecoParam()->MakeMoreTrackCandidates()) {
273 : // if (!MakeMoreTrackCandidates(clusterStore)) return;
274 : // }
275 : //
276 : //
277 : // // Follow tracks in stations(1..) 3, 2 and 1 (abort in case of failure)
278 : // if (!FollowTracks(clusterStore)) return;
279 : //
280 : // // Complement the reconstructed tracks
281 : // if (GetRecoParam()->ComplementTracks()) {
282 : // if (ComplementTracks(clusterStore)) RemoveIdenticalTracks();
283 : // }
284 : //
285 : // // Improve the reconstructed tracks
286 : // if (GetRecoParam()->ImproveTracks()) ImproveTracks();
287 : //
288 : // // Remove connected tracks
289 : // RemoveConnectedTracks(3, 4, kFALSE);
290 : // RemoveConnectedTracks(2, 2, kFALSE);
291 : // if (GetRecoParam()->RemoveConnectedTracksInSt12()) RemoveConnectedTracks(0, 1, kFALSE);
292 : //
293 : // // Fill AliMUONTrack data members
294 : // Finalize();
295 : // if (!GetRecoParam()->RemoveConnectedTracksInSt12()) TagConnectedTracks(0, 1, kTRUE);
296 : //
297 : // // Make sure there is no bad track left
298 : // RemoveBadTracks();
299 : //
300 : // // Refit the reconstructed tracks with a different resolution for mono-cathod clusters
301 : // if (GetRecoParam()->DiscardMonoCathodClusters()) DiscardMonoCathodClusters();
302 : //
303 : // // Add tracks to MUON data container
304 : // for (Int_t i=0; i<fNRecTracks; ++i)
305 : // {
306 : // AliMUONTrack * track = (AliMUONTrack*) fRecTracksPtr->At(i);
307 : // track->SetUniqueID(i+1);
308 : // trackStore.Add(*track);
309 : // }
310 :
311 0 : }
312 :
313 : //__________________________________________________________________________
314 : Double_t AliMFTTrackReconstructor::RunKalmanFilter(AliMFTTrackParam &trackParamAtCluster)
315 : {
316 : /// Compute new track parameters and their covariances including new cluster using kalman filter
317 : /// return the additional track chi2
318 0 : AliInfo("Enter RunKalmanFilter");
319 :
320 : // Get actual track parameters (p)
321 0 : TMatrixD param(trackParamAtCluster.GetParameters());
322 :
323 : // Get new cluster parameters (m)
324 0 : TMatrixD clusterParam(5,1);
325 0 : clusterParam.Zero();
326 0 : clusterParam(0,0) = trackParamAtCluster.GetClusterX();
327 0 : clusterParam(1,0) = trackParamAtCluster.GetClusterY();
328 :
329 0 : AliInfo(Form("Cluster X,Y : %f %f ",clusterParam(0,0), clusterParam(1,0)));
330 :
331 :
332 :
333 : // Compute the actual parameter weight (W)
334 0 : TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
335 0 : AliInfo("actual parameter covariance matrix");
336 0 : paramWeight.Print();
337 0 : AliInfo(Form("paramWeight.Determinant = %e ",paramWeight.Determinant()));
338 :
339 0 : if (paramWeight.Determinant() != 0) {
340 0 : paramWeight.Invert();
341 : } else {
342 0 : AliWarning(" paramWeight Determinant = 0");
343 0 : return 1.e6;
344 : }
345 0 : AliInfo("actual parameter weight");
346 0 : paramWeight.Print();
347 :
348 : // Compute the new cluster weight (U)
349 0 : TMatrixD clusterWeight(5,5);
350 0 : clusterWeight.Zero();
351 0 : Double_t errX2 = AliMFTConstants::kXPixelPitch * AliMFTConstants::kXPixelPitch / 12.;
352 0 : Double_t errY2 = AliMFTConstants::kYPixelPitch * AliMFTConstants::kYPixelPitch / 12.;
353 0 : clusterWeight(0,0) = 1./errX2;
354 0 : clusterWeight(1,1) = 1./errY2;
355 0 : AliInfo("new cluster weight");
356 0 : clusterWeight.Print();
357 : // Compute the new parameters covariance matrix ( (W+U)^-1 )
358 0 : TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
359 0 : AliInfo("new parameters weight");
360 0 : newParamCov.Print();
361 :
362 0 : if (newParamCov.Determinant() != 0) {
363 0 : newParamCov.Invert();
364 : } else {
365 0 : AliWarning(" newParamCov Determinant = 0");
366 0 : return 1.e6;
367 : }
368 :
369 : // Save the new parameters covariance matrix
370 0 : AliInfo("new parameters covariance matrix");
371 0 : newParamCov.Print();
372 :
373 0 : trackParamAtCluster.SetCovariances(newParamCov);
374 :
375 : // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
376 0 : TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
377 0 : TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
378 0 : TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
379 0 : AliInfo("delta parameters");
380 0 : newParam.Print();
381 0 : AliInfo("old parameters");
382 0 : param.Print();
383 :
384 0 : newParam += param; // ((W+U)^-1)U(m-p) + p
385 :
386 : // Save the new parameters
387 0 : trackParamAtCluster.SetParameters(newParam);
388 :
389 : // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
390 0 : tmp = newParam; // p'
391 0 : tmp -= param; // (p'-p)
392 0 : TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
393 0 : TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
394 0 : tmp = newParam; // p'
395 0 : tmp -= clusterParam; // (p'-m)
396 0 : TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
397 0 : addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
398 :
399 0 : return addChi2Track(0,0);
400 :
401 0 : }
402 :
|
