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