Config.C updated
[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 "AliMFT.h"
37 #include "AliMUONTrackExtrap.h"
38 #include "AliMUONTrack.h"
39 #include "AliMUONESDInterface.h"
40 #include "AliMuonForwardTrack.h"
41
42 ClassImp(AliMFTTrackerMU)
43
44 const Double_t AliMFTTrackerMU::fRadLengthSi = AliMFTConstants::fRadLengthSi;
45
46 //====================================================================================================================================================
47
48 AliMFTTrackerMU::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
89 AliMFTTrackerMU::~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
103 Int_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
122 void 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
138 Int_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
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());
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
292 void 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
312 Int_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
445 Double_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