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