3 #include "TClonesArray.h"
8 #include "TGeoManager.h"
10 #include "TParticle.h"
20 #include "TIterator.h"
25 #include "AliRunLoader.h"
26 #include "AliLoader.h"
27 #include "AliHeader.h"
31 #include "AliTracker.h"
32 #include "AliGRPObject.h"
33 #include "AliCDBEntry.h"
34 #include "AliCDBManager.h"
37 #include "AliMUONConstants.h"
38 #include "AliMUONTrack.h"
39 #include "AliMUONRecoCheck.h"
40 #include "AliMUONTrackParam.h"
41 #include "AliMUONTrackExtrap.h"
42 #include "AliMUONVTrackStore.h"
43 #include "AliMUONVCluster.h"
46 #include "AliMuonForwardTrack.h"
47 #include "AliMFTCluster.h"
49 #include "AliMFTSegmentation.h"
50 #include "AliMFTConstants.h"
52 #include "AliMuonForwardTrackFinder.h"
54 //====================================================================================================================================================
56 // Class for the creation of the "global muon tracks" built from the clusters in the
57 // muon spectrometer and the clusters of the Muon Forward Tracker. QA histograms are also created
59 // Contact author: antonio.uras@cern.ch
61 //====================================================================================================================================================
63 const Double_t AliMuonForwardTrackFinder::fRadLengthSi = AliMFTConstants::fRadLengthSi;
65 ClassImp(AliMuonForwardTrackFinder)
67 //=====================================================================================================
69 AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
74 fScaleSigmaClusterCut(1.),
75 fGlobalTrackingDiverged(kFALSE),
77 fSigmaSpectrometerCut(0),
81 fNFinalCandidatesCut(0),
86 fDistanceFromGoodClusterAndTrackAtLastPlane(-1),
87 fDistanceFromBestClusterAndTrackAtLastPlane(-1),
92 fNPlanesMFTAnalyzed(0),
93 fNMaxMissingMFTClusters(0),
98 fHistRadiusEndOfAbsorber(0),
99 fHistNGoodClustersForFinalTracks(0),
100 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane(0),
101 fHistDistanceGoodClusterFromTrackAtLastPlane(0),
103 fNtuFinalCandidates(0),
104 fNtuFinalBestCandidates(0),
109 fTxtTrackGoodClusters(0),
110 fTxtTrackFinalChi2(0),
111 fTxtTrackMomentum(0),
112 fTxtFinalCandidates(0),
115 fTxtClustGoodChi2(0),
119 fMrkClustGoodChi2(0),
123 fCountRealTracksAnalyzed(0),
124 fMaxNTracksToBeAnalyzed(99999999),
125 fCountRealTracksWithRefMC(0),
126 fCountRealTracksWithRefMC_andTrigger(0),
127 fCountRealTracksWithRefMC_andTrigger_andGoodPt(0),
128 fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta(0),
129 fCountRealTracksAnalyzedOfEvent(0),
130 fCountRealTracksAnalyzedWithFinalCandidates(0),
142 fFinalBestCandidate(0),
143 fIsCurrentMuonTrackable(0),
154 fMuonForwardTracks(0),
158 fBransonCorrection(kTRUE)
162 // Default constructor
164 for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
166 fHistNTracksAfterExtrapolation[iPlane] = 0;
167 fHistChi2Cluster_GoodCluster[iPlane] = 0;
168 fHistChi2Cluster_BadCluster[iPlane] = 0;
169 fHistResearchRadius[iPlane] = 0;
171 fIsGoodClusterInPlane[iPlane] = kFALSE;
173 fHistChi2Cluster_GoodCluster[iPlane] = 0;
174 fHistChi2Cluster_BadCluster[iPlane] = 0;
176 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = 0;
177 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = 0;
179 fZPlane[iPlane] = 0.;
180 fRPlaneMax[iPlane] = 0.;
181 fRPlaneMin[iPlane] = 0.;
183 for (Int_t i=0; i<4; i++) fGrMFTPlane[i][iPlane] = 0;
184 fCircleExt[iPlane] = 0;
185 fCircleInt[iPlane] = 0;
187 fTxtTrackChi2[iPlane] = 0;
189 fIsClusterCompatible[iPlane] = 0;
191 fMFTClusterArray[iPlane] = 0;
192 fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
193 fMFTClusterArrayBack[iPlane] = new TClonesArray("AliMFTCluster");
195 fIsPlaneMandatory[iPlane] = kFALSE;
197 fMinResearchRadiusAtPlane[iPlane] = 0.;
203 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane = 0;
204 fHistDistanceGoodClusterFromTrackAtLastPlane = 0;
207 fCandidateTracks = 0;
209 fOutputTreeFile = new TFile("MuonGlobalTracks.root", "recreate");
210 fOutputEventTree = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
211 fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack");
212 fOutputEventTree -> Branch("tracks", &fMuonForwardTracks);
216 //=====================================================================================================
218 AliMuonForwardTrackFinder::~AliMuonForwardTrackFinder() {
220 for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) {
222 delete fHistNTracksAfterExtrapolation[iPlane];
223 delete fHistChi2Cluster_GoodCluster[iPlane];
224 delete fHistChi2Cluster_BadCluster[iPlane];
225 delete fHistResearchRadius[iPlane];
227 delete fHistChi2Cluster_GoodCluster[iPlane];
228 delete fHistChi2Cluster_BadCluster[iPlane];
230 delete fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane];
231 delete fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane];
233 for (Int_t i=0; i<4; i++) delete fGrMFTPlane[i][iPlane];
234 delete fCircleExt[iPlane];
235 delete fCircleInt[iPlane];
237 delete fTxtTrackChi2[iPlane];
239 delete fMFTClusterArray[iPlane];
240 delete fMFTClusterArrayFront[iPlane];
241 delete fMFTClusterArrayBack[iPlane];
245 delete fNtuFinalCandidates;
246 delete fNtuFinalBestCandidates;
248 delete fHistRadiusEndOfAbsorber;
250 delete fHistNGoodClustersForFinalTracks;
251 delete fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane; //
252 delete fHistDistanceGoodClusterFromTrackAtLastPlane; //
256 delete fTxtMuonHistory;
257 delete fTxtTrackGoodClusters;
258 delete fTxtTrackFinalChi2;
259 delete fTxtTrackMomentum;
260 delete fTxtFinalCandidates;
263 delete fTxtClustGoodChi2;
265 delete fTxtClustOfTrack;
267 delete fMrkClustGoodChi2;
269 delete fMrkClustOfTrack;
277 delete fMuonRecoCheck;
279 delete fMFTClusterTree;
281 delete fMuonTrackReco;
282 delete fCurrentTrack;
283 delete fFinalBestCandidate;
285 delete fCandidateTracks;
288 delete fTrackRefStore;
295 delete fSegmentation;
297 delete fOutputTreeFile;
298 delete fOutputQAFile;
299 delete fOutputEventTree;
301 delete fMuonForwardTracks;
308 //=====================================================================================================
310 void AliMuonForwardTrackFinder::Init(Int_t nRun,
313 Int_t nEventsToAnalyze) {
316 AliInfo("WARNING: run already initialized!!\n");
323 AliInfo(Form("input dir = %s\n", fReadDir.Data()));
324 AliInfo(Form("output dir = %s\n", fOutDir.Data()));
326 // -------------------------- initializing files...
328 AliInfo(Form("initializing files for run %d...\n", fRun));
330 Char_t geoFileName[300];
331 Char_t esdFileName[300];
332 Char_t gAliceName[300];
333 Char_t clusterName[300];
335 snprintf(geoFileName , 300, "%s/geometry.root", fReadDir.Data());
336 snprintf(esdFileName , 300, "%s/AliESDs.root" , fReadDir.Data());
337 snprintf(gAliceName , 300, "%s/galice.root" , fReadDir.Data());
338 snprintf(clusterName , 300, "%s/MFT.RecPoints.root", fReadDir.Data());
340 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
342 TGeoManager::Import(geoFileName);
344 AliError(Form("getting geometry from file %s failed", geoFileName));
349 fFileESD = new TFile(esdFileName);
350 if (!fFileESD || !fFileESD->IsOpen()) return;
351 else AliInfo(Form("file %s successfully opened\n", fFileESD->GetName()));
353 fMuonRecoCheck = new AliMUONRecoCheck(esdFileName, Form("%s/generated/", fReadDir.Data())); // Utility class to check reconstruction
354 fFile_gAlice = new TFile(gAliceName);
355 if (!fFile_gAlice || !fFile_gAlice->IsOpen()) return;
356 else AliInfo(Form("file %s successfully opened\n", fFile_gAlice->GetName()));
358 fRunLoader = AliRunLoader::Open(gAliceName);
359 gAlice = fRunLoader->GetAliRun();
360 if (!gAlice) fRunLoader->LoadgAlice();
361 fMFT = (AliMFT*) gAlice->GetDetector("MFT");
362 fSegmentation = fMFT->GetSegmentation();
363 SetNPlanesMFT(fSegmentation->GetNPlanes());
365 if (!SetRunNumber()) return;
366 if (!InitGRP()) return;
367 AliMUONTrackExtrap::SetField(); // set the magnetic field for track extrapolations
369 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
370 fZPlane[iPlane] = fSegmentation->GetPlane(iPlane)->GetZCenter();
371 fRPlaneMax[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMaxSupport();
372 fRPlaneMin[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMinSupport();
375 // Loading MFT clusters
376 fMFTLoader = fRunLoader->GetDetectorLoader("MFT");
377 fMFTLoader->LoadRecPoints("READ");
379 fMFTClusterTree = fMFTLoader->TreeR();
381 Int_t nEventsInFile = fMuonRecoCheck->NumberOfEvents();
382 if (!nEventsInFile) {
383 AliError("no events available!!!\n");
386 if (nEventsInFile<nEventsToAnalyze || nEventsToAnalyze<0) fNEventsToAnalyze = nEventsInFile;
387 else fNEventsToAnalyze = nEventsToAnalyze;
389 fCandidateTracks = new TClonesArray("AliMuonForwardTrack",50000);
391 // -------------------------- initializing histograms...
393 AliInfo("\ninitializing histograms...\n");
396 AliInfo("... done!\n\n");
398 // -------------------------- initializing graphics...
400 AliInfo("initializing graphics...\n");
402 AliInfo("... done!\n\n");
404 SetSigmaSpectrometerCut(4.0);
405 SetSigmaClusterCut(4.5);
406 SetChi2GlobalCut(2.0);
407 SetNFinalCandidatesCut(10);
408 SetRAbsorberCut(26.4);
413 //======================================================================================================================================
415 Bool_t AliMuonForwardTrackFinder::LoadNextEvent() {
417 // load next reconstructed event from the tree
419 if (fEv) FillOutputTree();
421 if (fEv>=fNEventsToAnalyze) return kFALSE;
423 fCountRealTracksAnalyzedOfEvent = 0;
425 AliInfo(Form(" **** analyzing event # %d \n", fEv));
427 fTrackStore = fMuonRecoCheck->ReconstructedTracks(fEv);
428 if (fTrackStore->IsEmpty()) {
429 AliInfo("fTrackStore Is Empty: exiting NOW!");
432 AliInfo("fTrackStore contains tracks!");
434 AliDebug(2, Form("Getting fMuonRecoCheck->ReconstructibleTracks(%d)", fEv));
435 fTrackRefStore = fMuonRecoCheck->ReconstructibleTracks(fEv);
437 AliDebug(2, Form("Getting fRunLoader->GetEvent(%d)", fEv));
438 fRunLoader->GetEvent(fEv);
440 AliDebug(2, Form("fMFTLoader->TreeR() = %p",fMFTLoader->TreeR()));
441 if (!fMFTLoader->TreeR()->GetEvent()) return kFALSE;
442 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
443 AliDebug(1, Form("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries()));
444 fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
446 SeparateFrontBackClusters();
448 fRunLoader -> LoadKinematics();
449 fStack = fRunLoader->Stack();
450 fNextTrack = fTrackStore->CreateIterator();
451 fMuonForwardTracks->Delete();
459 //======================================================================================================================================
461 Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
463 fNPlanesMFTAnalyzed = 0;
465 // load next muon track from the reconstructed event
467 if (fCountRealTracksAnalyzed>=fMaxNTracksToBeAnalyzed) return kFALSE;
468 if (!fCountRealTracksAnalyzed) if (!LoadNextEvent()) return kFALSE;
470 if (!fGlobalTrackingDiverged) {
471 while ( !(fMuonTrackReco = static_cast<AliMUONTrack*>(fNextTrack->Next())) ) if (!LoadNextEvent()) return kFALSE;
472 fCountRealTracksAnalyzed++;
473 fCountRealTracksAnalyzedOfEvent++;
476 AliDebug(1, "**************************************************************************************\n");
477 AliDebug(1, Form("*************************** MUON TRACK %3d ***************************************\n", fCountRealTracksAnalyzedOfEvent));
478 AliDebug(1, "**************************************************************************************\n");
480 fCandidateTracks -> Delete();
483 fDistanceFromGoodClusterAndTrackAtLastPlane = -1.;
484 fDistanceFromBestClusterAndTrackAtLastPlane = -1.;
487 TIter nextTrackRef(fTrackRefStore->CreateIterator());
488 AliMUONTrack *trackRef=0;
490 // --------------------------------------- loop on MC generated tracks to find the MC reference...
492 while ( (trackRef = static_cast<AliMUONTrack*>(nextTrackRef())) ) {
493 // number of compatible clusters between trackReco and trackRef
494 Int_t nMatchCluster = fMuonTrackReco->FindCompatibleClusters(*trackRef, fSigmaSpectrometerCut, fIsClusterCompatible);
495 if ( (fIsClusterCompatible[0] || fIsClusterCompatible[1] || fIsClusterCompatible[2] || fIsClusterCompatible[3]) && // before the dipole
496 (fIsClusterCompatible[6] || fIsClusterCompatible[7] || fIsClusterCompatible[8] || fIsClusterCompatible[9]) && // after the dipole
497 2*nMatchCluster>fMuonTrackReco->GetNClusters() ) {
498 fMuonTrackReco->SetMCLabel(trackRef->GetUniqueID()); // MC reference has been found for trackReco!
503 // ------------------------------------- ...done!
505 fLabelMC = fMuonTrackReco->GetMCLabel();
509 if (!fGlobalTrackingDiverged) fCountRealTracksWithRefMC++;
510 if (fStack->Particle(fLabelMC)->GetFirstMother() != -1) {
511 motherPdg = fStack->Particle(fStack->Particle(fLabelMC)->GetFirstMother())->GetPdgCode();
515 CheckCurrentMuonTrackable();
517 if (!fGlobalTrackingDiverged) if (fMuonTrackReco->GetMatchTrigger()) fCountRealTracksWithRefMC_andTrigger++;
519 // the track we are going to build, starting from fMuonTrackReco and adding the MFT clusters
520 AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
521 track -> SetMUONTrack(new AliMUONTrack(*fMuonTrackReco));
522 if (fLabelMC>=0 && fStack->Particle(fLabelMC)) track->SetMCTrackRef(new TParticle(*(fStack->Particle(fLabelMC))));
523 track -> SetMCLabel(fMuonTrackReco->GetMCLabel());
524 track -> SetMatchTrigger(fMuonTrackReco->GetMatchTrigger());
527 Double_t xVtx=-999., yVtx=-999., zVtx=-999999.;
528 if (track->GetMCTrackRef()) {
529 xVtx = track->GetMCTrackRef()->Vx();
530 yVtx = track->GetMCTrackRef()->Vy();
531 zVtx = track->GetMCTrackRef()->Vz();
535 Double_t pt=-999., theta=-999., eta=-999.;
536 if (track->GetMCTrackRef()) {
537 pt = track->GetMCTrackRef()->Pt();
538 theta = track->GetMCTrackRef()->Theta();
539 if (theta<0.) theta += TMath::Pi();
540 eta = track->GetMCTrackRef()->Eta();
543 AliMUONTrackParam *param = (AliMUONTrackParam*) (fMuonTrackReco->GetTrackParamAtCluster()->First());
544 pt = TMath::Sqrt(param->Px()*param->Px() + param->Py()*param->Py());
545 theta = TMath::ATan(pt/param->Pz());
546 if (theta<0.) theta += TMath::Pi();
547 eta = -1.*TMath::Log(TMath::Tan(0.5*theta));
549 // if the transverse momentum is smaller than the threshold, skip to the next track
550 if (pt < fLowPtCut) return 3;
552 // track parameters linearly extrapolated from the first tracking station to the end of the absorber
553 AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
554 AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, -503.); // absorber extends from -90 to -503 cm
555 Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
556 Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
557 Double_t rAbsorber = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
558 fHistRadiusEndOfAbsorber -> Fill(rAbsorber);
560 // if the radial distance of the track at the end of the absorber is smaller than a given radius, skip to the next track
561 if (rAbsorber < fRAbsorberCut) return 4;
563 //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
565 for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) { // *** do not reverse the order of this cycle!!!
566 // *** this reflects the fact that the extrapolation is performed
567 // *** starting from the last MFT plane back to the origin
569 // --------- updating the array of tracks according to the clusters available in the i-th plane ---------
571 fNPlanesMFTAnalyzed++;
573 if (fMatchingMode==kRealMatching) {
574 Int_t nTracksToBeAnalyzed = fCandidateTracks->GetEntriesFast();
575 for (Int_t iTrack=0; iTrack<nTracksToBeAnalyzed; iTrack++) {
576 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
577 // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
578 // (several new tracks can be created for one old track)
579 if (FindClusterInPlane(iPlane) == 2) {
580 fGlobalTrackingDiverged = kTRUE;
581 if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
584 if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
585 fCandidateTracks->Remove(fCurrentTrack); // the old track is removed after the check;
588 fCandidateTracks->Compress();
589 if (fIsCurrentMuonTrackable) {
590 // fOutputQAFile->cd();
591 fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
595 else if (fMatchingMode==kIdealMatching && fIsCurrentMuonTrackable) {
596 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
597 AliDebug(2, Form("plane %02d: fCandidateTracks->GetEntriesFast() = %d fCandidateTracks->UncheckedAt(0) = %p fCurrentTrack = %p\n",
598 iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack));
599 AttachGoodClusterInPlane(iPlane);
604 // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
606 fGlobalTrackingDiverged = kFALSE;
607 fScaleSigmaClusterCut = 1.0;
609 AliDebug(1, "Finished cycle over planes");
611 Double_t momentum = pt * TMath::CosH(eta);
612 fTxtTrackMomentum = new TLatex(0.10, 0.70, Form("P_{spectro} = %3.1f GeV/c", momentum));
614 if (fMatchingMode==kIdealMatching) {
615 AliDebug(1, "Adding track to output tree...\n");
616 fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
617 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
618 new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
619 AliDebug(1, "...track added!\n");
620 fCandidateTracks->Delete();
621 fCountRealTracksAnalyzedOfEvent++;
622 fCountRealTracksAnalyzedWithFinalCandidates++;
623 PrintParticleHistory();
624 FillPlanesWithTrackHistory();
626 Double_t chi2AtPlane[fNMaxPlanes] = {0};
627 Int_t nGoodClusters = 0;
628 Int_t nMFTClusters = fFinalBestCandidate->GetNMFTClusters();
629 // Int_t nMUONClusters = fFinalBestCandidate->GetNMUONClusters();
631 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
632 while (!fFinalBestCandidate->PlaneExists(plane)) plane++;
633 AliMFTCluster *localCluster = fFinalBestCandidate->GetMFTCluster(iCluster);
634 chi2AtPlane[plane] = localCluster->GetLocalChi2();
635 if (IsCorrectMatch(localCluster)) nGoodClusters++;
636 // Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope
637 // Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
638 // chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
641 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
642 fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
643 0.90*fRPlaneMax[fNPlanesMFT-1],
644 Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
646 fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2AtPlane[0]));
648 if (fDrawOption) DrawPlanes();
652 // If we have several final tracks, we must find the best candidate:
654 Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
655 AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
657 if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++;
659 Double_t theVariable_Best = -1.; // variable defining the best candidate
660 Bool_t bestCandidateExists = kFALSE;
661 Int_t nGoodClustersBestCandidate = 0;
662 Int_t idBestCandidate = 0;
663 Double_t chi2HistoryForBestCandidate[fNMaxPlanes] = {0}; // chi2 on each plane, for the best candidate
664 Double_t nClustersPerPlane[fNMaxPlanes] = {0};
665 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
666 chi2HistoryForBestCandidate[iPlane] = -1.;
667 nClustersPerPlane[iPlane] = fMFTClusterArray[iPlane] -> GetEntries();
670 fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks));
672 Int_t nClustersMC = 0;
673 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) nClustersMC += fIsGoodClusterInPlane[iPlane];
675 for (Int_t iTrack=0; iTrack<nFinalTracks; iTrack++) {
677 AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
679 Double_t chi2AtPlane[fNMaxPlanes] = {0};
680 Int_t nGoodClusters = 0;
681 Int_t nMFTClusters = finalTrack->GetNMFTClusters();
682 // Int_t nMUONClusters = finalTrack->GetNMUONClusters();
685 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
686 while (!finalTrack->PlaneExists(plane)) plane++;
687 AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
688 chi2AtPlane[plane] = localCluster->GetLocalChi2();
689 if (IsCorrectMatch(localCluster)) nGoodClusters++;
690 // Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope
691 // Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
692 // chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
696 if (fIsCurrentMuonTrackable) {
697 // fOutputQAFile->cd();
698 fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
701 // fOutputQAFile->cd();
703 Float_t finalCandidatesInfo[] = {Double_t(fRun),
705 Double_t(fCountRealTracksAnalyzedOfEvent),
706 Double_t(nFinalTracks),
707 Double_t(fLabelMC>=0),
710 Double_t(fMuonTrackReco->GetMatchTrigger()),
711 Double_t(nClustersMC),
712 Double_t(nGoodClusters),
724 fNtuFinalCandidates -> Fill(finalCandidatesInfo);
726 // now comparing the tracks with various criteria, in order to find the best one
728 Double_t theVariable = 0.;
729 // theVariable = chi2AtPlane[0];
730 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) theVariable += chi2AtPlane[iCluster];
731 theVariable /= Double_t(nMFTClusters);
734 if (theVariable<theVariable_Best || theVariable_Best<0.) {
735 nGoodClustersBestCandidate = nGoodClusters;
736 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = chi2AtPlane[iPlane];
737 theVariable_Best = theVariable;
738 fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
739 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
740 fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
741 0.90*fRPlaneMax[fNPlanesMFT-1],
742 Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
744 idBestCandidate = iTrack;
745 bestCandidateExists=kTRUE;
748 // ----------------------------------------------------------
753 fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
754 AliInfo(Form("fFinalBestCandidate->GetNMFTClusters() = %d\n", fFinalBestCandidate->GetNMFTClusters()));
755 PrintParticleHistory();
756 FillPlanesWithTrackHistory();
757 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
758 newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
759 newTrack -> SetTrackMCId(fRun*100000+fEv*1000+fCountRealTracksAnalyzedOfEvent);
760 new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
763 // fOutputQAFile->cd();
765 Float_t finalBestCandidatesInfo[] = {Double_t(fRun),
767 Double_t(fCountRealTracksAnalyzedOfEvent),
768 Double_t(nFinalTracks),
769 Double_t(fLabelMC>=0),
772 Double_t(fMuonTrackReco->GetMatchTrigger()),
773 Double_t(nClustersMC),
774 Double_t(nGoodClustersBestCandidate),
776 chi2HistoryForBestCandidate[0],
777 chi2HistoryForBestCandidate[1],
778 chi2HistoryForBestCandidate[2],
779 chi2HistoryForBestCandidate[3],
780 chi2HistoryForBestCandidate[4],
781 chi2HistoryForBestCandidate[5],
782 chi2HistoryForBestCandidate[6],
783 chi2HistoryForBestCandidate[7],
784 chi2HistoryForBestCandidate[8],
785 nClustersPerPlane[0],
786 nClustersPerPlane[1],
787 nClustersPerPlane[2],
788 nClustersPerPlane[3],
789 nClustersPerPlane[4],
790 nClustersPerPlane[5],
791 nClustersPerPlane[6],
792 nClustersPerPlane[7],
793 nClustersPerPlane[8]};
795 fNtuFinalBestCandidates -> Fill(finalBestCandidatesInfo);
797 if (fDrawOption && bestCandidateExists) {
798 fTxtTrackGoodClusters = new TLatex(0.20, 0.51, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
802 // -------------------------------------------------------------------------------------------
804 fCandidateTracks->Delete();
805 fFinalBestCandidate = NULL;
811 //===========================================================================================================================================
813 Int_t AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
815 AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId));
817 // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
819 // propagate track to plane #planeId (both to front and back active sensors)
820 // look for compatible clusters
821 // update TrackParam at found cluster (if any) using Kalman Filter
823 AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
825 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
826 currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
827 currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
828 currentParamForResearchFront = currentParamFront;
829 currentParamForResearchBack = currentParamBack;
830 Double_t xExtrap = gRandom->Gaus(0,fVertexErrorX);
831 Double_t yExtrap = gRandom->Gaus(0,fVertexErrorY);
832 Double_t zExtrap = gRandom->Gaus(0,fVertexErrorZ);
833 if (fBransonCorrection) {
834 AliMUONTrackExtrap::ExtrapToVertex(¤tParamFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
835 AliMUONTrackExtrap::ExtrapToVertex(¤tParamBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
838 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamFront, zExtrap);
839 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamBack, zExtrap);
841 AliMUONTrackExtrap::ExtrapToVertex(¤tParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
842 AliMUONTrackExtrap::ExtrapToVertex(¤tParamForResearchBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
844 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
845 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
846 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
847 currentParamForResearchFront = currentParamFront;
848 currentParamForResearchBack = currentParamBack;
849 AliMUONTrackExtrap::AddMCSEffect(¤tParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
850 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
851 AliMUONTrackExtrap::AddMCSEffect(¤tParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
852 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
853 AliMUONTrackExtrap::AddMCSEffect(¤tParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
854 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
855 AliMUONTrackExtrap::AddMCSEffect(¤tParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
856 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
858 // for all planes: extrapolation to the Z of the plane
859 AliMUONTrackExtrap::ExtrapToZCov(¤tParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
860 AliMUONTrackExtrap::ExtrapToZCov(¤tParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
861 AliMUONTrackExtrap::ExtrapToZCov(¤tParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
862 AliMUONTrackExtrap::ExtrapToZCov(¤tParamForResearchBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
864 //---------------------------------------------------------------------------------------
866 TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
867 TMatrixD covBack(5,5); covBack = currentParamForResearchBack.GetCovariances();
869 Double_t squaredError_X_Front = covFront(0,0);
870 Double_t squaredError_Y_Front = covFront(2,2);
871 Double_t squaredError_X_Back = covBack(0,0);
872 Double_t squaredError_Y_Back = covBack(2,2);
874 Double_t corrFact = 1.0;
876 Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
877 Double_t researchRadiusBack = TMath::Sqrt(squaredError_X_Back + squaredError_Y_Back);
878 if (0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtPlane[planeId]) {
879 corrFact = fMinResearchRadiusAtPlane[planeId]/(0.5*(researchRadiusFront+researchRadiusBack));
881 if (fIsCurrentMuonTrackable) {
882 // fOutputQAFile->cd();
883 fHistResearchRadius[planeId] -> Fill(0.5*(researchRadiusFront+researchRadiusBack));
886 Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor();
887 Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor();
888 Double_t position_X_Back = currentParamForResearchBack.GetNonBendingCoor();
889 Double_t position_Y_Back = currentParamForResearchBack.GetBendingCoor();
890 Double_t radialPositionOfTrackFront = TMath::Sqrt(position_X_Front*position_X_Front + position_Y_Front*position_Y_Front);
891 Double_t radialPositionOfTrackBack = TMath::Sqrt(position_X_Back*position_X_Back + position_Y_Back*position_Y_Back);
893 //---------------------------------------------------------------------------------------
895 Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut; // depends on the number of variables (here, 2)
897 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
899 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
900 AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
902 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
904 Bool_t isGoodChi2 = kFALSE;
906 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster);
907 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster); // describes the compatibility between the track and the cluster
908 if (chi2<chi2cut) isGoodChi2 = kTRUE;
910 Double_t radialPositionOfClusterFront = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
911 if (planeId == fNPlanesMFT-1) {
912 if (TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront)<fDistanceFromBestClusterAndTrackAtLastPlane ||
913 fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
914 fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
916 if (IsCorrectMatch(cluster)) {
917 fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
921 if (fIsCurrentMuonTrackable) {
922 // fOutputQAFile->cd();
923 if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf
924 else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf
928 AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
929 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
930 if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return 2;
931 newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
932 AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
933 planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
934 newTrack->SetPlaneExists(planeId);
935 AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
936 if (fIsCurrentMuonTrackable) {
937 Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
938 AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
939 Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope
940 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
941 // fOutputQAFile->cd();
942 if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
943 else fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
945 fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
947 else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
951 // Analyizing the clusters: BACK ACTIVE ELEMENTS
953 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
954 AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
956 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
958 Bool_t isGoodChi2 = kFALSE;
960 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster);
961 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster); // describes the compatibility between the track and the cluster
962 if (chi2<chi2cut) isGoodChi2 = kTRUE;
964 Double_t radialPositionOfClusterBack = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
965 if (planeId == fNPlanesMFT-1) {
966 if (TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack)<fDistanceFromBestClusterAndTrackAtLastPlane ||
967 fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
968 fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
970 if (IsCorrectMatch(cluster)) {
971 fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
975 if (fIsCurrentMuonTrackable) {
976 // fOutputQAFile->cd();
977 if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf
978 else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf
982 AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
983 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
984 if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return 2;
985 newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
986 AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
987 planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
988 newTrack->SetPlaneExists(planeId);
989 AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
990 if (fIsCurrentMuonTrackable) {
991 Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
992 AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
993 Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope
994 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
995 // fOutputQAFile->cd();
996 if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
997 else fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
999 fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
1001 else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
1005 //---------------------------------------------------------------------------------------------
1007 if (planeId == fNPlanesMFT-1) {
1008 if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) {
1009 // fOutputQAFile->cd();
1010 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane-
1011 fDistanceFromGoodClusterAndTrackAtLastPlane));
1012 fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane);
1020 //==========================================================================================================================================
1022 void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
1024 AliDebug(1, Form(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId));
1026 AliMUONTrackParam currentParamFront, currentParamBack;
1028 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
1029 currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
1030 currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
1031 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamFront, 0.);
1032 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamBack, 0.);
1034 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
1035 AliDebug(2, Form("fCurrentTrack = %p\n", fCurrentTrack));
1036 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
1037 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
1038 AliMUONTrackExtrap::AddMCSEffect(¤tParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
1039 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
1040 AliMUONTrackExtrap::AddMCSEffect(¤tParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
1041 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
1043 // for all planes: linear extrapolation to the Z of the plane
1044 AliMUONTrackExtrap::ExtrapToZCov(¤tParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
1045 AliMUONTrackExtrap::ExtrapToZCov(¤tParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
1047 Bool_t goodClusterFound = kFALSE;
1049 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
1051 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
1053 AliDebug(1, Form("nClustersFront = %d\n", nClustersFront));
1054 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
1055 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->UncheckedAt(iCluster);
1056 AliDebug(2, Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack));
1057 if (IsCorrectMatch(cluster)) {
1058 fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
1059 fCurrentTrack->SetPlaneExists(planeId);
1060 goodClusterFound = kTRUE;
1065 if (goodClusterFound) return;
1067 // Analyizing the clusters: BACK ACTIVE ELEMENTS
1069 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
1071 AliDebug(1, Form("nClustersBack = %d\n", nClustersBack));
1072 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
1073 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->UncheckedAt(iCluster);
1074 AliDebug(2,Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack));
1075 if (IsCorrectMatch(cluster)) {
1076 fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
1077 fCurrentTrack->SetPlaneExists(planeId);
1078 goodClusterFound = kTRUE;
1085 //==========================================================================================================================================
1087 void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() {
1089 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1090 fIsGoodClusterInPlane[iPlane] = kFALSE;
1091 Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
1092 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1093 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
1094 for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
1095 if (cluster->GetMCLabel(iTrack)==fLabelMC) {
1096 fIsGoodClusterInPlane[iPlane] = kTRUE;
1103 fIsCurrentMuonTrackable = kTRUE;
1104 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) fIsCurrentMuonTrackable = (fIsCurrentMuonTrackable&&fIsGoodClusterInPlane[iPlane]);
1108 //==========================================================================================================================================
1110 void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory() {
1112 // Fill planes with the clusters
1115 AliDebug(2, Form("fFinalBestCandidate->GetNMFTClusters() = %d\n", fFinalBestCandidate->GetNMFTClusters()));
1116 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1117 if (fFinalBestCandidate->PlaneExists(iPlane)) {
1118 AliMFTCluster *trackCluster = fFinalBestCandidate->GetMFTCluster(cluster++);
1119 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetPoint(fGrMFTPlane[kClusterOfTrack][iPlane]->GetN(), trackCluster->GetX(), trackCluster->GetY());
1121 Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
1122 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1123 AliMFTCluster *myCluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster);
1124 fGrMFTPlane[kAllClusters][iPlane] -> SetPoint(fGrMFTPlane[kAllClusters][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
1125 if (IsCorrectMatch(myCluster)) {
1126 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetPoint(fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
1133 //======================================================================================================================================
1135 Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) {
1137 Bool_t result = kFALSE;
1139 // check if the cluster belongs to the correct MC track
1141 for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
1142 if (cluster->GetMCLabel(iTrack)==fLabelMC) {
1148 AliDebug(2,Form("returning %d\n", result));
1154 //======================================================================================================================================
1156 Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
1158 // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
1159 // return the corresponding Chi2
1160 // assume the track parameters are given at the Z of the cluster
1162 // Set differences between trackParam and cluster in the bending and non bending directions
1163 Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
1164 Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
1165 AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
1167 // Calculate errors and covariances
1168 const TMatrixD& kParamCov = trackParam.GetCovariances();
1169 Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
1170 Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
1171 AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
1172 Double_t covXY = kParamCov(0,2);
1173 Double_t det = sigmaX2 * sigmaY2 - covXY * covXY;
1176 if (det==0.) return 1.e10;
1177 return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
1181 //=========================================================================================================================================
1183 void AliMuonForwardTrackFinder::SeparateFrontBackClusters() {
1185 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1186 fMFTClusterArrayFront[iPlane]->Delete();
1187 fMFTClusterArrayBack[iPlane] ->Delete();
1188 for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
1189 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
1190 if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
1191 new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
1194 new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
1201 //=========================================================================================================================================
1203 Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) {
1205 // the same definition as in AliMUONTrack is implemented, since here we just add more clusters to the Muon track
1207 Int_t ndf = 2 * nClusters - 5;
1208 return (ndf > 0) ? ndf : 0;
1212 //============================================================================================================================================
1214 void AliMuonForwardTrackFinder::BookHistos() {
1216 const Int_t nMaxNewTracks[] = {150, 200, 250, 600, 1000};
1217 const Double_t radiusPlane[] = {0.010, 0.010, 0.050, 0.5, 1.5};
1219 fHistRadiusEndOfAbsorber = new TH1D("hRadiusEndOfAbsorber", "Track radial distance at the end of the absorber", 1000, 0, 100.);
1221 fHistNGoodClustersForFinalTracks = new TH1D("hNGoodClustersForFinalTracks", "Number of Good Clusters per Final Track", 20, -0.25, 9.75);
1223 fHistDistanceGoodClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackAtLastPlane",
1224 "Distance of MC Good Cluster from Track in last MFT plane", 200, 0., 2.);
1226 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane =
1227 new TH1D("hDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane",
1228 "Good Cluster distance from track - Best Cluster distance from track in last MFT plane", 200, 0., 2.);
1230 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1232 fHistNTracksAfterExtrapolation[iPlane] = new TH1D(Form("hNTracksAfterExtrapolation_pl%02d", iPlane),
1233 Form("Number of Candidates after analysis of MFT plane %02d", iPlane),
1234 nMaxNewTracks[iPlane], -0.5, nMaxNewTracks[iPlane]-0.5);
1236 fHistResearchRadius[iPlane] = new TH1D(Form("hResearchRadius_pl%02d", iPlane),
1237 Form("Research Radius for candidate clusters in MFT plane %02d", iPlane),
1238 1000, 0., radiusPlane[iPlane]);
1240 fHistChi2Cluster_GoodCluster[iPlane] = new TH1D(Form("hChi2Cluster_GoodCluster_pl%02d", iPlane),
1241 Form("#chi^{2}_{clust} for Good clusters in MFT plane %02d", iPlane),
1244 fHistChi2Cluster_BadCluster[iPlane] = new TH1D(Form("hChi2Cluster_BadCluster_pl%02d", iPlane),
1245 Form("#chi^{2}_{clust} for Bad clusters in MFT plane %02d", iPlane),
1248 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1249 Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
1252 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1253 Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
1258 //------------------------------------------
1260 fHistRadiusEndOfAbsorber -> Sumw2();
1261 fHistNGoodClustersForFinalTracks -> Sumw2();
1263 fHistDistanceGoodClusterFromTrackAtLastPlane -> Sumw2();
1264 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Sumw2();
1266 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1268 fHistNTracksAfterExtrapolation[iPlane] -> Sumw2();
1269 fHistResearchRadius[iPlane] -> Sumw2();
1271 fHistChi2Cluster_GoodCluster[iPlane] -> Sumw2();
1272 fHistChi2Cluster_BadCluster[iPlane] -> Sumw2();
1274 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1275 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1279 fNtuFinalCandidates = new TNtuple("ntuFinalCandidates", "Final Candidates (ALL)", "run:event:muonTrack:nFinalCandidates:MCTrackRefExists:xVtx:yVtx:zVtx:motherPdg:triggerMatch:nClustersMC:nGoodClusters:pt:theta:eta:chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8");
1281 fNtuFinalBestCandidates = new TNtuple("ntuFinalBestCandidates", "Final Best Candidates", "run:event:muonTrack:nFinalCandidates:MCTrackRefExists:xVtx:yVtx:zVtx:motherPdg:triggerMatch:nClustersMC:nGoodClusters:pt:theta:eta:chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:nClustersAtPlane0:nClustersAtPlane1:nClustersAtPlane2:nClustersAtPlane3:nClustersAtPlane4:nClustersAtPlane5:nClustersAtPlane6:nClustersAtPlane7:nClustersAtPlane8");
1285 //============================================================================================================================================
1287 void AliMuonForwardTrackFinder::SetTitleHistos() {
1289 fHistRadiusEndOfAbsorber -> SetXTitle("R_{abs} [cm]");
1290 fHistNGoodClustersForFinalTracks -> SetXTitle("N_{GoodClusters}");
1292 fHistDistanceGoodClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]");
1293 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]");
1296 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1298 fHistNTracksAfterExtrapolation[iPlane] -> SetXTitle("N_{tracks}");
1299 fHistResearchRadius[iPlane] -> SetXTitle("Research Radius [cm]");
1301 fHistChi2Cluster_GoodCluster[iPlane] -> SetXTitle("#chi^{2}/ndf");
1302 fHistChi2Cluster_BadCluster[iPlane] -> SetXTitle("#chi^{2}/ndf");
1304 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1305 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1311 //===========================================================================================================================================
1313 void AliMuonForwardTrackFinder::BookPlanes() {
1315 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1316 fGrMFTPlane[kAllClusters][iPlane] = new TGraph();
1317 fGrMFTPlane[kAllClusters][iPlane] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
1318 fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerStyle(20);
1319 // fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.5);
1320 // fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.3);
1321 fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.2);
1324 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1325 fGrMFTPlane[kClustersGoodChi2][iPlane] = new TGraph();
1326 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
1327 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerStyle(20);
1328 // fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.8);
1329 // fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.4);
1330 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.3);
1331 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerColor(kBlue);
1334 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1335 fGrMFTPlane[kClusterOfTrack][iPlane] = new TGraph();
1336 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
1337 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerStyle(25);
1338 // fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(1.2);
1339 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(0.9);
1340 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerColor(kRed);
1341 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
1344 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1345 fGrMFTPlane[kClusterCorrectMC][iPlane] = new TGraph();
1346 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
1347 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerStyle(20);
1348 // fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.8);
1349 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.5);
1350 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerColor(kGreen);
1353 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1354 fCircleExt[iPlane] = new TEllipse(0., 0., fRPlaneMax[iPlane], fRPlaneMax[iPlane]);
1355 fCircleInt[iPlane] = new TEllipse(0., 0., fRPlaneMin[iPlane], fRPlaneMin[iPlane]);
1358 fTxtDummy = new TLatex(0.10, 0.59, "Best Candidate:");
1360 //---------------------------------------------------
1362 fMrkAllClust = new TMarker(0.10, 0.32, 20);
1363 fMrkAllClust -> SetMarkerSize(0.5);
1365 fMrkClustGoodChi2 = new TMarker(0.10, 0.26, 20);
1366 fMrkClustGoodChi2 -> SetMarkerSize(0.8);
1367 fMrkClustGoodChi2 -> SetMarkerColor(kBlue);
1369 fMrkClustMC = new TMarker(0.10, 0.20, 20);
1370 fMrkClustMC -> SetMarkerSize(0.8);
1371 fMrkClustMC -> SetMarkerColor(kGreen);
1373 fMrkClustOfTrack = new TMarker(0.10, 0.14, 25);
1374 fMrkClustOfTrack -> SetMarkerSize(1.2);
1375 fMrkClustOfTrack -> SetMarkerColor(kRed);
1377 fTxtAllClust = new TLatex(0.15, 0.30, "All Clusters");
1378 fTxtAllClust -> SetTextSize(0.040);
1380 fTxtClustGoodChi2 = new TLatex(0.15, 0.24, "Clusters involved in the research");
1381 fTxtClustGoodChi2 -> SetTextSize(0.040);
1383 fTxtClustMC = new TLatex(0.15, 0.18, "MC good clusters");
1384 fTxtClustMC -> SetTextSize(0.040);
1386 fTxtClustOfTrack = new TLatex(0.15, 0.12, "Clusters of the best candidate");
1387 fTxtClustOfTrack -> SetTextSize(0.040);
1391 //===========================================================================================================================================
1393 void AliMuonForwardTrackFinder::ResetPlanes() {
1395 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1396 for (Int_t iGr=0; iGr<4; iGr++) {
1397 Int_t nOldClusters = fGrMFTPlane[iGr][iPlane]->GetN();
1398 for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iGr][iPlane]->RemovePoint(iPoint);
1404 //===========================================================================================================================================
1406 void AliMuonForwardTrackFinder::PrintParticleHistory() {
1408 AliDebug(1, "Entering");
1410 TString history = "";
1412 TParticle *part = 0;
1413 if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
1415 AliDebug(1, Form("fStack->Particle(%d) = %p", fLabelMC, part));
1416 AliDebug(1, Form("fStack->Particle(%d)->GetPdgCode() = %d", fLabelMC, part->GetPdgCode()));
1419 if (part->GetFirstMother() != -1) {
1420 TParticle *partMother = fStack->Particle(part->GetFirstMother());
1421 AliDebug(1, Form("fStack->Particle(%d) = %p", part->GetFirstMother(), partMother));
1423 Char_t newName[100];
1424 if (partMother->GetFirstMother() != -1) history += "... #rightarrow ";
1425 PDGNameConverter(partMother->GetName(), newName);
1426 history += Form("%s #rightarrow ", newName);
1429 Char_t newName[100];
1430 AliDebug(1, Form("fStack->Particle(%d)->GetPdgCode() = %d", fLabelMC, part->GetPdgCode()));
1431 AliDebug(1, Form("fStack->Particle(%d)->GetName() = %s", fLabelMC, part->GetName()));
1432 PDGNameConverter(part->GetName(), newName);
1433 history += Form("%s at z = %5.1f cm", newName, part->Vz());
1434 // printf("%s", history.Data());
1436 else history += "NO AVAILABLE HISTORY";
1438 fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data());
1440 // Filling particle history in the fFinalBestCandidate
1443 for (Int_t iParent=0; iParent<AliMuonForwardTrack::fgkNParentsMax; iParent++) {
1444 if (part->GetFirstMother() == -1) break;
1445 if (!(fStack->Particle(part->GetFirstMother()))) break;
1446 AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", fStack->Particle(part->GetFirstMother())));
1447 fFinalBestCandidate->SetParentMCLabel(iParent, part->GetFirstMother());
1448 fFinalBestCandidate->SetParentPDGCode(iParent, fStack->Particle(part->GetFirstMother())->GetPdgCode());
1449 part = fStack->Particle(part->GetFirstMother());
1455 //===========================================================================================================================================
1457 Bool_t AliMuonForwardTrackFinder::IsMother(const Char_t *nameMother) {
1459 Bool_t result = kFALSE;
1461 TParticle *part = 0;
1462 if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
1465 if (part->GetFirstMother() != -1) {
1466 TParticle *partMother = fStack->Particle(part->GetFirstMother());
1468 if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
1477 //===========================================================================================================================================
1479 void AliMuonForwardTrackFinder::DrawPlanes() {
1482 if (fNPlanesMFT <= 5) fCanvas -> Divide(3,2);
1483 else if (fNPlanesMFT <= 11) fCanvas -> Divide(4,3);
1484 else if (fNPlanesMFT <= 19) fCanvas -> Divide(5,4);
1486 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1488 fCanvas->cd(fNPlanesMFT-iPlane+1);
1490 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1491 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1492 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetTitle("X [cm]");
1493 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetTitle("Y [cm]");
1494 fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("ap");
1496 fCircleExt[iPlane] -> Draw("same");
1497 fCircleInt[iPlane] -> Draw("same");
1499 if (fGrMFTPlane[kAllClusters][iPlane]->GetN()) fGrMFTPlane[kAllClusters][iPlane] -> Draw("psame");
1500 if (fGrMFTPlane[kClustersGoodChi2][iPlane]->GetN()) fGrMFTPlane[kClustersGoodChi2][iPlane] -> Draw("psame");
1501 if (fGrMFTPlane[kClusterOfTrack][iPlane]->GetN()) fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("psame");
1502 if (fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN()) fGrMFTPlane[kClusterCorrectMC][iPlane] -> Draw("psame");
1504 fTxtTrackChi2[iPlane] -> Draw("same");
1509 fTxtMuonHistory -> Draw();
1510 fTxtDummy -> Draw("same");
1511 if (fMatchingMode==kRealMatching) fTxtTrackGoodClusters -> Draw("same");
1512 fTxtTrackFinalChi2 -> Draw("same");
1513 fTxtTrackMomentum -> Draw("same");
1514 if (fMatchingMode==kRealMatching) fTxtFinalCandidates -> Draw("same");
1516 fMrkAllClust -> Draw("same");
1517 fMrkClustGoodChi2 -> Draw("same");
1518 fMrkClustMC -> Draw("same");
1519 fMrkClustOfTrack -> Draw("same");
1521 fTxtAllClust -> Draw("same");
1522 fTxtClustGoodChi2 -> Draw("same");
1523 fTxtClustMC -> Draw("same");
1524 fTxtClustOfTrack -> Draw("same");
1526 // fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1527 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1528 if (IsMother("phi")) {
1529 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1530 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1532 if (IsMother("J/psi")) {
1533 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1534 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1539 //===========================================================================================================================================
1541 void AliMuonForwardTrackFinder::Terminate() {
1544 AliInfo("---------------------------------------------------------------------------------------------------------------");
1545 AliInfo(Form("%8d tracks analyzed", fCountRealTracksAnalyzed));
1546 AliInfo(Form("%8d tracks with MC ref", fCountRealTracksWithRefMC));
1547 AliInfo(Form("%8d tracks with MC ref & trigger match", fCountRealTracksWithRefMC_andTrigger));
1548 if (fMatchingMode==kRealMatching) {
1549 AliInfo(Form("%8d tracks analyzed with final candidates", fCountRealTracksAnalyzedWithFinalCandidates));
1552 AliInfo(Form("%8d tracks matched with their MC clusters", fCountRealTracksAnalyzedWithFinalCandidates));
1554 // printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c", fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
1555 // printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
1556 AliInfo("---------------------------------------------------------------------------------------------------------------");
1563 //==========================================================================================================================================
1565 void AliMuonForwardTrackFinder::FillOutputTree() {
1567 if (!fMuonForwardTracks || !fOutputEventTree) return;
1569 AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p",
1570 fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0)));
1572 // fOutputTreeFile->cd();
1573 fOutputEventTree->Fill();
1574 AliDebug(1, Form("\nFilled Tree: nEvents = %d!!!!\n", Int_t(fOutputEventTree->GetEntries())));
1578 //==========================================================================================================================================
1580 void AliMuonForwardTrackFinder::WriteOutputTree() {
1582 if (!fOutputEventTree || !fOutputTreeFile) return;
1584 fOutputTreeFile -> cd();
1586 fOutputEventTree -> Write();
1587 fOutputTreeFile -> Close();
1591 //==========================================================================================================================================
1593 void AliMuonForwardTrackFinder::WriteHistos() {
1595 fOutputQAFile = new TFile(Form("MuonGlobalTracking.QA.run%d.root", fRun), "recreate");
1596 fOutputQAFile -> cd();
1598 fHistRadiusEndOfAbsorber -> Write();
1599 fHistNGoodClustersForFinalTracks -> Write();
1601 fHistDistanceGoodClusterFromTrackAtLastPlane -> Write();
1602 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Write();
1604 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1606 fHistNTracksAfterExtrapolation[iPlane] -> Write();
1607 fHistResearchRadius[iPlane] -> Write();
1609 fHistChi2Cluster_GoodCluster[iPlane] -> Write();
1610 fHistChi2Cluster_BadCluster[iPlane] -> Write();
1612 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
1613 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Write();
1617 fNtuFinalCandidates -> Write();
1618 fNtuFinalBestCandidates -> Write();
1620 fOutputQAFile -> Close();
1624 //===========================================================================================================================================
1626 void AliMuonForwardTrackFinder::PDGNameConverter(const Char_t *nameIn, Char_t *nameOut) {
1628 if (!strcmp(nameIn, "mu+")) snprintf(nameOut, 50, "#mu^{+}");
1629 else if (!strcmp(nameIn, "mu-")) snprintf(nameOut, 50, "#mu^{-}");
1630 else if (!strcmp(nameIn, "pi+")) snprintf(nameOut, 50, "#pi^{+}");
1631 else if (!strcmp(nameIn, "pi-")) snprintf(nameOut, 50, "#pi^{-}");
1632 else if (!strcmp(nameIn, "K+")) snprintf(nameOut, 50, "K^{+}");
1633 else if (!strcmp(nameIn, "K-")) snprintf(nameOut, 50, "K^{-}");
1634 else if (!strcmp(nameIn, "K*+")) snprintf(nameOut, 50, "K^{*+}");
1635 else if (!strcmp(nameIn, "K*-")) snprintf(nameOut, 50, "K^{*-}");
1636 else if (!strcmp(nameIn, "K_S0")) snprintf(nameOut, 50, "K_{S}^{0}");
1637 else if (!strcmp(nameIn, "K_L0")) snprintf(nameOut, 50, "K_{L}^{0}");
1638 else if (!strcmp(nameIn, "K0")) snprintf(nameOut, 50, "K^{0}");
1639 else if (!strcmp(nameIn, "K0_bar")) snprintf(nameOut, 50, "#bar{K}^{0}");
1640 else if (!strcmp(nameIn, "K*0")) snprintf(nameOut, 50, "K^{*0}");
1641 else if (!strcmp(nameIn, "K*0_bar")) snprintf(nameOut, 50, "#bar{K}^{*0}");
1642 else if (!strcmp(nameIn, "rho0")) snprintf(nameOut, 50, "#rho^{0}");
1643 else if (!strcmp(nameIn, "rho+")) snprintf(nameOut, 50, "#rho^{+}");
1644 else if (!strcmp(nameIn, "rho-")) snprintf(nameOut, 50, "#rho^{-}");
1645 else if (!strcmp(nameIn, "omega")) snprintf(nameOut, 50, "#omega");
1646 else if (!strcmp(nameIn, "eta'")) snprintf(nameOut, 50, "#eta'");
1647 else if (!strcmp(nameIn, "phi")) snprintf(nameOut, 50, "#phi");
1649 else if (!strcmp(nameIn, "D-")) snprintf(nameOut, 50, "D^{-}");
1650 else if (!strcmp(nameIn, "D+")) snprintf(nameOut, 50, "D^{+}");
1651 else if (!strcmp(nameIn, "D0")) snprintf(nameOut, 50, "D^{0}");
1652 else if (!strcmp(nameIn, "D0_bar")) snprintf(nameOut, 50, "#bar{D}^{0}");
1653 else if (!strcmp(nameIn, "D*-")) snprintf(nameOut, 50, "D^{*-}");
1654 else if (!strcmp(nameIn, "D*+")) snprintf(nameOut, 50, "D^{*+}");
1655 else if (!strcmp(nameIn, "D_s+")) snprintf(nameOut, 50, "D_{s}^{+}");
1656 else if (!strcmp(nameIn, "D*_s+")) snprintf(nameOut, 50, "D_{s}^{*+}");
1658 else if (!strcmp(nameIn, "B-")) snprintf(nameOut, 50, "B^{-}");
1659 else if (!strcmp(nameIn, "B+")) snprintf(nameOut, 50, "B^{+}");
1660 else if (!strcmp(nameIn, "B_s0_bar")) snprintf(nameOut, 50, "#bar{B}_{s}^{0}");
1662 else if (!strcmp(nameIn, "antiproton")) snprintf(nameOut, 50, "#bar{p}");
1663 else if (!strcmp(nameIn, "proton")) snprintf(nameOut, 50, "p");
1664 else if (!strcmp(nameIn, "neutron")) snprintf(nameOut, 50, "n");
1665 else if (!strcmp(nameIn, "Sigma+")) snprintf(nameOut, 50, "#Sigma^{+}");
1666 else if (!strcmp(nameIn, "Delta+")) snprintf(nameOut, 50, "#Delta{+}");
1667 else if (!strcmp(nameIn, "Delta--")) snprintf(nameOut, 50, "#Delta{--}");
1668 else if (!strcmp(nameIn, "Lambda0")) snprintf(nameOut, 50, "#Lambda_0");
1669 else if (!strcmp(nameIn, "Lambda0_bar")) snprintf(nameOut, 50, "#bar{Lambda}_0");
1671 else snprintf(nameOut, 50, "%s", nameIn);
1675 //===========================================================================================================================================
1677 void AliMuonForwardTrackFinder::SetDraw(Bool_t drawOption) {
1679 fDrawOption = drawOption;
1682 fCanvas = new TCanvas("tracking", "tracking", 1200, 800);
1683 fCanvas -> Divide(3,2);
1688 //===========================================================================================================================================
1690 Bool_t AliMuonForwardTrackFinder::InitGRP() {
1692 //------------------------------------
1693 // Initialization of the GRP entry
1694 //------------------------------------
1696 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
1700 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
1703 AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
1705 fGRPData = new AliGRPObject();
1706 fGRPData->ReadValuesFromMap(m);
1710 AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
1711 fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
1715 // FIX ME: The unloading of GRP entry is temporarily disabled
1716 // because ZDC and VZERO are using it in order to initialize
1717 // their reconstructor objects. In the future one has to think
1718 // of propagating AliRunInfo to the reconstructors.
1719 // AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
1723 AliError("No GRP entry found in OCDB!");
1727 TString lhcState = fGRPData->GetLHCState();
1728 if (lhcState==AliGRPObject::GetInvalidString()) {
1729 AliError("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN");
1730 lhcState = "UNKNOWN";
1733 TString beamType = fGRPData->GetBeamType();
1734 if (beamType==AliGRPObject::GetInvalidString()) {
1735 AliError("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN");
1736 beamType = "UNKNOWN";
1739 Float_t beamEnergy = fGRPData->GetBeamEnergy();
1740 if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
1741 AliError("GRP/GRP/Data entry: missing value for the beam energy ! Using 0");
1745 TString runType = fGRPData->GetRunType();
1746 if (runType==AliGRPObject::GetInvalidString()) {
1747 AliError("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN");
1748 runType = "UNKNOWN";
1751 Int_t activeDetectors = fGRPData->GetDetectorMask();
1752 if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
1753 AliError("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399");
1754 activeDetectors = 1074790399;
1756 AliDebug(1, Form("activeDetectors = %d", activeDetectors));
1758 fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
1761 // *** Dealing with the magnetic field map
1763 if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
1764 if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
1765 AliInfo("ExpertMode!!! GRP information will be ignored !");
1766 AliInfo("ExpertMode!!! Running with the externally locked B field !");
1769 AliInfo("Destroying existing B field instance!");
1770 delete TGeoGlobalMagField::Instance();
1773 if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
1774 // Construct the field map out of the information retrieved from GRP.
1777 Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
1778 if (l3Current == AliGRPObject::GetInvalidFloat()) {
1779 AliError("GRP/GRP/Data entry: missing value for the L3 current !");
1783 Char_t l3Polarity = fGRPData->GetL3Polarity();
1784 if (l3Polarity == AliGRPObject::GetInvalidChar()) {
1785 AliError("GRP/GRP/Data entry: missing value for the L3 polarity !");
1790 Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
1791 if (diCurrent == AliGRPObject::GetInvalidFloat()) {
1792 AliError("GRP/GRP/Data entry: missing value for the dipole current !");
1796 Char_t diPolarity = fGRPData->GetDipolePolarity();
1797 if (diPolarity == AliGRPObject::GetInvalidChar()) {
1798 AliError("GRP/GRP/Data entry: missing value for the dipole polarity !");
1802 // read special bits for the polarity convention and map type
1803 Int_t polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
1804 Bool_t uniformB = fGRPData->IsUniformBMap();
1807 AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
1808 TMath::Abs(diCurrent) * (diPolarity ? -1:1),
1809 polConvention,uniformB,beamEnergy, beamType.Data());
1811 TGeoGlobalMagField::Instance()->SetField( fld );
1812 TGeoGlobalMagField::Instance()->Lock();
1813 AliInfo("Running with the B field constructed out of GRP !");
1815 else AliFatal("Failed to create a B field map !");
1817 else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1823 //====================================================================================================================================================
1825 Bool_t AliMuonForwardTrackFinder::SetRunNumber() {
1827 AliCDBManager *man = AliCDBManager::Instance();
1830 AliError("No run loader found!");
1834 fRunLoader->LoadHeader();
1835 // read run number from gAlice
1836 if (fRunLoader->GetHeader()) {
1837 man->SetRun(fRunLoader->GetHeader()->GetRun());
1838 fRunLoader->UnloadHeader();
1841 AliError("No run-loader header found!");
1850 //====================================================================================================================================================