]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MFT/AliMFTTrackerMU.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
MFT methods for muon and dimuon offset modified
[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
581 AliRunLoader *fRunLoader = AliRunLoader::Open("galice.root");
582 if (!fRunLoader) {
583 AliError("no run loader found in file galice.root");
584 return;
585 }
586
587 fRunLoader->CdGAFile();
588 fRunLoader->LoadgAlice();
589 fRunLoader->LoadHeader();
590 fRunLoader->GetEvent(gAlice->GetEvNumber());
591
592 TArrayF vtx(3);
593 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vtx);
594
595 fXExtrapVertex = gRandom->Gaus(vtx[0], AliMFTConstants::fPrimaryVertexResX);
596 fYExtrapVertex = gRandom->Gaus(vtx[1], AliMFTConstants::fPrimaryVertexResY);
597 fZExtrapVertex = gRandom->Gaus(vtx[2], AliMFTConstants::fPrimaryVertexResZ);
598 fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
599 fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
600 AliInfo(Form("Set ESD vertex from MC (%f +/- %f, %f +/- %f, %f)",
601 fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
602
603}
604
605//======================================================================================================================================
606
7bafb008 607void AliMFTTrackerMU::AddClustersFromUnderlyingEvent() {
608
609 AliInfo("Adding clusters from underlying event");
610
611 if (!fMFT) return;
612
613 TGrid::Connect("alien://");
614
615 TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForUnderlyingEvent());
616 if (!fileWithClustersToAdd) return;
617 if (!(fileWithClustersToAdd->IsOpen())) return;
618 if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetUnderlyingEventID())))) return;
619
620 TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
621 TTree *treeIn = 0;
622
623 for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
624
625 treeIn = (TTree*) gDirectory->Get("TreeR");
626
627 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
628 if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
629 for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
630 return;
631 }
632 else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
633 }
634
635 treeIn -> GetEntry(0);
636
637 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
638 printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
639 Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
640 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
641 AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
642 for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
643 new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
644 }
645 printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
646 }
647
648 for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
649
650}
651
652//======================================================================================================================================
653
654void AliMFTTrackerMU::AddClustersFromPileUpEvents() {
655
656 AliInfo("Adding clusters from pile-up event(s)");
657
658 if (!fMFT) return;
659
660 TGrid::Connect("alien://");
661
662 TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForPileUpEvents());
663 if (!fileWithClustersToAdd) return;
664 if (!(fileWithClustersToAdd->IsOpen())) return;
665
666 TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
667 TTree *treeIn = 0;
668
669 for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) {
670
671 if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetPileUpEventID(iPileUp))))) continue;
672
673 for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
674
675 treeIn = (TTree*) gDirectory->Get("TreeR");
676
677 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
678 if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
679 for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
680 return;
681 }
682 else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
683 }
684
685 treeIn -> GetEntry(0);
686
687 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
688 printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
689 Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
690 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
691 AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
692 for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
693 new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
694 }
695 printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
696 }
697
698 for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
699
700 }
701
702}
703
704//======================================================================================================================================
705