MFT Digitization algorithm modified
[u/mrichter/AliRoot.git] / MFT / AliMFTTrackerMU.cxx
CommitLineData
4cc511f4
AU
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//====================================================================================================================================================
17//
18// MFT tracker
19//
20// Class for the creation of the "global muon tracks" built from the tracks reconstructed in the
21// muon spectrometer and the clusters of the Muon Forward Tracker
22//
23// Contact author: antonio.uras@cern.ch
24//
25//====================================================================================================================================================
26
27#include "TTree.h"
28#include "AliLog.h"
29#include "AliGeomManager.h"
30#include "AliESDEvent.h"
31#include "AliESDMuonTrack.h"
33e3eaf3 32// #include "AliESDMuonGlobalTrack.h"
4cc511f4
AU
33#include "AliMFTTrackerMU.h"
34#include "TMath.h"
35#include "AliRun.h"
36#include "AliMFT.h"
37#include "AliMUONTrackExtrap.h"
38#include "AliMUONTrack.h"
39#include "AliMUONESDInterface.h"
40#include "AliMuonForwardTrack.h"
41
42ClassImp(AliMFTTrackerMU)
43
44const Double_t AliMFTTrackerMU::fRadLengthSi = AliMFTConstants::fRadLengthSi;
45
46//====================================================================================================================================================
47
48AliMFTTrackerMU::AliMFTTrackerMU() :
49 AliTracker(),
50 fESD(0),
51 fMFT(0),
52 fSegmentation(0),
53 fNPlanesMFT(0),
54 fNPlanesMFTAnalyzed(0),
55 fSigmaClusterCut(0),
56 fScaleSigmaClusterCut(1.),
57 fNMaxMissingMFTClusters(0),
58 fGlobalTrackingDiverged(kFALSE),
59 fCandidateTracks(0),
60 fMUONTrack(0),
61 fCurrentTrack(0),
62 fFinalBestCandidate(0),
63 fVertexErrorX(0.015),
64 fVertexErrorY(0.015),
65 fVertexErrorZ(0.010),
66 fBransonCorrection(kFALSE)
67{
68
69 //--------------------------------------------------------------------
70 // This is the AliMFTTrackerMU constructor
71 //--------------------------------------------------------------------
72
73 fMFT = (AliMFT*) gAlice->GetDetector("MFT");
74 fSegmentation = fMFT->GetSegmentation();
75 SetNPlanesMFT(fSegmentation->GetNPlanes());
76 AliMUONTrackExtrap::SetField(); // set the magnetic field for track extrapolations
77
78 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
79 fMFTClusterArray[iPlane] = 0;
80 fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
81 fMFTClusterArrayBack[iPlane] = new TClonesArray("AliMFTCluster");
82 fIsPlaneMandatory[iPlane] = kFALSE;
83 }
84
85}
86
87//====================================================================================================================================================
88
89AliMFTTrackerMU::~AliMFTTrackerMU() {
90
91 // destructor
92
93 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
94 delete fMFTClusterArray[iPlane];
95 delete fMFTClusterArrayFront[iPlane];
96 delete fMFTClusterArrayBack[iPlane];
97 }
98
99}
100
101//====================================================================================================================================================
102
103Int_t AliMFTTrackerMU::LoadClusters(TTree *cTree) {
104
105 //--------------------------------------------------------------------
106 // This function loads the MFT clusters
107 //--------------------------------------------------------------------
108
109 if (!cTree->GetEvent()) return kFALSE;
110 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
111 AliDebug(1, Form("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries()));
112 fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
113 }
114 SeparateFrontBackClusters();
115
116 return 0;
117
118}
119
120//====================================================================================================================================================
121
122void AliMFTTrackerMU::UnloadClusters() {
123
124 //--------------------------------------------------------------------
125 // This function unloads MFT clusters
126 //--------------------------------------------------------------------
127
128 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
129 fMFTClusterArray[iPlane] -> Clear("C");
130 fMFTClusterArrayFront[iPlane] -> Clear("C");
131 fMFTClusterArrayBack[iPlane] -> Clear("C");
132 }
133
134}
135
136//====================================================================================================================================================
137
138Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
139
140 //--------------------------------------------------------------------
141 // This functions reconstructs the Muon Forward Tracks
142 // The clusters must be already loaded !
143 //--------------------------------------------------------------------
144
145 TTree *outputTreeMuonGlobalTracks = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
146 TClonesArray *muonForwardTracks = new TClonesArray("AliMuonForwardTrack");
147 outputTreeMuonGlobalTracks -> Branch("tracks", &muonForwardTracks);
148
149 //--------------------------------------------------------------------
150
151 fESD = event;
152
153 //----------- Read ESD MUON tracks -------------------
154
155 Int_t nTracksMUON = event->GetNumberOfMuonTracks();
156
157 AliInfo(Form("Number of ESD MUON tracks: %d\n", nTracksMUON));
158
159 Int_t iTrack=0;
160 while (iTrack<nTracksMUON) {
161
162 fNPlanesMFTAnalyzed = 0;
163
164 AliDebug(1, "**************************************************************************************\n");
165 AliDebug(1, Form("*************************** MUON TRACK %3d ***************************************\n", iTrack));
166 AliDebug(1, "**************************************************************************************\n");
167
168 fCandidateTracks -> Delete();
169
170 fNPlanesMFTAnalyzed = 0;
171
172 const AliESDMuonTrack *esdTrack = event->GetMuonTrack(iTrack);
173 fMUONTrack = NULL;
174 AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kFALSE);
175
176 // the track we are going to build, starting from fMUONTrack and adding the MFT clusters
177 AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
178 track -> SetMUONTrack(new AliMUONTrack(*fMUONTrack));
179 track -> SetMCLabel(fMUONTrack->GetMCLabel());
180 track -> SetMatchTrigger(fMUONTrack->GetMatchTrigger());
181
182 //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
183
184 for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) { /* *** do not reverse the order of this cycle!!!
185 *** this reflects the fact that the extrapolation is performed
186 *** starting from the last MFT plane back to the origin */
187
188 // --------- updating the array of candidates according to the clusters available in the i-th plane ---------
189
190 fNPlanesMFTAnalyzed++;
191
192 Int_t nCandidates = fCandidateTracks->GetEntriesFast();
193 for (Int_t iCandidate=0; iCandidate<nCandidates; iCandidate++) {
194 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iCandidate);
195 // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
196 // (several new tracks can be created for one old track)
197 if (FindClusterInPlane(iPlane) == kDiverged) {
198 fGlobalTrackingDiverged = kTRUE;
199 break;
200 }
201
202 if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
203 fCandidateTracks->Remove(fCurrentTrack); // the old track is removed after the check;
204 }
205 }
206 if (fGlobalTrackingDiverged) {
207 if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
208 continue;
209 }
210
211 fCandidateTracks->Compress();
212
213 }
214
215 // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
216
217 fGlobalTrackingDiverged = kFALSE;
218 fScaleSigmaClusterCut = 1.0;
219
220 AliDebug(1, "Finished cycle over planes");
221
222 iTrack++;
223
224 // If we have several final tracks, we must find the best candidate:
225
226 Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
227 AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
228
229 Int_t nGoodClustersBestCandidate = 0;
230 Int_t idBestCandidate = 0;
231 Double_t bestChi2 = -1.; // variable defining the best candidate
232
233 for (Int_t iFinalCandidate=0; iFinalCandidate<nFinalTracks; iFinalCandidate++) {
234
235 AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iFinalCandidate);
236 Int_t nMFTClusters = finalTrack->GetNMFTClusters();
237
238 Double_t chi2 = 0;
239 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
240 AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
241 chi2 += localCluster->GetLocalChi2();
242 }
243 chi2 /= nMFTClusters;
244
245 // now comparing the tracks in order to find the best one
246
247 if (chi2<bestChi2 || bestChi2<0) {
248 bestChi2 = chi2;
249 idBestCandidate = iFinalCandidate;
250 }
251
252 }
253
254 if (nFinalTracks) {
255
256 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
257 newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
258
259 //----------------------- Save the information to the AliESDMuonForwardTrack object
260
4608b85e
AU
261// AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
262// myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
263// myESDTrack -> SetChi2(newTrack->GetGlobalChi2());
264// myESDTrack -> SetCharge(newTrack->GetCharge());
265// myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
4cc511f4
AU
266
267 //---------------------------------------------------------------------------------
268
269 new ((*muonForwardTracks)[muonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
270
271 }
272
273 fCandidateTracks->Delete();
274 fFinalBestCandidate = NULL;
275
276 }
277
278 Int_t myEventID = 0;
279 TFile *outputFileMuonGlobalTracks = new TFile("MuonGlobalTracks.root", "update");
280 while (outputFileMuonGlobalTracks->cd(Form("Event%d",myEventID))) myEventID++;
281 outputFileMuonGlobalTracks -> mkdir(Form("Event%d",myEventID));
282 outputFileMuonGlobalTracks -> cd(Form("Event%d",myEventID));
283 outputTreeMuonGlobalTracks -> Write();
284 outputFileMuonGlobalTracks -> Close();
285
286 return 0;
287
288}
289
290//=========================================================================================================================================
291
292void AliMFTTrackerMU::SeparateFrontBackClusters() {
293
294 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
295 fMFTClusterArrayFront[iPlane]->Delete();
296 fMFTClusterArrayBack[iPlane] ->Delete();
297 for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
298 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
299 if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
300 new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
301 }
302 else {
303 new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
304 }
305 }
306 }
307
308}
309
310//==========================================================================================================================================
311
312Int_t AliMFTTrackerMU::FindClusterInPlane(Int_t planeId) {
313
314 AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId));
315
316 // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
317
318 // propagate track to plane #planeId (both to front and back active sensors)
319 // look for compatible clusters
320 // update TrackParam at found cluster (if any) using Kalman Filter
321
322 AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
323
324 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
325 currentParamFront = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
326 currentParamBack = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
327 currentParamForResearchFront = currentParamFront;
328 currentParamForResearchBack = currentParamBack;
329 Double_t xExtrap = gRandom->Gaus(0,fVertexErrorX);
330 Double_t yExtrap = gRandom->Gaus(0,fVertexErrorY);
331 Double_t zExtrap = gRandom->Gaus(0,fVertexErrorZ);
332 if (fBransonCorrection) {
333 AliMUONTrackExtrap::ExtrapToVertex(&currentParamFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
334 AliMUONTrackExtrap::ExtrapToVertex(&currentParamBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
335 }
336 else {
337 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, zExtrap);
338 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack, zExtrap);
339 }
340 AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
341 AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
342 }
343 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
344 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
345 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
346 currentParamForResearchFront = currentParamFront;
347 currentParamForResearchBack = currentParamBack;
348 AliMUONTrackExtrap::AddMCSEffect(&currentParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
349 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
350 AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
351 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
352 AliMUONTrackExtrap::AddMCSEffect(&currentParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
353 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
354 AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
355 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
356 }
357 // for all planes: extrapolation to the Z of the plane
358 AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
359 AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
360 AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
361 AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
362
363 //---------------------------------------------------------------------------------------
364
365 TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
366 TMatrixD covBack(5,5); covBack = currentParamForResearchBack.GetCovariances();
367
368 Double_t squaredError_X_Front = covFront(0,0);
369 Double_t squaredError_Y_Front = covFront(2,2);
370 Double_t squaredError_X_Back = covBack(0,0);
371 Double_t squaredError_Y_Back = covBack(2,2);
372
373 Double_t corrFact = 1.0;
374
375 Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
376 Double_t researchRadiusBack = TMath::Sqrt(squaredError_X_Back + squaredError_Y_Back);
377 if (0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtPlane[planeId]) {
378 corrFact = fMinResearchRadiusAtPlane[planeId]/(0.5*(researchRadiusFront+researchRadiusBack));
379 }
380
381 //---------------------------------------------------------------------------------------
382
383 Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut; // depends on the number of variables (here, 2)
384
385 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
386
387 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
388 AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
389
390 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
391
392 Bool_t isGoodChi2 = kFALSE;
393
394 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster);
395 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster); // describes the compatibility between the track and the cluster
396 if (chi2<chi2cut) isGoodChi2 = kTRUE;
397
398 if (isGoodChi2) {
399 AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
400 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
401 if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
402 newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
403 AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
404 planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
405 newTrack->SetPlaneExists(planeId);
406 }
407 else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
408
409 }
410
411 // Analyizing the clusters: BACK ACTIVE ELEMENTS
412
413 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
414 AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
415
416 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
417
418 Bool_t isGoodChi2 = kFALSE;
419
420 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster);
421 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster); // describes the compatibility between the track and the cluster
422 if (chi2<chi2cut) isGoodChi2 = kTRUE;
423
424 if (isGoodChi2) {
425 AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
426 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
427 if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
428 newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
429 AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
430 planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
431 newTrack->SetPlaneExists(planeId);
432 }
433 else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
434
435 }
436
437 //---------------------------------------------------------------------------------------------
438
439 return kConverged;
440
441}
442
443//==========================================================================================================================================
444
445Double_t AliMFTTrackerMU::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
446
447 // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
448 // return the corresponding Chi2
449 // assume the track parameters are given at the Z of the cluster
450
451 // Set differences between trackParam and cluster in the bending and non bending directions
452 Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
453 Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
454 AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
455
456 // Calculate errors and covariances
457 const TMatrixD& kParamCov = trackParam.GetCovariances();
458 Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
459 Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
460 AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
461 Double_t covXY = kParamCov(0,2);
462 Double_t det = sigmaX2 * sigmaY2 - covXY * covXY;
463
464 // Compute chi2
465 if (det==0.) return 1.e10;
466 return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
467
468}
469
470//=========================================================================================================================================
471