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);
559 track -> SetRAtAbsorberEnd(rAbsorber);
561 // if the radial distance of the track at the end of the absorber is smaller than a given radius, skip to the next track
562 if (rAbsorber < fRAbsorberCut) return 4;
564 //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
566 for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) { // *** do not reverse the order of this cycle!!!
567 // *** this reflects the fact that the extrapolation is performed
568 // *** starting from the last MFT plane back to the origin
570 // --------- updating the array of tracks according to the clusters available in the i-th plane ---------
572 fNPlanesMFTAnalyzed++;
574 if (fMatchingMode==kRealMatching) {
575 Int_t nTracksToBeAnalyzed = fCandidateTracks->GetEntriesFast();
576 for (Int_t iTrack=0; iTrack<nTracksToBeAnalyzed; iTrack++) {
577 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
578 // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
579 // (several new tracks can be created for one old track)
580 if (FindClusterInPlane(iPlane) == kDiverged) {
581 fGlobalTrackingDiverged = kTRUE;
582 if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1;
585 if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
586 fCandidateTracks->Remove(fCurrentTrack); // the old track is removed after the check;
589 fCandidateTracks->Compress();
590 if (fIsCurrentMuonTrackable) {
591 // fOutputQAFile->cd();
592 fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
596 else if (fMatchingMode==kIdealMatching && fIsCurrentMuonTrackable) {
597 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
598 AliDebug(2, Form("plane %02d: fCandidateTracks->GetEntriesFast() = %d fCandidateTracks->UncheckedAt(0) = %p fCurrentTrack = %p\n",
599 iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack));
600 AttachGoodClusterInPlane(iPlane);
605 // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
607 fGlobalTrackingDiverged = kFALSE;
608 fScaleSigmaClusterCut = 1.0;
610 AliDebug(1, "Finished cycle over planes");
612 Double_t momentum = pt * TMath::CosH(eta);
613 fTxtTrackMomentum = new TLatex(0.10, 0.70, Form("P_{spectro} = %3.1f GeV/c", momentum));
615 if (fMatchingMode==kIdealMatching) {
616 AliDebug(1, "Adding track to output tree...\n");
617 fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
618 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
619 new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
620 AliDebug(1, "...track added!\n");
621 fCandidateTracks->Delete();
622 fCountRealTracksAnalyzedOfEvent++;
623 fCountRealTracksAnalyzedWithFinalCandidates++;
624 PrintParticleHistory();
625 FillPlanesWithTrackHistory();
627 Double_t chi2AtPlane[fNMaxPlanes] = {0};
628 Int_t nGoodClusters = 0;
629 Int_t nMFTClusters = fFinalBestCandidate->GetNMFTClusters();
630 // Int_t nMUONClusters = fFinalBestCandidate->GetNMUONClusters();
632 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
633 while (!fFinalBestCandidate->PlaneExists(plane)) plane++;
634 AliMFTCluster *localCluster = fFinalBestCandidate->GetMFTCluster(iCluster);
635 chi2AtPlane[plane] = localCluster->GetLocalChi2();
636 if (IsCorrectMatch(localCluster)) nGoodClusters++;
637 // Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope
638 // Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
639 // chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
642 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
643 fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
644 0.90*fRPlaneMax[fNPlanesMFT-1],
645 Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
647 fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2AtPlane[0]));
649 if (fDrawOption) DrawPlanes();
653 // If we have several final tracks, we must find the best candidate:
655 Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
656 AliDebug(1, Form("nFinalTracks = %d", nFinalTracks));
658 if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++;
660 Double_t theVariable_Best = -1.; // variable defining the best candidate
661 Bool_t bestCandidateExists = kFALSE;
662 Int_t nGoodClustersBestCandidate = 0;
663 Int_t idBestCandidate = 0;
664 Double_t chi2HistoryForBestCandidate[fNMaxPlanes] = {0}; // chi2 on each plane, for the best candidate
665 Double_t nClustersPerPlane[fNMaxPlanes] = {0};
666 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
667 chi2HistoryForBestCandidate[iPlane] = -1.;
668 nClustersPerPlane[iPlane] = fMFTClusterArray[iPlane] -> GetEntries();
671 fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks));
673 Int_t nClustersMC = 0;
674 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) nClustersMC += fIsGoodClusterInPlane[iPlane];
676 for (Int_t iTrack=0; iTrack<nFinalTracks; iTrack++) {
678 AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
680 Double_t chi2AtPlane[fNMaxPlanes] = {0};
681 Int_t nGoodClusters = 0;
682 Int_t nMFTClusters = finalTrack->GetNMFTClusters();
683 // Int_t nMUONClusters = finalTrack->GetNMUONClusters();
686 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
687 while (!finalTrack->PlaneExists(plane)) plane++;
688 AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
689 chi2AtPlane[plane] = localCluster->GetLocalChi2();
690 if (IsCorrectMatch(localCluster)) nGoodClusters++;
691 // Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope
692 // Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
693 // chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
697 if (fIsCurrentMuonTrackable) {
698 // fOutputQAFile->cd();
699 fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
702 // fOutputQAFile->cd();
704 Float_t finalCandidatesInfo[] = {static_cast<Float_t>(fRun),
705 static_cast<Float_t>(fEv),
706 static_cast<Float_t>(fCountRealTracksAnalyzedOfEvent),
707 static_cast<Float_t>(nFinalTracks),
708 static_cast<Float_t>(fLabelMC>=0),
709 static_cast<Float_t>(xVtx), static_cast<Float_t>(yVtx), static_cast<Float_t>(zVtx),
710 static_cast<Float_t>(motherPdg),
711 static_cast<Float_t>(fMuonTrackReco->GetMatchTrigger()),
712 static_cast<Float_t>(nClustersMC),
713 static_cast<Float_t>(nGoodClusters),
714 static_cast<Float_t>(pt), static_cast<Float_t>(theta), static_cast<Float_t>(eta),
715 static_cast<Float_t>(chi2AtPlane[0]),
716 static_cast<Float_t>(chi2AtPlane[1]),
717 static_cast<Float_t>(chi2AtPlane[2]),
718 static_cast<Float_t>(chi2AtPlane[3]),
719 static_cast<Float_t>(chi2AtPlane[4]),
720 static_cast<Float_t>(chi2AtPlane[5]),
721 static_cast<Float_t>(chi2AtPlane[6]),
722 static_cast<Float_t>(chi2AtPlane[7]),
723 static_cast<Float_t>(chi2AtPlane[8])};
725 fNtuFinalCandidates -> Fill(finalCandidatesInfo);
727 // now comparing the tracks with various criteria, in order to find the best one
729 Double_t theVariable = 0.;
730 // theVariable = chi2AtPlane[0];
731 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) theVariable += chi2AtPlane[iCluster];
732 theVariable /= Double_t(nMFTClusters);
735 if (theVariable<theVariable_Best || theVariable_Best<0.) {
736 nGoodClustersBestCandidate = nGoodClusters;
737 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = chi2AtPlane[iPlane];
738 theVariable_Best = theVariable;
739 fTxtTrackFinalChi2 = new TLatex(0.20, 0.44, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
740 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
741 fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
742 0.90*fRPlaneMax[fNPlanesMFT-1],
743 Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
745 idBestCandidate = iTrack;
746 bestCandidateExists=kTRUE;
749 // ----------------------------------------------------------
754 fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
755 AliInfo(Form("fFinalBestCandidate->GetNMFTClusters() = %d\n", fFinalBestCandidate->GetNMFTClusters()));
756 PrintParticleHistory();
757 FillPlanesWithTrackHistory();
758 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
759 newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate);
760 newTrack -> SetTrackMCId(fRun*100000+fEv*1000+fCountRealTracksAnalyzedOfEvent);
761 new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack);
764 // fOutputQAFile->cd();
766 Float_t finalBestCandidatesInfo[] = {static_cast<Float_t>(fRun),
767 static_cast<Float_t>(fEv),
768 static_cast<Float_t>(fCountRealTracksAnalyzedOfEvent),
769 static_cast<Float_t>(nFinalTracks),
770 static_cast<Float_t>(fLabelMC>=0),
771 static_cast<Float_t>(xVtx), static_cast<Float_t>(yVtx), static_cast<Float_t>(zVtx),
772 static_cast<Float_t>(motherPdg),
773 static_cast<Float_t>(fMuonTrackReco->GetMatchTrigger()),
774 static_cast<Float_t>(nClustersMC),
775 static_cast<Float_t>(nGoodClustersBestCandidate),
776 static_cast<Float_t>(pt), static_cast<Float_t>(theta), static_cast<Float_t>(eta),
777 static_cast<Float_t>(chi2HistoryForBestCandidate[0]),
778 static_cast<Float_t>(chi2HistoryForBestCandidate[1]),
779 static_cast<Float_t>(chi2HistoryForBestCandidate[2]),
780 static_cast<Float_t>(chi2HistoryForBestCandidate[3]),
781 static_cast<Float_t>(chi2HistoryForBestCandidate[4]),
782 static_cast<Float_t>(chi2HistoryForBestCandidate[5]),
783 static_cast<Float_t>(chi2HistoryForBestCandidate[6]),
784 static_cast<Float_t>(chi2HistoryForBestCandidate[7]),
785 static_cast<Float_t>(chi2HistoryForBestCandidate[8]),
786 static_cast<Float_t>(nClustersPerPlane[0]),
787 static_cast<Float_t>(nClustersPerPlane[1]),
788 static_cast<Float_t>(nClustersPerPlane[2]),
789 static_cast<Float_t>(nClustersPerPlane[3]),
790 static_cast<Float_t>(nClustersPerPlane[4]),
791 static_cast<Float_t>(nClustersPerPlane[5]),
792 static_cast<Float_t>(nClustersPerPlane[6]),
793 static_cast<Float_t>(nClustersPerPlane[7]),
794 static_cast<Float_t>(nClustersPerPlane[8])};
796 fNtuFinalBestCandidates -> Fill(finalBestCandidatesInfo);
798 if (fDrawOption && bestCandidateExists) {
799 fTxtTrackGoodClusters = new TLatex(0.20, 0.51, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
803 // -------------------------------------------------------------------------------------------
805 fCandidateTracks->Delete();
806 fFinalBestCandidate = NULL;
812 //===========================================================================================================================================
814 Int_t AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
816 AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId));
818 // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
820 // propagate track to plane #planeId (both to front and back active sensors)
821 // look for compatible clusters
822 // update TrackParam at found cluster (if any) using Kalman Filter
824 AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
826 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
827 currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
828 currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
829 currentParamForResearchFront = currentParamFront;
830 currentParamForResearchBack = currentParamBack;
831 Double_t xExtrap = gRandom->Gaus(0,fVertexErrorX);
832 Double_t yExtrap = gRandom->Gaus(0,fVertexErrorY);
833 Double_t zExtrap = gRandom->Gaus(0,fVertexErrorZ);
834 if (fBransonCorrection) {
835 AliMUONTrackExtrap::ExtrapToVertex(¤tParamFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
836 AliMUONTrackExtrap::ExtrapToVertex(¤tParamBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
839 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamFront, zExtrap);
840 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamBack, zExtrap);
842 AliMUONTrackExtrap::ExtrapToVertex(¤tParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
843 AliMUONTrackExtrap::ExtrapToVertex(¤tParamForResearchBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY);
845 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
846 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
847 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
848 currentParamForResearchFront = currentParamFront;
849 currentParamForResearchBack = currentParamBack;
850 AliMUONTrackExtrap::AddMCSEffect(¤tParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
851 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
852 AliMUONTrackExtrap::AddMCSEffect(¤tParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
853 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
854 AliMUONTrackExtrap::AddMCSEffect(¤tParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
855 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
856 AliMUONTrackExtrap::AddMCSEffect(¤tParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
857 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
859 // for all planes: extrapolation to the Z of the plane
860 AliMUONTrackExtrap::ExtrapToZCov(¤tParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
861 AliMUONTrackExtrap::ExtrapToZCov(¤tParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
862 AliMUONTrackExtrap::ExtrapToZCov(¤tParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
863 AliMUONTrackExtrap::ExtrapToZCov(¤tParamForResearchBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
865 //---------------------------------------------------------------------------------------
867 TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
868 TMatrixD covBack(5,5); covBack = currentParamForResearchBack.GetCovariances();
870 Double_t squaredError_X_Front = covFront(0,0);
871 Double_t squaredError_Y_Front = covFront(2,2);
872 Double_t squaredError_X_Back = covBack(0,0);
873 Double_t squaredError_Y_Back = covBack(2,2);
875 Double_t corrFact = 1.0;
877 Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
878 Double_t researchRadiusBack = TMath::Sqrt(squaredError_X_Back + squaredError_Y_Back);
879 if (0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtPlane[planeId]) {
880 corrFact = fMinResearchRadiusAtPlane[planeId]/(0.5*(researchRadiusFront+researchRadiusBack));
882 if (fIsCurrentMuonTrackable) {
883 // fOutputQAFile->cd();
884 fHistResearchRadius[planeId] -> Fill(0.5*(researchRadiusFront+researchRadiusBack));
887 Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor();
888 Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor();
889 Double_t position_X_Back = currentParamForResearchBack.GetNonBendingCoor();
890 Double_t position_Y_Back = currentParamForResearchBack.GetBendingCoor();
891 Double_t radialPositionOfTrackFront = TMath::Sqrt(position_X_Front*position_X_Front + position_Y_Front*position_Y_Front);
892 Double_t radialPositionOfTrackBack = TMath::Sqrt(position_X_Back*position_X_Back + position_Y_Back*position_Y_Back);
894 //---------------------------------------------------------------------------------------
896 Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut; // depends on the number of variables (here, 2)
898 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
900 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
901 AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId));
903 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
905 Bool_t isGoodChi2 = kFALSE;
907 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster);
908 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster); // describes the compatibility between the track and the cluster
909 if (chi2<chi2cut) isGoodChi2 = kTRUE;
911 Double_t radialPositionOfClusterFront = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
912 if (planeId == fNPlanesMFT-1) {
913 if (TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront)<fDistanceFromBestClusterAndTrackAtLastPlane ||
914 fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
915 fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
917 if (IsCorrectMatch(cluster)) {
918 fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
922 if (fIsCurrentMuonTrackable) {
923 // fOutputQAFile->cd();
924 if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf
925 else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf
929 AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
930 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
931 if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
932 newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
933 AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
934 planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
935 newTrack->SetPlaneExists(planeId);
936 AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
937 if (fIsCurrentMuonTrackable) {
938 Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
939 AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
940 Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope
941 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
942 // fOutputQAFile->cd();
943 if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
944 else fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
946 fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
948 else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
952 // Analyizing the clusters: BACK ACTIVE ELEMENTS
954 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
955 AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId));
957 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
959 Bool_t isGoodChi2 = kFALSE;
961 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster);
962 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster); // describes the compatibility between the track and the cluster
963 if (chi2<chi2cut) isGoodChi2 = kTRUE;
965 Double_t radialPositionOfClusterBack = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
966 if (planeId == fNPlanesMFT-1) {
967 if (TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack)<fDistanceFromBestClusterAndTrackAtLastPlane ||
968 fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
969 fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
971 if (IsCorrectMatch(cluster)) {
972 fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
976 if (fIsCurrentMuonTrackable) {
977 // fOutputQAFile->cd();
978 if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf
979 else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf
983 AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
984 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
985 if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return kDiverged;
986 newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
987 AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)",
988 planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters()));
989 newTrack->SetPlaneExists(planeId);
990 AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable));
991 if (fIsCurrentMuonTrackable) {
992 Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
993 AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2()));
994 Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope
995 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
996 // fOutputQAFile->cd();
997 if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
998 else fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
1000 fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY());
1002 else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut));
1006 //---------------------------------------------------------------------------------------------
1008 if (planeId == fNPlanesMFT-1) {
1009 if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) {
1010 // fOutputQAFile->cd();
1011 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane-
1012 fDistanceFromGoodClusterAndTrackAtLastPlane));
1013 fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane);
1021 //==========================================================================================================================================
1023 void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
1025 AliDebug(1, Form(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId));
1027 AliMUONTrackParam currentParamFront, currentParamBack;
1029 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
1030 currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
1031 currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
1032 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamFront, 0.);
1033 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamBack, 0.);
1035 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
1036 AliDebug(2, Form("fCurrentTrack = %p\n", fCurrentTrack));
1037 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
1038 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
1039 AliMUONTrackExtrap::AddMCSEffect(¤tParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
1040 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.);
1041 AliMUONTrackExtrap::AddMCSEffect(¤tParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
1042 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.);
1044 // for all planes: linear extrapolation to the Z of the plane
1045 AliMUONTrackExtrap::ExtrapToZCov(¤tParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
1046 AliMUONTrackExtrap::ExtrapToZCov(¤tParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
1048 Bool_t goodClusterFound = kFALSE;
1050 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
1052 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
1054 AliDebug(1, Form("nClustersFront = %d\n", nClustersFront));
1055 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
1056 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->UncheckedAt(iCluster);
1057 AliDebug(2, Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack));
1058 if (IsCorrectMatch(cluster)) {
1059 fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
1060 fCurrentTrack->SetPlaneExists(planeId);
1061 goodClusterFound = kTRUE;
1066 if (goodClusterFound) return;
1068 // Analyizing the clusters: BACK ACTIVE ELEMENTS
1070 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
1072 AliDebug(1, Form("nClustersBack = %d\n", nClustersBack));
1073 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
1074 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->UncheckedAt(iCluster);
1075 AliDebug(2,Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack));
1076 if (IsCorrectMatch(cluster)) {
1077 fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
1078 fCurrentTrack->SetPlaneExists(planeId);
1079 goodClusterFound = kTRUE;
1086 //==========================================================================================================================================
1088 void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() {
1090 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1091 fIsGoodClusterInPlane[iPlane] = kFALSE;
1092 Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
1093 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1094 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
1095 for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
1096 if (cluster->GetMCLabel(iTrack)==fLabelMC) {
1097 fIsGoodClusterInPlane[iPlane] = kTRUE;
1104 fIsCurrentMuonTrackable = kTRUE;
1105 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) fIsCurrentMuonTrackable = (fIsCurrentMuonTrackable&&fIsGoodClusterInPlane[iPlane]);
1109 //==========================================================================================================================================
1111 void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory() {
1113 // Fill planes with the clusters
1116 AliDebug(2, Form("fFinalBestCandidate->GetNMFTClusters() = %d\n", fFinalBestCandidate->GetNMFTClusters()));
1117 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1118 if (fFinalBestCandidate->PlaneExists(iPlane)) {
1119 AliMFTCluster *trackCluster = fFinalBestCandidate->GetMFTCluster(cluster++);
1120 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetPoint(fGrMFTPlane[kClusterOfTrack][iPlane]->GetN(), trackCluster->GetX(), trackCluster->GetY());
1122 Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
1123 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
1124 AliMFTCluster *myCluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster);
1125 fGrMFTPlane[kAllClusters][iPlane] -> SetPoint(fGrMFTPlane[kAllClusters][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
1126 if (IsCorrectMatch(myCluster)) {
1127 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetPoint(fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY());
1134 //======================================================================================================================================
1136 Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) {
1138 Bool_t result = kFALSE;
1140 // check if the cluster belongs to the correct MC track
1142 for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
1143 if (cluster->GetMCLabel(iTrack)==fLabelMC) {
1149 AliDebug(2,Form("returning %d\n", result));
1155 //======================================================================================================================================
1157 Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
1159 // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
1160 // return the corresponding Chi2
1161 // assume the track parameters are given at the Z of the cluster
1163 // Set differences between trackParam and cluster in the bending and non bending directions
1164 Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
1165 Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
1166 AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY));
1168 // Calculate errors and covariances
1169 const TMatrixD& kParamCov = trackParam.GetCovariances();
1170 Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
1171 Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
1172 AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2));
1173 Double_t covXY = kParamCov(0,2);
1174 Double_t det = sigmaX2 * sigmaY2 - covXY * covXY;
1177 if (det==0.) return 1.e10;
1178 return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
1182 //=========================================================================================================================================
1184 void AliMuonForwardTrackFinder::SeparateFrontBackClusters() {
1186 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1187 printf("Separating front/back clusters\n");
1188 fMFTClusterArrayFront[iPlane]->Delete();
1189 fMFTClusterArrayBack[iPlane] ->Delete();
1190 for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
1191 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
1192 if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
1193 new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
1196 new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
1203 //=========================================================================================================================================
1205 Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) {
1207 // the same definition as in AliMUONTrack is implemented, since here we just add more clusters to the Muon track
1209 Int_t ndf = 2 * nClusters - 5;
1210 return (ndf > 0) ? ndf : 0;
1214 //============================================================================================================================================
1216 void AliMuonForwardTrackFinder::BookHistos() {
1218 const Int_t nMaxNewTracks[] = {150, 200, 250, 600, 1000};
1219 const Double_t radiusPlane[] = {0.010, 0.010, 0.050, 0.5, 1.5};
1221 fHistRadiusEndOfAbsorber = new TH1D("hRadiusEndOfAbsorber", "Track radial distance at the end of the absorber", 1000, 0, 100.);
1223 fHistNGoodClustersForFinalTracks = new TH1D("hNGoodClustersForFinalTracks", "Number of Good Clusters per Final Track", 20, -0.25, 9.75);
1225 fHistDistanceGoodClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackAtLastPlane",
1226 "Distance of MC Good Cluster from Track in last MFT plane", 200, 0., 2.);
1228 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane =
1229 new TH1D("hDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane",
1230 "Good Cluster distance from track - Best Cluster distance from track in last MFT plane", 200, 0., 2.);
1232 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1234 fHistNTracksAfterExtrapolation[iPlane] = new TH1D(Form("hNTracksAfterExtrapolation_pl%02d", iPlane),
1235 Form("Number of Candidates after analysis of MFT plane %02d", iPlane),
1236 nMaxNewTracks[iPlane], -0.5, nMaxNewTracks[iPlane]-0.5);
1238 fHistResearchRadius[iPlane] = new TH1D(Form("hResearchRadius_pl%02d", iPlane),
1239 Form("Research Radius for candidate clusters in MFT plane %02d", iPlane),
1240 1000, 0., radiusPlane[iPlane]);
1242 fHistChi2Cluster_GoodCluster[iPlane] = new TH1D(Form("hChi2Cluster_GoodCluster_pl%02d", iPlane),
1243 Form("#chi^{2}_{clust} for Good clusters in MFT plane %02d", iPlane),
1246 fHistChi2Cluster_BadCluster[iPlane] = new TH1D(Form("hChi2Cluster_BadCluster_pl%02d", iPlane),
1247 Form("#chi^{2}_{clust} for Bad clusters in MFT plane %02d", iPlane),
1250 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1251 Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
1254 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1255 Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
1260 //------------------------------------------
1262 fHistRadiusEndOfAbsorber -> Sumw2();
1263 fHistNGoodClustersForFinalTracks -> Sumw2();
1265 fHistDistanceGoodClusterFromTrackAtLastPlane -> Sumw2();
1266 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Sumw2();
1268 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1270 fHistNTracksAfterExtrapolation[iPlane] -> Sumw2();
1271 fHistResearchRadius[iPlane] -> Sumw2();
1273 fHistChi2Cluster_GoodCluster[iPlane] -> Sumw2();
1274 fHistChi2Cluster_BadCluster[iPlane] -> Sumw2();
1276 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1277 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1281 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");
1283 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");
1287 //============================================================================================================================================
1289 void AliMuonForwardTrackFinder::SetTitleHistos() {
1291 fHistRadiusEndOfAbsorber -> SetXTitle("R_{abs} [cm]");
1292 fHistNGoodClustersForFinalTracks -> SetXTitle("N_{GoodClusters}");
1294 fHistDistanceGoodClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]");
1295 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]");
1298 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1300 fHistNTracksAfterExtrapolation[iPlane] -> SetXTitle("N_{tracks}");
1301 fHistResearchRadius[iPlane] -> SetXTitle("Research Radius [cm]");
1303 fHistChi2Cluster_GoodCluster[iPlane] -> SetXTitle("#chi^{2}/ndf");
1304 fHistChi2Cluster_BadCluster[iPlane] -> SetXTitle("#chi^{2}/ndf");
1306 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1307 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1313 //===========================================================================================================================================
1315 void AliMuonForwardTrackFinder::BookPlanes() {
1317 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1318 fGrMFTPlane[kAllClusters][iPlane] = new TGraph();
1319 fGrMFTPlane[kAllClusters][iPlane] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
1320 fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerStyle(20);
1321 // fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.5);
1322 // fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.3);
1323 fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.2);
1326 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1327 fGrMFTPlane[kClustersGoodChi2][iPlane] = new TGraph();
1328 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
1329 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerStyle(20);
1330 // fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.8);
1331 // fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.4);
1332 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.3);
1333 fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerColor(kBlue);
1336 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1337 fGrMFTPlane[kClusterOfTrack][iPlane] = new TGraph();
1338 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
1339 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerStyle(25);
1340 // fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(1.2);
1341 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(0.9);
1342 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerColor(kRed);
1343 fGrMFTPlane[kClusterOfTrack][iPlane] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
1346 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1347 fGrMFTPlane[kClusterCorrectMC][iPlane] = new TGraph();
1348 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
1349 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerStyle(20);
1350 // fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.8);
1351 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.5);
1352 fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerColor(kGreen);
1355 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1356 fCircleExt[iPlane] = new TEllipse(0., 0., fRPlaneMax[iPlane], fRPlaneMax[iPlane]);
1357 fCircleInt[iPlane] = new TEllipse(0., 0., fRPlaneMin[iPlane], fRPlaneMin[iPlane]);
1360 fTxtDummy = new TLatex(0.10, 0.59, "Best Candidate:");
1362 //---------------------------------------------------
1364 fMrkAllClust = new TMarker(0.10, 0.32, 20);
1365 fMrkAllClust -> SetMarkerSize(0.5);
1367 fMrkClustGoodChi2 = new TMarker(0.10, 0.26, 20);
1368 fMrkClustGoodChi2 -> SetMarkerSize(0.8);
1369 fMrkClustGoodChi2 -> SetMarkerColor(kBlue);
1371 fMrkClustMC = new TMarker(0.10, 0.20, 20);
1372 fMrkClustMC -> SetMarkerSize(0.8);
1373 fMrkClustMC -> SetMarkerColor(kGreen);
1375 fMrkClustOfTrack = new TMarker(0.10, 0.14, 25);
1376 fMrkClustOfTrack -> SetMarkerSize(1.2);
1377 fMrkClustOfTrack -> SetMarkerColor(kRed);
1379 fTxtAllClust = new TLatex(0.15, 0.30, "All Clusters");
1380 fTxtAllClust -> SetTextSize(0.040);
1382 fTxtClustGoodChi2 = new TLatex(0.15, 0.24, "Clusters involved in the research");
1383 fTxtClustGoodChi2 -> SetTextSize(0.040);
1385 fTxtClustMC = new TLatex(0.15, 0.18, "MC good clusters");
1386 fTxtClustMC -> SetTextSize(0.040);
1388 fTxtClustOfTrack = new TLatex(0.15, 0.12, "Clusters of the best candidate");
1389 fTxtClustOfTrack -> SetTextSize(0.040);
1393 //===========================================================================================================================================
1395 void AliMuonForwardTrackFinder::ResetPlanes() {
1397 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1398 for (Int_t iGr=0; iGr<4; iGr++) {
1399 Int_t nOldClusters = fGrMFTPlane[iGr][iPlane]->GetN();
1400 for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iGr][iPlane]->RemovePoint(iPoint);
1406 //===========================================================================================================================================
1408 void AliMuonForwardTrackFinder::PrintParticleHistory() {
1410 AliDebug(1, "Entering");
1412 TString history = "";
1414 TParticle *part = 0;
1415 if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
1417 AliDebug(1, Form("fStack->Particle(%d) = %p", fLabelMC, part));
1420 AliDebug(1, Form("fStack->Particle(%d)->GetPdgCode() = %d", fLabelMC, part->GetPdgCode()));
1421 if (part->GetFirstMother() != -1) {
1422 TParticle *partMother = fStack->Particle(part->GetFirstMother());
1423 AliDebug(1, Form("fStack->Particle(%d) = %p", part->GetFirstMother(), partMother));
1425 Char_t newName[100];
1426 if (partMother->GetFirstMother() != -1) history += "... #rightarrow ";
1427 PDGNameConverter(partMother->GetName(), newName);
1428 history += Form("%s #rightarrow ", newName);
1431 Char_t newName[100];
1432 AliDebug(1, Form("fStack->Particle(%d)->GetPdgCode() = %d", fLabelMC, part->GetPdgCode()));
1433 AliDebug(1, Form("fStack->Particle(%d)->GetName() = %s", fLabelMC, part->GetName()));
1434 PDGNameConverter(part->GetName(), newName);
1435 history += Form("%s at z = %5.1f cm", newName, part->Vz());
1436 // printf("%s", history.Data());
1438 else history += "NO AVAILABLE HISTORY";
1440 fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data());
1442 // Filling particle history in the fFinalBestCandidate
1445 for (Int_t iParent=0; iParent<AliMuonForwardTrack::fgkNParentsMax; iParent++) {
1446 if (part->GetFirstMother() == -1) break;
1447 if (!(fStack->Particle(part->GetFirstMother()))) break;
1448 AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", fStack->Particle(part->GetFirstMother())));
1449 fFinalBestCandidate->SetParentMCLabel(iParent, part->GetFirstMother());
1450 fFinalBestCandidate->SetParentPDGCode(iParent, fStack->Particle(part->GetFirstMother())->GetPdgCode());
1451 part = fStack->Particle(part->GetFirstMother());
1457 //===========================================================================================================================================
1459 Bool_t AliMuonForwardTrackFinder::IsMother(const Char_t *nameMother) {
1461 Bool_t result = kFALSE;
1463 TParticle *part = 0;
1464 if (fLabelMC>=0) part = fStack->Particle(fLabelMC);
1467 if (part->GetFirstMother() != -1) {
1468 TParticle *partMother = fStack->Particle(part->GetFirstMother());
1470 if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
1479 //===========================================================================================================================================
1481 void AliMuonForwardTrackFinder::DrawPlanes() {
1484 if (fNPlanesMFT <= 5) fCanvas -> Divide(3,2);
1485 else if (fNPlanesMFT <= 11) fCanvas -> Divide(4,3);
1486 else if (fNPlanesMFT <= 19) fCanvas -> Divide(5,4);
1488 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1490 fCanvas->cd(fNPlanesMFT-iPlane+1);
1492 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1493 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1494 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetTitle("X [cm]");
1495 fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetTitle("Y [cm]");
1496 fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("ap");
1498 fCircleExt[iPlane] -> Draw("same");
1499 fCircleInt[iPlane] -> Draw("same");
1501 if (fGrMFTPlane[kAllClusters][iPlane]->GetN()) fGrMFTPlane[kAllClusters][iPlane] -> Draw("psame");
1502 if (fGrMFTPlane[kClustersGoodChi2][iPlane]->GetN()) fGrMFTPlane[kClustersGoodChi2][iPlane] -> Draw("psame");
1503 if (fGrMFTPlane[kClusterOfTrack][iPlane]->GetN()) fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("psame");
1504 if (fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN()) fGrMFTPlane[kClusterCorrectMC][iPlane] -> Draw("psame");
1506 fTxtTrackChi2[iPlane] -> Draw("same");
1511 fTxtMuonHistory -> Draw();
1512 fTxtDummy -> Draw("same");
1513 if (fMatchingMode==kRealMatching) fTxtTrackGoodClusters -> Draw("same");
1514 fTxtTrackFinalChi2 -> Draw("same");
1515 fTxtTrackMomentum -> Draw("same");
1516 if (fMatchingMode==kRealMatching) fTxtFinalCandidates -> Draw("same");
1518 fMrkAllClust -> Draw("same");
1519 fMrkClustGoodChi2 -> Draw("same");
1520 fMrkClustMC -> Draw("same");
1521 fMrkClustOfTrack -> Draw("same");
1523 fTxtAllClust -> Draw("same");
1524 fTxtClustGoodChi2 -> Draw("same");
1525 fTxtClustMC -> Draw("same");
1526 fTxtClustOfTrack -> Draw("same");
1528 // fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1529 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1530 if (IsMother("phi")) {
1531 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1532 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1534 if (IsMother("J/psi")) {
1535 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1536 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1541 //===========================================================================================================================================
1543 void AliMuonForwardTrackFinder::Terminate() {
1546 AliInfo("---------------------------------------------------------------------------------------------------------------");
1547 AliInfo(Form("%8d tracks analyzed", fCountRealTracksAnalyzed));
1548 AliInfo(Form("%8d tracks with MC ref", fCountRealTracksWithRefMC));
1549 AliInfo(Form("%8d tracks with MC ref & trigger match", fCountRealTracksWithRefMC_andTrigger));
1550 if (fMatchingMode==kRealMatching) {
1551 AliInfo(Form("%8d tracks analyzed with final candidates", fCountRealTracksAnalyzedWithFinalCandidates));
1554 AliInfo(Form("%8d tracks matched with their MC clusters", fCountRealTracksAnalyzedWithFinalCandidates));
1556 // printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c", fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
1557 // printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
1558 AliInfo("---------------------------------------------------------------------------------------------------------------");
1565 //==========================================================================================================================================
1567 void AliMuonForwardTrackFinder::FillOutputTree() {
1569 if (!fMuonForwardTracks || !fOutputEventTree) return;
1571 AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p",
1572 fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0)));
1574 // fOutputTreeFile->cd();
1575 fOutputEventTree->Fill();
1576 AliDebug(1, Form("\nFilled Tree: nEvents = %d!!!!\n", Int_t(fOutputEventTree->GetEntries())));
1580 //==========================================================================================================================================
1582 void AliMuonForwardTrackFinder::WriteOutputTree() {
1584 if (!fOutputEventTree || !fOutputTreeFile) return;
1586 fOutputTreeFile -> cd();
1588 fOutputEventTree -> Write();
1589 fOutputTreeFile -> Close();
1593 //==========================================================================================================================================
1595 void AliMuonForwardTrackFinder::WriteHistos() {
1597 fOutputQAFile = new TFile(Form("MuonGlobalTracking.QA.run%d.root", fRun), "recreate");
1598 fOutputQAFile -> cd();
1600 fHistRadiusEndOfAbsorber -> Write();
1601 fHistNGoodClustersForFinalTracks -> Write();
1603 fHistDistanceGoodClusterFromTrackAtLastPlane -> Write();
1604 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Write();
1606 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1608 fHistNTracksAfterExtrapolation[iPlane] -> Write();
1609 fHistResearchRadius[iPlane] -> Write();
1611 fHistChi2Cluster_GoodCluster[iPlane] -> Write();
1612 fHistChi2Cluster_BadCluster[iPlane] -> Write();
1614 fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
1615 fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Write();
1619 fNtuFinalCandidates -> Write();
1620 fNtuFinalBestCandidates -> Write();
1622 fOutputQAFile -> Close();
1626 //===========================================================================================================================================
1628 void AliMuonForwardTrackFinder::PDGNameConverter(const Char_t *nameIn, Char_t *nameOut) {
1630 if (!strcmp(nameIn, "mu+")) snprintf(nameOut, 50, "#mu^{+}");
1631 else if (!strcmp(nameIn, "mu-")) snprintf(nameOut, 50, "#mu^{-}");
1632 else if (!strcmp(nameIn, "pi+")) snprintf(nameOut, 50, "#pi^{+}");
1633 else if (!strcmp(nameIn, "pi-")) snprintf(nameOut, 50, "#pi^{-}");
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*+")) snprintf(nameOut, 50, "K^{*+}");
1637 else if (!strcmp(nameIn, "K*-")) snprintf(nameOut, 50, "K^{*-}");
1638 else if (!strcmp(nameIn, "K_S0")) snprintf(nameOut, 50, "K_{S}^{0}");
1639 else if (!strcmp(nameIn, "K_L0")) snprintf(nameOut, 50, "K_{L}^{0}");
1640 else if (!strcmp(nameIn, "K0")) snprintf(nameOut, 50, "K^{0}");
1641 else if (!strcmp(nameIn, "K0_bar")) snprintf(nameOut, 50, "#bar{K}^{0}");
1642 else if (!strcmp(nameIn, "K*0")) snprintf(nameOut, 50, "K^{*0}");
1643 else if (!strcmp(nameIn, "K*0_bar")) snprintf(nameOut, 50, "#bar{K}^{*0}");
1644 else if (!strcmp(nameIn, "rho0")) snprintf(nameOut, 50, "#rho^{0}");
1645 else if (!strcmp(nameIn, "rho+")) snprintf(nameOut, 50, "#rho^{+}");
1646 else if (!strcmp(nameIn, "rho-")) snprintf(nameOut, 50, "#rho^{-}");
1647 else if (!strcmp(nameIn, "omega")) snprintf(nameOut, 50, "#omega");
1648 else if (!strcmp(nameIn, "eta'")) snprintf(nameOut, 50, "#eta'");
1649 else if (!strcmp(nameIn, "phi")) snprintf(nameOut, 50, "#phi");
1651 else if (!strcmp(nameIn, "D-")) snprintf(nameOut, 50, "D^{-}");
1652 else if (!strcmp(nameIn, "D+")) snprintf(nameOut, 50, "D^{+}");
1653 else if (!strcmp(nameIn, "D0")) snprintf(nameOut, 50, "D^{0}");
1654 else if (!strcmp(nameIn, "D0_bar")) snprintf(nameOut, 50, "#bar{D}^{0}");
1655 else if (!strcmp(nameIn, "D*-")) snprintf(nameOut, 50, "D^{*-}");
1656 else if (!strcmp(nameIn, "D*+")) snprintf(nameOut, 50, "D^{*+}");
1657 else if (!strcmp(nameIn, "D_s+")) snprintf(nameOut, 50, "D_{s}^{+}");
1658 else if (!strcmp(nameIn, "D*_s+")) snprintf(nameOut, 50, "D_{s}^{*+}");
1660 else if (!strcmp(nameIn, "B-")) snprintf(nameOut, 50, "B^{-}");
1661 else if (!strcmp(nameIn, "B+")) snprintf(nameOut, 50, "B^{+}");
1662 else if (!strcmp(nameIn, "B_s0_bar")) snprintf(nameOut, 50, "#bar{B}_{s}^{0}");
1664 else if (!strcmp(nameIn, "antiproton")) snprintf(nameOut, 50, "#bar{p}");
1665 else if (!strcmp(nameIn, "proton")) snprintf(nameOut, 50, "p");
1666 else if (!strcmp(nameIn, "neutron")) snprintf(nameOut, 50, "n");
1667 else if (!strcmp(nameIn, "Sigma+")) snprintf(nameOut, 50, "#Sigma^{+}");
1668 else if (!strcmp(nameIn, "Delta+")) snprintf(nameOut, 50, "#Delta{+}");
1669 else if (!strcmp(nameIn, "Delta--")) snprintf(nameOut, 50, "#Delta{--}");
1670 else if (!strcmp(nameIn, "Lambda0")) snprintf(nameOut, 50, "#Lambda_0");
1671 else if (!strcmp(nameIn, "Lambda0_bar")) snprintf(nameOut, 50, "#bar{Lambda}_0");
1673 else snprintf(nameOut, 50, "%s", nameIn);
1677 //===========================================================================================================================================
1679 void AliMuonForwardTrackFinder::SetDraw(Bool_t drawOption) {
1681 fDrawOption = drawOption;
1684 fCanvas = new TCanvas("tracking", "tracking", 1200, 800);
1685 fCanvas -> Divide(3,2);
1690 //===========================================================================================================================================
1692 Bool_t AliMuonForwardTrackFinder::InitGRP() {
1694 //------------------------------------
1695 // Initialization of the GRP entry
1696 //------------------------------------
1698 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
1702 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
1705 AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
1707 fGRPData = new AliGRPObject();
1708 fGRPData->ReadValuesFromMap(m);
1712 AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
1713 fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
1717 // FIX ME: The unloading of GRP entry is temporarily disabled
1718 // because ZDC and VZERO are using it in order to initialize
1719 // their reconstructor objects. In the future one has to think
1720 // of propagating AliRunInfo to the reconstructors.
1721 // AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
1725 AliError("No GRP entry found in OCDB!");
1729 TString lhcState = fGRPData->GetLHCState();
1730 if (lhcState==AliGRPObject::GetInvalidString()) {
1731 AliError("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN");
1732 lhcState = "UNKNOWN";
1735 TString beamType = fGRPData->GetBeamType();
1736 if (beamType==AliGRPObject::GetInvalidString()) {
1737 AliError("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN");
1738 beamType = "UNKNOWN";
1741 Float_t beamEnergy = fGRPData->GetBeamEnergy();
1742 if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
1743 AliError("GRP/GRP/Data entry: missing value for the beam energy ! Using 0");
1747 TString runType = fGRPData->GetRunType();
1748 if (runType==AliGRPObject::GetInvalidString()) {
1749 AliError("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN");
1750 runType = "UNKNOWN";
1753 Int_t activeDetectors = fGRPData->GetDetectorMask();
1754 if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
1755 AliError("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399");
1756 activeDetectors = 1074790399;
1758 AliDebug(1, Form("activeDetectors = %d", activeDetectors));
1760 fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
1763 // *** Dealing with the magnetic field map
1765 if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
1766 if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
1767 AliInfo("ExpertMode!!! GRP information will be ignored !");
1768 AliInfo("ExpertMode!!! Running with the externally locked B field !");
1771 AliInfo("Destroying existing B field instance!");
1772 delete TGeoGlobalMagField::Instance();
1775 if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
1776 // Construct the field map out of the information retrieved from GRP.
1779 Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
1780 if (l3Current == AliGRPObject::GetInvalidFloat()) {
1781 AliError("GRP/GRP/Data entry: missing value for the L3 current !");
1785 Char_t l3Polarity = fGRPData->GetL3Polarity();
1786 if (l3Polarity == AliGRPObject::GetInvalidChar()) {
1787 AliError("GRP/GRP/Data entry: missing value for the L3 polarity !");
1792 Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
1793 if (diCurrent == AliGRPObject::GetInvalidFloat()) {
1794 AliError("GRP/GRP/Data entry: missing value for the dipole current !");
1798 Char_t diPolarity = fGRPData->GetDipolePolarity();
1799 if (diPolarity == AliGRPObject::GetInvalidChar()) {
1800 AliError("GRP/GRP/Data entry: missing value for the dipole polarity !");
1804 // read special bits for the polarity convention and map type
1805 Int_t polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
1806 Bool_t uniformB = fGRPData->IsUniformBMap();
1809 AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
1810 TMath::Abs(diCurrent) * (diPolarity ? -1:1),
1811 polConvention,uniformB,beamEnergy, beamType.Data());
1813 TGeoGlobalMagField::Instance()->SetField( fld );
1814 TGeoGlobalMagField::Instance()->Lock();
1815 AliInfo("Running with the B field constructed out of GRP !");
1817 else AliFatal("Failed to create a B field map !");
1819 else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1825 //====================================================================================================================================================
1827 Bool_t AliMuonForwardTrackFinder::SetRunNumber() {
1829 AliCDBManager *man = AliCDBManager::Instance();
1832 AliError("No run loader found!");
1836 fRunLoader->LoadHeader();
1837 // read run number from gAlice
1838 if (fRunLoader->GetHeader()) {
1839 man->SetRun(fRunLoader->GetHeader()->GetRun());
1840 fRunLoader->UnloadHeader();
1843 AliError("No run-loader header found!");
1852 //====================================================================================================================================================