AliRoot Core  ee782a0 (ee782a0)
AliMFTTrackReconstructor.cxx
Go to the documentation of this file.
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 
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 
33 ClassImp(AliMFTTrackReconstructor); // Class implementation in ROOT context
35 
36 
37 //=============================================================================================
38 
40 {
42 
43  // set the magnetic field for track extrapolations
45 
46 }
47 
48 
49 //=============================================================================================
50 
51 
53 
54 }
55 //__________________________________________________________________________
57 
58  Bool_t extrapStatus = kTRUE;
59  Double_t addChi2TrackAtCluster;
60 
61  AliMFTCATrack * currentCATrack = currentTrack->GetCATrack();
62  Int_t nCells = currentCATrack->GetNcells();
63 
64  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  AliMFTCACell * caCell = currentCATrack->GetCell(0);
68  caCell->PrintCell("");
69 
70  // Evaluate the track sign and Pt
71  currentCATrack->EvalSignedPt();
72 
73  Double_t *caHit1, *caHit2;
74  caHit1 = caCell->GetHit1();
75  caHit2 = caCell->GetHit2();
76 
77  Double_t dX = caHit1[0] - caHit2[0];
78  Double_t dY = caHit1[1] - caHit2[1];
79  Double_t dZ = caHit1[2] - caHit2[2];
80  Double_t dr = TMath::Sqrt(dX*dX+dY*dY);
81  Double_t slopeX_Z = dX / dZ;
82  Double_t slopeY_Z = dY / dZ;
83  Double_t slopeY_R = dY / dr;
84  Double_t slope2 = slopeX_Z*slopeX_Z + slopeY_Z*slopeY_Z;
85 
86  // Inverse momentum
87  Double_t inversePt = 1./currentCATrack->GetPt(); // Signed Pt estimation
88 
89  // Set track parameters at second cluster
90  AliMFTTrackParam* trackParamAtLastCluster = (AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
91 
92  trackParamAtLastCluster->SetX(caHit2[0]);
93  trackParamAtLastCluster->SetY(caHit2[1]);
94  trackParamAtLastCluster->SetZ(caHit2[2]);
95  trackParamAtLastCluster->SetSlopeX(slopeX_Z);
96  trackParamAtLastCluster->SetSlopeY(slopeY_Z);
97  trackParamAtLastCluster->SetInverseTransverseMomentum(inversePt);
98 
99  Double_t inverseMomentum;
100 
101  // Compute and set track parameters covariances at first cluster
102  TMatrixD paramCov(5,5);
103  paramCov.Zero();
106 
107 
108  paramCov(0,0) = errX2;
109  paramCov(1,1) = errY2;
110  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  paramCov(2,0) = -errX2 / dZ;
112  paramCov(0,2) = paramCov(2,0);
113 
114  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  paramCov(3,1) = -errY2 / dZ;
116  paramCov(1,3) = paramCov(3,1);
117 
118 
119 
120  //take 100% error on inverse momentum
121 
122  paramCov(4,4) = inversePt*inversePt
123  *(0.1*0.1+ (slopeX_Z*slopeX_Z *errX2 *1001. + slopeY_Z*slopeY_Z *errY2 *1001.)/dZ/dZ / slope2 / slope2);
124  paramCov(4,0) = inversePt / dZ / slope2 * errX2 * slopeX_Z;
125  paramCov(4,1) = inversePt / dZ / slope2 * errY2 * slopeY_Z;
126  paramCov(0,4) = paramCov(4,0);
127  paramCov(1,4) = paramCov(4,1);
128  paramCov(4,2) = - inversePt / slope2 * slopeX_Z * paramCov(2,2);
129  paramCov(4,3) = - inversePt / slope2 * slopeY_Z * paramCov(3,3);
130  paramCov(2,4) = paramCov(4,2);
131  paramCov(3,4) = paramCov(4,3);
132 
133  trackParamAtLastCluster->SetCovariances(paramCov);
134  AliInfo("Starting Covariance Matrix");
135  paramCov.Print();
136  // Reset the track chi2
137  trackParamAtLastCluster->SetTrackChi2(0.);
138  trackParamAtLastCluster->Print("FULL");
139 
140 
141  // Follow the track going upstream
142  AliMFTTrackParam * startingTrackParam = NULL;
143  startingTrackParam = trackParamAtLastCluster;
144  AliMFTTrackParam * trackParamAtCluster = (AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->After(startingTrackParam);
145  currentTrack->SetMCLabel(currentCATrack->GetMCindex());
146 
147  for (Int_t iCell = 0 ; iCell < nCells; iCell++) {
148  caCell = currentCATrack->GetCell(iCell);
149  caHit1 = caCell->GetHit1();
150  caHit2 = caCell->GetHit2();
151 
152  caCell->PrintCell("MC");
153 
154  // reset track parameters and their covariances
155  trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
156  trackParamAtCluster->SetZ(startingTrackParam->GetZ());
157  trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
158 
159  // add MCS effect
160  extrapStatus = AddMCSEffect(caCell, trackParamAtCluster);
161 
162  // extrapolation to the plane of the cluster attached to the current trackParamAtCluster (update the propagator)
163  if (!AliMFTTrackExtrap::ExtrapToZCov(trackParamAtCluster, caHit1[2],
164  kFALSE)) extrapStatus = kFALSE;
165 
166  // Compute new track parameters using kalman filter
167  addChi2TrackAtCluster = RunKalmanFilter(*trackParamAtCluster);
168 
169  // Update the track chi2
170  currentTrack->SetChi2(currentTrack->GetChi2() + addChi2TrackAtCluster);
171 
172  trackParamAtCluster->SetTrackChi2(currentTrack->GetChi2());
173  trackParamAtCluster->SetLocalChi2(addChi2TrackAtCluster);
174  AliInfo("Param at cluster");
175  trackParamAtCluster->Print("FULL");
176 
177 
178  // prepare next step
179  startingTrackParam = trackParamAtCluster;
180  trackParamAtCluster = (AliMFTTrackParam*) (currentTrack->GetTrackParamAtCluster()->After(startingTrackParam));
181 
182  if( (iCell == nCells-1) && trackParamAtCluster )
183  AliWarning("Reaching last cell but still some AliMFTTrackParameter objects in the array ...");
184 
185  }
186 
187  currentTrack->SetChi2(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->GetTrackChi2()/(nCells+1));
188  currentTrack->SetPhi(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->GetPhi());
189  currentTrack->SetTheta(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->GetTheta());
190  currentTrack->SetP(((AliMFTTrackParam*) currentTrack->GetTrackParamAtCluster()->Last())->P());
191  AliInfo(Form("Input Pt %f ",currentCATrack->GetPt()));
192 
193  currentTrack->SetPt(currentCATrack->GetPt());
194  currentTrack->Print();
195  return extrapStatus;
196 }
197 //__________________________________________________________________________
199 {
200  Bool_t returnStatus = kTRUE;
201  Int_t *planeIDs = currentCell->GetLayers();
202  Int_t startingPlaneID = planeIDs[1];
203  Int_t startingDiskID = startingPlaneID/2;
204  AliInfo(Form("Layer ID = %d ; Length = %d ",startingPlaneID,currentCell->GetLength()));
205 
206  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  if(startingPlaneID%2==1){
209  AliInfo(Form("Add MCS of Disk %d ",startingDiskID));
210 
212  }
213  return kTRUE;
214  }
215 
216 
217  // Applying MCS taking into account missing planes if any
218  Int_t currentPlaneID = startingPlaneID;
219  AliInfo(Form("Layer IDs %d - %d ",planeIDs[0], planeIDs[1]));
220  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  if(currentPlaneID%2==1){
224  AliInfo(Form("Add MCS of Disk %d ",currentPlaneID/2));
225  AliMFTTrackExtrap::AddMCSEffect(trackParam,-1.4,AliMFTConstants::DiskThicknessInX0(currentPlaneID/2));
226 
227  }
228  AliInfo(Form("Extrapolate to Plane ID %d ",currentPlaneID-1));
229  if (!AliMFTTrackExtrap::ExtrapToZCov(trackParam, AliMFTConstants::DefaultPlaneZ(currentPlaneID-1), kFALSE)) returnStatus = kFALSE;
230  AliInfo(Form("After extrapolation to Z = %f",trackParam->GetZ()));
231  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  currentPlaneID--;
238  }
239 
240  return returnStatus;
241 
242 }
243 //__________________________________________________________________________
244 void AliMFTTrackReconstructor::EventReconstruct(TClonesArray *fMFTTracks )
245 {
247  AliDebug(1,"");
248  AliCodeTimerAuto("",0);
249 
250  // Stop tracking if no track candidate found
251 
252  Int_t nTracks = fMFTTracks->GetEntriesFast();
253  if (nTracks == 0) return;
254 
255 
256  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
257  AliInfo(Form("Treating Track %d",iTrack));
258  TraceTrack((AliMFTTrack *) fMFTTracks->UncheckedAt(iTrack));
259  ((AliMFTTrack *) fMFTTracks->UncheckedAt(iTrack))->Print("");
260  }
261 
262 
263 
264  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 }
312 
313 //__________________________________________________________________________
315 {
318  AliInfo("Enter RunKalmanFilter");
319 
320  // Get actual track parameters (p)
321  TMatrixD param(trackParamAtCluster.GetParameters());
322 
323  // Get new cluster parameters (m)
324  TMatrixD clusterParam(5,1);
325  clusterParam.Zero();
326  clusterParam(0,0) = trackParamAtCluster.GetClusterX();
327  clusterParam(1,0) = trackParamAtCluster.GetClusterY();
328 
329  AliInfo(Form("Cluster X,Y : %f %f ",clusterParam(0,0), clusterParam(1,0)));
330 
331 
332 
333  // Compute the actual parameter weight (W)
334  TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
335  AliInfo("actual parameter covariance matrix");
336  paramWeight.Print();
337  AliInfo(Form("paramWeight.Determinant = %e ",paramWeight.Determinant()));
338 
339  if (paramWeight.Determinant() != 0) {
340  paramWeight.Invert();
341  } else {
342  AliWarning(" paramWeight Determinant = 0");
343  return 1.e6;
344  }
345  AliInfo("actual parameter weight");
346  paramWeight.Print();
347 
348  // Compute the new cluster weight (U)
349  TMatrixD clusterWeight(5,5);
350  clusterWeight.Zero();
353  clusterWeight(0,0) = 1./errX2;
354  clusterWeight(1,1) = 1./errY2;
355  AliInfo("new cluster weight");
356  clusterWeight.Print();
357  // Compute the new parameters covariance matrix ( (W+U)^-1 )
358  TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
359  AliInfo("new parameters weight");
360  newParamCov.Print();
361 
362  if (newParamCov.Determinant() != 0) {
363  newParamCov.Invert();
364  } else {
365  AliWarning(" newParamCov Determinant = 0");
366  return 1.e6;
367  }
368 
369  // Save the new parameters covariance matrix
370  AliInfo("new parameters covariance matrix");
371  newParamCov.Print();
372 
373  trackParamAtCluster.SetCovariances(newParamCov);
374 
375  // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
376  TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
377  TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
378  TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
379  AliInfo("delta parameters");
380  newParam.Print();
381  AliInfo("old parameters");
382  param.Print();
383 
384  newParam += param; // ((W+U)^-1)U(m-p) + p
385 
386  // Save the new parameters
387  trackParamAtCluster.SetParameters(newParam);
388 
389  // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
390  tmp = newParam; // p'
391  tmp -= param; // (p'-p)
392  TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
393  TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
394  tmp = newParam; // p'
395  tmp -= clusterParam; // (p'-m)
396  TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
397  addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
398 
399  return addChi2Track(0,0);
400 
401 }
402 
static const Double_t kYPixelPitch
Pixel pitch along Y.
const TMatrixD & GetParameters() const
return track parameters
const Int_t GetLength()
Definition: AliMFTCACell.h:59
void SetPhi(double val)
Definition: AliMFTTrack.h:67
void SetZ(Double_t z)
set Z coordinate (cm)
static Double_t DefaultPlaneZ(Int_t Id)
Return Plane Z position.
Class Doing MFT Track reconstruction.
Double_t GetClusterX() const
return cluster X coordinate (cm)
void dZ()
Definition: CalibAlign.C:517
Int_t * GetLayers()
Definition: AliMFTCACell.h:29
void SetLocalChi2(Double_t chi2)
set the local chi2 of the associated cluster with respect to the track
const TMatrixD & GetCovariances() const
Double_t * GetHit1()
Definition: AliMFTCACell.h:25
const Int_t GetNcells() const
Definition: AliMFTCATrack.h:26
void SetInverseTransverseMomentum(Double_t val)
set Inverse Momentum
void PrintCell(Option_t *opt)
Definition: AliMFTCACell.h:65
Bool_t AddMCSEffect(AliMFTCACell *currentCell, AliMFTTrackParam *trackParam)
void SetCovariances(const TMatrixD &covariances)
const Double_t GetPt()
Definition: AliMFTCATrack.h:55
#define AliWarning(message)
Definition: AliLog.h:541
void SetX(Double_t val)
set X coordinate (cm)
static const Double_t kXPixelPitch
Pixel pitch along X.
void SetTheta(double val)
Definition: AliMFTTrack.h:66
void Print(const char *method, TStopwatch &timer, Int_t n)
Description of an ALICE Standalone MFT track.
Definition: AliMFTTrack.h:24
void dY()
Definition: CalibAlign.C:547
void SetChi2(Double_t chi2)
set the minimum value of the function minimized by the fit
Definition: AliMFTTrack.h:37
void SetTrackChi2(Double_t chi2)
set the chi2 of the track when the associated cluster was attached
Double_t * GetHit2()
Definition: AliMFTCACell.h:26
Double_t GetChi2() const
return the minimum value of the function minimized by the fit
Definition: AliMFTTrack.h:35
void SetMCLabel(Int_t label)
set the corresponding MC track number
Definition: AliMFTTrack.h:57
virtual void Print(Option_t *opt="") const
void EvalSignedPt()
Estimate the charge sign.
#define AliInfo(message)
Definition: AliLog.h:484
#define AliCodeTimerAuto(message, counter)
Definition: AliCodeTimer.h:137
AliMFTCACell * GetCell(Int_t ic) const
Definition: AliMFTCATrack.h:28
void EventReconstruct(TClonesArray *fMFTTracks)
Double_t GetClusterY() const
return cluster Y coordinate (cm)
static void AddMCSEffect(AliMFTTrackParam *param, Double_t dZ, Double_t x0)
TObjArray * GetTrackParamAtCluster() const
Double_t GetZ() const
return Z coordinate (cm)
void SetSlopeY(Double_t val)
set Y slope
void SetP(double val)
Definition: AliMFTTrack.h:64
const Int_t GetMCindex()
Definition: AliMFTCATrack.h:50
void SetY(Double_t val)
set Y coordinate (cm)
Bool_t TraceTrack(AliMFTTrack *fMFTTracks)
#define AliDebug(logLevel, message)
Definition: AliLog.h:300
static Double_t DiskThicknessInX0(Int_t Id)
Return Disk thickness in X0.
Class holding the parameter of a MFT Standalone Track.
Double_t RunKalmanFilter(AliMFTTrackParam &trackParamAtCluster)
static void SetField()
virtual void Print(Option_t *opt="") const
static Bool_t ExtrapToZCov(AliMFTTrackParam *trackParam, Double_t zEnd, Bool_t updatePropagator=kFALSE)
void SetPt(double val)
Definition: AliMFTTrack.h:65
void SetParameters(const TMatrixD &parameters)
set track parameters
AliMFTCATrack * GetCATrack() const
return pointer to track found by Track Finder (includes clusters)
Definition: AliMFTTrack.h:50
void SetSlopeX(Double_t val)
set X slope
class TMatrixT< Double_t > TMatrixD