- Add simple QA to the trigger patch maker and set the default to
[u/mrichter/AliRoot.git] / MFT / AliMFTTrackerMU.cxx
... / ...
CommitLineData
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
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),
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
104AliMFTTrackerMU::~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
121Int_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
148void 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
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
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
388void 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
408Int_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
536Double_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
563Bool_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
582void 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
611void 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
658void 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