Memory leak fixed in AliMFTAnalysisTools
[u/mrichter/AliRoot.git] / MFT / AliMFTTrackerMU.cxx
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"
32 #include "AliESDMuonGlobalTrack.h"
33 #include "AliMFTTrackerMU.h"
34 #include "TMath.h"
35 #include "AliRun.h"
36 #include "AliRunLoader.h"
37 #include "AliHeader.h"
38 #include "AliGenEventHeader.h"
39 #include "TArrayF.h"
40 #include "AliMFT.h"
41 #include "AliMUONTrackExtrap.h"
42 #include "AliMUONTrack.h"
43 #include "AliMUONESDInterface.h"
44 #include "AliMuonForwardTrack.h"
45 #include "AliMUONConstants.h"
46
47 ClassImp(AliMFTTrackerMU)
48
49 const Double_t AliMFTTrackerMU::fRadLengthSi = AliMFTConstants::fRadLengthSi;
50
51 //====================================================================================================================================================
52
53 AliMFTTrackerMU::AliMFTTrackerMU() : 
54   AliTracker(),
55   fESD(0),
56   fMFT(0),
57   fSegmentation(0),
58   fNPlanesMFT(0),
59   fNPlanesMFTAnalyzed(0),
60   fSigmaClusterCut(2),
61   fScaleSigmaClusterCut(1.),
62   fNMaxMissingMFTClusters(0),
63   fGlobalTrackingDiverged(kFALSE),
64   fCandidateTracks(0),
65   fMUONTrack(0),
66   fCurrentTrack(0),
67   fFinalBestCandidate(0),
68   fXExtrapVertex(0),
69   fYExtrapVertex(0),
70   fZExtrapVertex(0),
71   fXExtrapVertexError(0),
72   fYExtrapVertexError(0),
73   fBransonCorrection(kFALSE)
74 {
75
76   //--------------------------------------------------------------------
77   // This is the AliMFTTrackerMU constructor
78   //--------------------------------------------------------------------
79
80   fMFT = (AliMFT*) gAlice->GetDetector("MFT");
81   fSegmentation = fMFT->GetSegmentation();
82   SetNPlanesMFT(fSegmentation->GetNPlanes());
83   AliMUONTrackExtrap::SetField();                 // set the magnetic field for track extrapolations
84
85   AliInfo(Form("fMFT = %p, fSegmentation = %p", fMFT, fSegmentation));
86
87   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
88     fMFTClusterArray[iPlane]      = new TClonesArray("AliMFTCluster");
89     fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
90     fMFTClusterArrayBack[iPlane]  = new TClonesArray("AliMFTCluster");
91     fMFTClusterArray[iPlane]      -> SetOwner(kTRUE);
92     fMFTClusterArrayFront[iPlane] -> SetOwner(kTRUE);
93     fMFTClusterArrayBack[iPlane]  -> SetOwner(kTRUE);
94     fMinResearchRadiusAtPlane[iPlane] = 0.;
95   }
96
97   fCandidateTracks = new TClonesArray("AliMuonForwardTrack",50000);
98
99 }
100
101 //====================================================================================================================================================
102
103 AliMFTTrackerMU::~AliMFTTrackerMU() {
104
105   // destructor
106
107   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
108     fMFTClusterArray[iPlane] -> Delete();
109     delete fMFTClusterArray[iPlane];
110     delete fMFTClusterArrayFront[iPlane];
111     delete fMFTClusterArrayBack[iPlane];
112   }
113
114   delete fCandidateTracks;
115
116 }
117
118 //====================================================================================================================================================
119
120 Int_t AliMFTTrackerMU::LoadClusters(TTree *cTree) {
121
122   //--------------------------------------------------------------------
123   // This function loads the MFT clusters
124   //--------------------------------------------------------------------
125  
126   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
127     AliDebug(1, Form("Setting Address for Branch Plane_%02d", iPlane)); 
128     cTree->SetBranchAddress(Form("Plane_%02d",iPlane), &fMFTClusterArray[iPlane]);
129   }
130  
131   if (!cTree->GetEvent()) return kFALSE;
132   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
133     AliInfo(Form("plane %02d: nClusters = %d\n", iPlane, fMFTClusterArray[iPlane]->GetEntries()));
134   }
135   SeparateFrontBackClusters();
136
137   return 0;
138
139 }
140
141 //====================================================================================================================================================
142
143 void AliMFTTrackerMU::UnloadClusters() {
144   
145   //--------------------------------------------------------------------
146   // This function unloads MFT clusters
147   //--------------------------------------------------------------------
148
149   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
150     fMFTClusterArray[iPlane]      -> Clear("C");
151     fMFTClusterArrayFront[iPlane] -> Clear("C");
152     fMFTClusterArrayBack[iPlane]  -> Clear("C");
153   }
154
155 }
156
157 //====================================================================================================================================================
158
159 Int_t AliMFTTrackerMU::Clusters2Tracks(AliESDEvent *event) {
160
161   //--------------------------------------------------------------------
162   // This functions reconstructs the Muon Forward Tracks
163   // The clusters must be already loaded !
164   //--------------------------------------------------------------------
165
166   // Tree of AliMuonForwardTrack objects. Created outside the ESD framework for cross-check purposes
167
168   TFile *outputFileMuonGlobalTracks = new TFile("MuonGlobalTracks.root", "update");
169   TTree *outputTreeMuonGlobalTracks = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
170   TClonesArray *muonForwardTracks = new TClonesArray("AliMuonForwardTrack");
171   outputTreeMuonGlobalTracks -> Branch("tracks", &muonForwardTracks);
172  
173   //--------------------------------------------------------------------
174
175   fESD = event;
176
177   fXExtrapVertex = 0;
178   fYExtrapVertex = 0;
179   fZExtrapVertex = 0;
180   fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
181   fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
182
183   // Imposing a fixed primary vertex, in order to keep memory of its position. Taking the primary vertex would imply the risk
184   // of loosing memory of its position when passing from ESD to AOD, due to possible refitting
185   const AliESDVertex* esdVert = fESD->GetVertex(); 
186   if (esdVert->GetNContributors() > 0 || !strcmp(esdVert->GetTitle(),"vertexer: smearMC")) {
187     fXExtrapVertex = esdVert->GetX();
188     fYExtrapVertex = esdVert->GetY();
189     fZExtrapVertex = esdVert->GetZ();
190     fXExtrapVertexError = TMath::Max(AliMFTConstants::fXVertexTolerance, esdVert->GetXRes());
191     fYExtrapVertexError = TMath::Max(AliMFTConstants::fYVertexTolerance, esdVert->GetYRes());
192     AliInfo(Form("Found ESD vertex from %d contributors (%f +/- %f,  %f +/- %f,  %f)", 
193                  esdVert->GetNContributors(),fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
194   }
195   else GetVertexFromMC();    
196
197   //----------- Read ESD MUON tracks -------------------
198
199   Int_t nTracksMUON = event->GetNumberOfMuonTracks();
200
201   AliInfo(Form("Number of ESD MUON tracks: %d\n", nTracksMUON));
202
203   Int_t iTrack=0;
204   while (iTrack<nTracksMUON) {
205
206     fNPlanesMFTAnalyzed = 0;
207
208     AliInfo("****************************************************************************************");
209     AliInfo(Form("***************************   MUON TRACK %3d/%d   ***************************************", iTrack, nTracksMUON));
210     AliInfo("****************************************************************************************");
211     
212     fCandidateTracks -> Delete();
213     
214     fNPlanesMFTAnalyzed = 0;
215     
216     const AliESDMuonTrack *esdTrack = event->GetMuonTrack(iTrack);
217     if (fMUONTrack) delete fMUONTrack;
218     fMUONTrack = new AliMUONTrack();
219
220     AliMUONESDInterface::ESDToMUON(*esdTrack, *fMUONTrack, kFALSE);
221
222     if (!fMUONTrack->GetTrackParamAtCluster()->First()) {
223       AliInfo("Skipping track, no parameters available!!!");
224       iTrack++;
225       continue;
226     }
227
228     // the track we are going to build, starting from fMUONTrack and adding the MFT clusters
229     AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
230     track -> SetMUONTrack(new AliMUONTrack(*fMUONTrack));
231     track -> SetMCLabel(fMUONTrack->GetMCLabel());
232     track -> SetMatchTrigger(fMUONTrack->GetMatchTrigger());
233
234     // track parameters linearly extrapolated from the first tracking station to the end of the absorber
235     AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
236     AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, AliMUONConstants::AbsZEnd());   // absorber extends from -90 to -503 cm
237     Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
238     Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
239     Double_t rAbsorber      = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
240     track -> SetRAtAbsorberEnd(rAbsorber);
241
242     //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
243
244     for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) {   /* *** do not reverse the order of this cycle!!! 
245                                                                  *** this reflects the fact that the extrapolation is performed 
246                                                                  *** starting from the last MFT plane back to the origin */
247       
248       // --------- updating the array of candidates according to the clusters available in the i-th plane ---------
249       
250       fNPlanesMFTAnalyzed++;
251       
252       Int_t nCandidates = fCandidateTracks->GetEntriesFast();
253       for (Int_t iCandidate=0; iCandidate<nCandidates; iCandidate++) {
254         fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iCandidate);
255
256         // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array 
257         // (several new tracks can be created for one old track)
258         if (FindClusterInPlane(iPlane) == kDiverged) {
259           fGlobalTrackingDiverged = kTRUE;
260           break;
261         }
262
263         if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
264           fCandidateTracks->Remove(fCurrentTrack);     // the old track is removed after the check;
265         }
266       }
267       if (fGlobalTrackingDiverged) {
268         if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
269         continue;
270       }
271
272       fCandidateTracks->Compress();
273       
274     }      
275
276     // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
277     
278     fGlobalTrackingDiverged = kFALSE;
279     fScaleSigmaClusterCut = 1.0;
280     
281     AliDebug(1, "Finished cycle over planes");
282     
283     iTrack++;
284
285     // If we have several final tracks, we must find the best candidate:
286     
287     Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
288     AliInfo(Form("nFinalTracks = %d", nFinalTracks));
289
290     Int_t nGoodClustersBestCandidate =  0;
291     Int_t idBestCandidate            =  0;
292     Double_t bestChi2                = -1.;  // variable defining the best candidate
293
294     for (Int_t iFinalCandidate=0; iFinalCandidate<nFinalTracks; iFinalCandidate++) {
295       
296       AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iFinalCandidate);
297       Int_t nMFTClusters  = finalTrack->GetNMFTClusters();
298
299       Double_t chi2 = 0;
300       for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
301         AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
302         chi2 += localCluster->GetLocalChi2();
303       }
304       chi2 /= nMFTClusters;
305
306       // now comparing the tracks in order to find the best one
307       
308       if (chi2<bestChi2 || bestChi2<0) {
309         bestChi2 = chi2;
310         idBestCandidate = iFinalCandidate;
311       }
312       
313     }
314     
315     if (nFinalTracks) {
316
317       AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
318       newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
319
320       new ((*muonForwardTracks)[muonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
321
322       //----------------------- Save the information to the AliESDMuonGlobalTrack object
323
324       newTrack -> EvalKinem(AliMFTConstants::fZEvalKinem);
325
326       AliDebug(2,"Creating a new Muon Global Track");
327       AliESDMuonGlobalTrack *myESDTrack = event->NewMuonGlobalTrack();
328       myESDTrack -> SetPxPyPz(newTrack->Px(), newTrack->Py(), newTrack->Pz());
329       
330       myESDTrack -> SetLabel(newTrack->GetMCLabel());
331       myESDTrack -> SetChi2OverNdf(newTrack->GetChi2OverNdf());
332       myESDTrack -> SetCharge(newTrack->GetCharge());
333       myESDTrack -> SetMatchTrigger(newTrack->GetMatchTrigger());
334       myESDTrack -> SetNMFTClusters(newTrack->GetNMFTClusters());
335       myESDTrack -> SetNWrongMFTClustersMC(newTrack->GetNWrongClustersMC());
336       myESDTrack -> SetFirstTrackingPoint(newTrack->GetMFTCluster(0)->GetX(), newTrack->GetMFTCluster(0)->GetY(), newTrack->GetMFTCluster(0)->GetZ());
337       myESDTrack -> SetXYAtVertex(newTrack->GetOffsetX(0., AliMFTConstants::fZEvalKinem), newTrack->GetOffsetY(0., AliMFTConstants::fZEvalKinem));
338       myESDTrack -> SetRAtAbsorberEnd(newTrack->GetRAtAbsorberEnd());
339       myESDTrack -> SetCovariances(newTrack->GetTrackParamAtMFTCluster(0)->GetCovariances());
340       myESDTrack -> SetChi2MatchTrigger(esdTrack->GetChi2MatchTrigger());
341       myESDTrack -> SetMuonClusterMap(esdTrack->GetMuonClusterMap());
342       myESDTrack -> SetHitsPatternInTrigCh(esdTrack->GetHitsPatternInTrigCh());
343       myESDTrack -> SetHitsPatternInTrigChTrk(esdTrack->GetHitsPatternInTrigChTrk());
344       myESDTrack -> Connected(esdTrack->IsConnected());
345
346       ULong_t mftClusterPattern = 0;
347       for (Int_t iCluster=0; iCluster<newTrack->GetNMFTClusters(); iCluster++) {
348         AliMFTCluster *localCluster = newTrack->GetMFTCluster(iCluster);
349         mftClusterPattern |= 1 << localCluster->GetPlane();
350         mftClusterPattern |= IsCorrectMatch(localCluster, newTrack->GetMCLabel()) << (fNMaxPlanes + localCluster->GetPlane());
351       }
352       myESDTrack -> SetMFTClusterPattern(mftClusterPattern);
353       
354       //---------------------------------------------------------------------------------
355
356     }
357
358     fFinalBestCandidate = NULL;
359    
360   }
361
362   outputTreeMuonGlobalTracks -> Fill();
363
364   Int_t myEventID = 0;
365   while (outputFileMuonGlobalTracks->cd(Form("Event%d",myEventID))) myEventID++;
366   outputFileMuonGlobalTracks -> mkdir(Form("Event%d",myEventID));
367   outputFileMuonGlobalTracks -> cd(Form("Event%d",myEventID));
368   outputTreeMuonGlobalTracks -> Write();
369   outputFileMuonGlobalTracks -> Close();
370
371   muonForwardTracks -> Delete();
372   delete muonForwardTracks;
373
374   return 0;
375
376 }
377
378 //=========================================================================================================================================
379
380 void AliMFTTrackerMU::SeparateFrontBackClusters() {
381
382   for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
383     fMFTClusterArrayFront[iPlane]->Delete();
384     fMFTClusterArrayBack[iPlane] ->Delete();
385     for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
386       AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
387       if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
388         new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
389       }
390       else {
391         new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
392       }
393     }
394   }
395
396 }
397
398 //==========================================================================================================================================
399
400 Int_t AliMFTTrackerMU::FindClusterInPlane(Int_t planeId) { 
401   
402   // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
403
404   // propagate track to plane #planeId (both to front and back active sensors)
405   // look for compatible clusters
406   // update TrackParam at found cluster (if any) using Kalman Filter
407
408   AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
409
410   if (planeId == fNPlanesMFT-1) {      // last plane of the telecope
411     currentParamFront = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
412     currentParamBack  = (*((AliMUONTrackParam*)(fMUONTrack->GetTrackParamAtCluster()->First())));
413     currentParamForResearchFront = currentParamFront;
414     currentParamForResearchBack  = currentParamBack;
415     if (fBransonCorrection) {
416       AliMUONTrackExtrap::ExtrapToVertex(&currentParamFront, fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
417       AliMUONTrackExtrap::ExtrapToVertex(&currentParamBack,  fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
418     }
419     else {
420       AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, fZExtrapVertex);
421       AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack,  fZExtrapVertex);
422     }
423     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
424     AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack,  fXExtrapVertex, fYExtrapVertex, fZExtrapVertex, fXExtrapVertexError, fYExtrapVertexError); 
425   }
426   else {          // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
427     currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
428     currentParamBack  = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
429     currentParamForResearchFront = currentParamFront;
430     currentParamForResearchBack  = currentParamBack;
431     AliMUONTrackExtrap::AddMCSEffect(&currentParamFront,           (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
432                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
433     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
434                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
435     AliMUONTrackExtrap::AddMCSEffect(&currentParamBack,            (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
436                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
437     AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
438                                                                     fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
439   }
440   // for all planes: extrapolation to the Z of the plane
441   AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront,            -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());   
442   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
443   AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack,             -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());   
444   AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack,  -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
445
446   //---------------------------------------------------------------------------------------
447
448   TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
449   TMatrixD covBack(5,5);  covBack  = currentParamForResearchBack.GetCovariances();
450   
451   Double_t squaredError_X_Front = covFront(0,0);
452   Double_t squaredError_Y_Front = covFront(2,2);
453   Double_t squaredError_X_Back  = covBack(0,0);
454   Double_t squaredError_Y_Back  = covBack(2,2);
455
456   Double_t corrFact = 1.0;
457
458   Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
459   Double_t researchRadiusBack  = TMath::Sqrt(squaredError_X_Back  + squaredError_Y_Back);
460   if (0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtPlane[planeId]) {
461     corrFact = fMinResearchRadiusAtPlane[planeId]/(0.5*(researchRadiusFront+researchRadiusBack));
462   }
463
464   //---------------------------------------------------------------------------------------
465
466   Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut;     // depends on the number of variables (here, 2)
467   
468   // Analyizing the clusters: FRONT ACTIVE ELEMENTS
469   
470   Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
471   AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
472   
473   for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
474
475     Bool_t isGoodChi2 = kFALSE;
476
477     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster); 
478     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster);     // describes the compatibility between the track and the cluster
479     if (chi2<chi2cut) isGoodChi2 = kTRUE;
480
481     if (isGoodChi2) {
482       AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
483       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
484       if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
485       newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster);    // creating new track param and attaching the cluster
486       AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
487                        planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
488       newTrack->SetPlaneExists(planeId);
489     }
490     else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
491
492   }
493
494   // Analyizing the clusters: BACK ACTIVE ELEMENTS
495   
496   Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
497   AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
498   
499   for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
500
501     Bool_t isGoodChi2 = kFALSE;
502
503     AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster); 
504     Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster);     // describes the compatibility between the track and the cluster
505     if (chi2<chi2cut) isGoodChi2 = kTRUE;
506
507     if (isGoodChi2) {
508       AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
509       AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
510       if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
511       newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster);    // creating new track param and attaching the cluster
512       AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", 
513                        planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
514       newTrack->SetPlaneExists(planeId);
515     }
516     else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
517
518   }
519
520   //---------------------------------------------------------------------------------------------
521
522   return kConverged;
523   
524 }
525
526 //==========================================================================================================================================
527
528 Double_t AliMFTTrackerMU::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
529
530   // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
531   // return the corresponding Chi2
532   // assume the track parameters are given at the Z of the cluster
533   
534   // Set differences between trackParam and cluster in the bending and non bending directions
535   Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
536   Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
537   AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
538   
539   // Calculate errors and covariances
540   const TMatrixD& kParamCov = trackParam.GetCovariances();
541   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
542   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
543   AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
544   Double_t covXY   = kParamCov(0,2);
545   Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
546   
547   // Compute chi2
548   if (det==0.) return 1.e10;
549   return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
550   
551 }
552
553 //=========================================================================================================================================
554
555 Bool_t AliMFTTrackerMU::IsCorrectMatch(AliMFTCluster *cluster, Int_t labelMC) {
556
557   Bool_t result = kFALSE;
558
559   // check if the cluster belongs to the correct MC track
560
561   for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
562     if (cluster->GetMCLabel(iTrack)==labelMC) {
563       result = kTRUE;
564       break;
565     }
566   }
567
568   return result;
569
570 }
571
572 //======================================================================================================================================
573
574 void AliMFTTrackerMU::GetVertexFromMC() {
575
576   AliRunLoader *fRunLoader = AliRunLoader::Open("galice.root");
577   if (!fRunLoader) {
578     AliError("no run loader found in file galice.root");
579     return;
580   }
581
582   fRunLoader->CdGAFile();
583   fRunLoader->LoadgAlice();
584   fRunLoader->LoadHeader();
585   fRunLoader->GetEvent(gAlice->GetEvNumber());
586   
587   TArrayF vtx(3);
588   fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vtx);
589
590   fXExtrapVertex = gRandom->Gaus(vtx[0], AliMFTConstants::fPrimaryVertexResX);
591   fYExtrapVertex = gRandom->Gaus(vtx[1], AliMFTConstants::fPrimaryVertexResY);
592   fZExtrapVertex = gRandom->Gaus(vtx[2], AliMFTConstants::fPrimaryVertexResZ);
593   fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
594   fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
595   AliInfo(Form("Set ESD vertex from MC (%f +/- %f,  %f +/- %f,  %f)", 
596                fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
597   
598 }
599
600 //======================================================================================================================================
601