// ROOT includes #include "TObject.h" #include "TClonesArray.h" #include "TObjArray.h" #include "TH2D.h" #include "TH1D.h" #include "TFile.h" #include "TGeoManager.h" #include "TMatrixD.h" #include "TParticle.h" #include "TMath.h" #include "TGraph.h" #include "TEllipse.h" #include "TCanvas.h" #include "TString.h" #include "TLatex.h" #include "TMarker.h" #include "TNtuple.h" #include "TRandom.h" #include "TIterator.h" // STEER includes #include "AliLog.h" #include "AliRun.h" #include "AliRunLoader.h" #include "AliLoader.h" #include "AliHeader.h" #include "AliMC.h" #include "AliStack.h" #include "AliMagF.h" #include "AliTracker.h" #include "AliGRPObject.h" #include "AliCDBEntry.h" #include "AliCDBManager.h" // MUON includes #include "AliMUONConstants.h" #include "AliMUONTrack.h" #include "AliMUONRecoCheck.h" #include "AliMUONTrackParam.h" #include "AliMUONTrackExtrap.h" #include "AliMUONVTrackStore.h" #include "AliMUONVCluster.h" // MFT includes #include "AliMuonForwardTrack.h" #include "AliMFTCluster.h" #include "AliMFT.h" #include "AliMFTSegmentation.h" #include "AliMFTConstants.h" #include "AliMuonForwardTrackFinder.h" //==================================================================================================================================================== // // Class for the creation of the "global muon tracks" built from the clusters in the // muon spectrometer and the clusters of the Muon Forward Tracker. QA histograms are also created // // Contact author: antonio.uras@cern.ch // //==================================================================================================================================================== const Double_t AliMuonForwardTrackFinder::fRadLengthSi = AliMFTConstants::fRadLengthSi; ClassImp(AliMuonForwardTrackFinder) //===================================================================================================== AliMuonForwardTrackFinder::AliMuonForwardTrackFinder(): TObject(), fRun(0), fNEventsToAnalyze(0), fSigmaClusterCut(0), fScaleSigmaClusterCut(1.), fGlobalTrackingDiverged(kFALSE), fChi2GlobalCut(0), fSigmaSpectrometerCut(0), fVertexErrorX(0.015), fVertexErrorY(0.015), fVertexErrorZ(0.010), fNFinalCandidatesCut(0), fReadDir(0), fOutDir(0), fDrawOption(0), fDistanceFromGoodClusterAndTrackAtLastPlane(-1), fDistanceFromBestClusterAndTrackAtLastPlane(-1), fRAbsorberCut(0), fLowPtCut(0), fNPlanesMFT(0), fNPlanesMFTAnalyzed(0), fNMaxMissingMFTClusters(0), fEv(0), fLabelMC(0), fHistRadiusEndOfAbsorber(0), fHistNGoodClustersForFinalTracks(0), fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane(0), fHistDistanceGoodClusterFromTrackAtLastPlane(0), fNtuFinalCandidates(0), fNtuFinalBestCandidates(0), fCanvas(0), fTxtMuonHistory(0), fTxtTrackGoodClusters(0), fTxtTrackFinalChi2(0), fTxtTrackMomentum(0), fTxtFinalCandidates(0), fTxtDummy(0), fTxtAllClust(0), fTxtClustGoodChi2(0), fTxtClustMC(0), fTxtClustOfTrack(0), fMrkAllClust(0), fMrkClustGoodChi2(0), fMrkClustMC(0), fMrkClustOfTrack(0), fCountRealTracksAnalyzed(0), fMaxNTracksToBeAnalyzed(99999999), fCountRealTracksWithRefMC(0), fCountRealTracksWithRefMC_andTrigger(0), fCountRealTracksWithRefMC_andTrigger_andGoodPt(0), fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta(0), fCountRealTracksAnalyzedOfEvent(0), fCountRealTracksAnalyzedWithFinalCandidates(0), fFileCluster(0), fFileESD(0), fFile_gAlice(0), fRunLoader(0), fMFTLoader(0), fMuonRecoCheck(0), fMFTClusterTree(0), fMuonTrackReco(0), fCurrentTrack(0), fFinalBestCandidate(0), fIsCurrentMuonTrackable(0), fCandidateTracks(0), fTrackStore(0), fTrackRefStore(0), fNextTrack(0), fStack(0), fMFT(0), fSegmentation(0), fOutputTreeFile(0), fOutputQAFile(0), fOutputEventTree(0), fMuonForwardTracks(0), fMatchingMode(-1), fGRPData(0), fRunInfo(0), fBransonCorrection(kTRUE) { // Default constructor for (Int_t iPlane=0; iPlane Branch("tracks", &fMuonForwardTracks); } //===================================================================================================== AliMuonForwardTrackFinder::~AliMuonForwardTrackFinder() { for (Int_t iPlane=0; iPlaneIsOpen()) return; else AliInfo(Form("file %s successfully opened\n", fFileESD->GetName())); fMuonRecoCheck = new AliMUONRecoCheck(esdFileName, Form("%s/generated/", fReadDir.Data())); // Utility class to check reconstruction fFile_gAlice = new TFile(gAliceName); if (!fFile_gAlice || !fFile_gAlice->IsOpen()) return; else AliInfo(Form("file %s successfully opened\n", fFile_gAlice->GetName())); fRunLoader = AliRunLoader::Open(gAliceName); gAlice = fRunLoader->GetAliRun(); if (!gAlice) fRunLoader->LoadgAlice(); fMFT = (AliMFT*) gAlice->GetDetector("MFT"); fSegmentation = fMFT->GetSegmentation(); SetNPlanesMFT(fSegmentation->GetNPlanes()); if (!SetRunNumber()) return; if (!InitGRP()) return; AliMUONTrackExtrap::SetField(); // set the magnetic field for track extrapolations for (Int_t iPlane=0; iPlaneGetPlane(iPlane)->GetZCenter(); fRPlaneMax[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMaxSupport(); fRPlaneMin[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMinSupport(); } // Loading MFT clusters fMFTLoader = fRunLoader->GetDetectorLoader("MFT"); fMFTLoader->LoadRecPoints("READ"); fMFTClusterTree = fMFTLoader->TreeR(); Int_t nEventsInFile = fMuonRecoCheck->NumberOfEvents(); if (!nEventsInFile) { AliError("no events available!!!\n"); return; } if (nEventsInFile=fNEventsToAnalyze) return kFALSE; fCountRealTracksAnalyzedOfEvent = 0; AliInfo(Form(" **** analyzing event # %d \n", fEv)); fTrackStore = fMuonRecoCheck->ReconstructedTracks(fEv); if (fTrackStore->IsEmpty()) { AliInfo("fTrackStore Is Empty: exiting NOW!"); return kFALSE; } AliInfo("fTrackStore contains tracks!"); AliDebug(2, Form("Getting fMuonRecoCheck->ReconstructibleTracks(%d)", fEv)); fTrackRefStore = fMuonRecoCheck->ReconstructibleTracks(fEv); AliDebug(2, Form("Getting fRunLoader->GetEvent(%d)", fEv)); fRunLoader->GetEvent(fEv); AliDebug(2, Form("fMFTLoader->TreeR() = %p",fMFTLoader->TreeR())); if (!fMFTLoader->TreeR()->GetEvent()) return kFALSE; for (Int_t iPlane=0; iPlaneGetRecPointsList(iPlane))->GetEntries())); fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane); } SeparateFrontBackClusters(); fRunLoader -> LoadKinematics(); fStack = fRunLoader->Stack(); fNextTrack = fTrackStore->CreateIterator(); fMuonForwardTracks->Delete(); fEv++; return kTRUE; } //====================================================================================================================================== Int_t AliMuonForwardTrackFinder::LoadNextTrack() { fNPlanesMFTAnalyzed = 0; // load next muon track from the reconstructed event if (fCountRealTracksAnalyzed>=fMaxNTracksToBeAnalyzed) return kFALSE; if (!fCountRealTracksAnalyzed) if (!LoadNextEvent()) return kFALSE; if (!fGlobalTrackingDiverged) { while ( !(fMuonTrackReco = static_cast(fNextTrack->Next())) ) if (!LoadNextEvent()) return kFALSE; fCountRealTracksAnalyzed++; fCountRealTracksAnalyzedOfEvent++; } AliDebug(1, "**************************************************************************************\n"); AliDebug(1, Form("*************************** MUON TRACK %3d ***************************************\n", fCountRealTracksAnalyzedOfEvent)); AliDebug(1, "**************************************************************************************\n"); fCandidateTracks -> Delete(); fLabelMC = -1; fDistanceFromGoodClusterAndTrackAtLastPlane = -1.; fDistanceFromBestClusterAndTrackAtLastPlane = -1.; ResetPlanes(); TIter nextTrackRef(fTrackRefStore->CreateIterator()); AliMUONTrack *trackRef=0; // --------------------------------------- loop on MC generated tracks to find the MC reference... while ( (trackRef = static_cast(nextTrackRef())) ) { // number of compatible clusters between trackReco and trackRef Int_t nMatchCluster = fMuonTrackReco->FindCompatibleClusters(*trackRef, fSigmaSpectrometerCut, fIsClusterCompatible); if ( (fIsClusterCompatible[0] || fIsClusterCompatible[1] || fIsClusterCompatible[2] || fIsClusterCompatible[3]) && // before the dipole (fIsClusterCompatible[6] || fIsClusterCompatible[7] || fIsClusterCompatible[8] || fIsClusterCompatible[9]) && // after the dipole 2*nMatchCluster>fMuonTrackReco->GetNClusters() ) { fMuonTrackReco->SetMCLabel(trackRef->GetUniqueID()); // MC reference has been found for trackReco! break; } } // ------------------------------------- ...done! fLabelMC = fMuonTrackReco->GetMCLabel(); Int_t motherPdg=0; if (fLabelMC>=0) { if (!fGlobalTrackingDiverged) fCountRealTracksWithRefMC++; if (fStack->Particle(fLabelMC)->GetFirstMother() != -1) { motherPdg = fStack->Particle(fStack->Particle(fLabelMC)->GetFirstMother())->GetPdgCode(); } } CheckCurrentMuonTrackable(); if (!fGlobalTrackingDiverged) if (fMuonTrackReco->GetMatchTrigger()) fCountRealTracksWithRefMC_andTrigger++; // the track we are going to build, starting from fMuonTrackReco and adding the MFT clusters AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack(); track -> SetMUONTrack(new AliMUONTrack(*fMuonTrackReco)); if (fLabelMC>=0 && fStack->Particle(fLabelMC)) track->SetMCTrackRef(new TParticle(*(fStack->Particle(fLabelMC)))); track -> SetMCLabel(fMuonTrackReco->GetMCLabel()); track -> SetMatchTrigger(fMuonTrackReco->GetMatchTrigger()); // track origin Double_t xVtx=-999., yVtx=-999., zVtx=-999999.; if (track->GetMCTrackRef()) { xVtx = track->GetMCTrackRef()->Vx(); yVtx = track->GetMCTrackRef()->Vy(); zVtx = track->GetMCTrackRef()->Vz(); } // track kinematics Double_t pt=-999., theta=-999., eta=-999.; if (track->GetMCTrackRef()) { pt = track->GetMCTrackRef()->Pt(); theta = track->GetMCTrackRef()->Theta(); if (theta<0.) theta += TMath::Pi(); eta = track->GetMCTrackRef()->Eta(); } else { AliMUONTrackParam *param = (AliMUONTrackParam*) (fMuonTrackReco->GetTrackParamAtCluster()->First()); pt = TMath::Sqrt(param->Px()*param->Px() + param->Py()*param->Py()); theta = TMath::ATan(pt/param->Pz()); if (theta<0.) theta += TMath::Pi(); eta = -1.*TMath::Log(TMath::Tan(0.5*theta)); } // if the transverse momentum is smaller than the threshold, skip to the next track if (pt < fLowPtCut) return 3; // track parameters linearly extrapolated from the first tracking station to the end of the absorber AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First()))); AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, -503.); // absorber extends from -90 to -503 cm Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor(); Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor(); Double_t rAbsorber = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber); fHistRadiusEndOfAbsorber -> Fill(rAbsorber); // if the radial distance of the track at the end of the absorber is smaller than a given radius, skip to the next track if (rAbsorber < fRAbsorberCut) return 4; //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS --------------------------------------- for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) { // *** do not reverse the order of this cycle!!! // *** this reflects the fact that the extrapolation is performed // *** starting from the last MFT plane back to the origin // --------- updating the array of tracks according to the clusters available in the i-th plane --------- fNPlanesMFTAnalyzed++; if (fMatchingMode==kRealMatching) { Int_t nTracksToBeAnalyzed = fCandidateTracks->GetEntriesFast(); for (Int_t iTrack=0; iTrackUncheckedAt(iTrack); // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array // (several new tracks can be created for one old track) if (FindClusterInPlane(iPlane) == 2) { fGlobalTrackingDiverged = kTRUE; if (fScaleSigmaClusterCut>0) fScaleSigmaClusterCut -= 0.1; return 6; } if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) { fCandidateTracks->Remove(fCurrentTrack); // the old track is removed after the check; } } fCandidateTracks->Compress(); if (fIsCurrentMuonTrackable) { // fOutputQAFile->cd(); fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast()); } } else if (fMatchingMode==kIdealMatching && fIsCurrentMuonTrackable) { fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0); AliDebug(2, Form("plane %02d: fCandidateTracks->GetEntriesFast() = %d fCandidateTracks->UncheckedAt(0) = %p fCurrentTrack = %p\n", iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack)); AttachGoodClusterInPlane(iPlane); } } // -------------------------- END OF THE CYCLE OVER THE MFT PLANES -------------------------------------------- fGlobalTrackingDiverged = kFALSE; fScaleSigmaClusterCut = 1.0; AliDebug(1, "Finished cycle over planes"); Double_t momentum = pt * TMath::CosH(eta); fTxtTrackMomentum = new TLatex(0.10, 0.70, Form("P_{spectro} = %3.1f GeV/c", momentum)); if (fMatchingMode==kIdealMatching) { AliDebug(1, "Adding track to output tree...\n"); fFinalBestCandidate = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0); AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0); new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack); AliDebug(1, "...track added!\n"); fCandidateTracks->Delete(); fCountRealTracksAnalyzedOfEvent++; fCountRealTracksAnalyzedWithFinalCandidates++; PrintParticleHistory(); FillPlanesWithTrackHistory(); Double_t chi2AtPlane[fNMaxPlanes] = {0}; Int_t nGoodClusters = 0; Int_t nMFTClusters = fFinalBestCandidate->GetNMFTClusters(); // Int_t nMUONClusters = fFinalBestCandidate->GetNMUONClusters(); Int_t plane = 0; for (Int_t iCluster=0; iClusterPlaneExists(plane)) plane++; AliMFTCluster *localCluster = fFinalBestCandidate->GetMFTCluster(iCluster); chi2AtPlane[plane] = localCluster->GetLocalChi2(); if (IsCorrectMatch(localCluster)) nGoodClusters++; // Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope // Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack); // chi2AtPlane[plane] /= Double_t(ndfGlobalTrack); plane++; } for (Int_t iPlane=0; iPlaneGetEntriesFast(); AliDebug(1, Form("nFinalTracks = %d", nFinalTracks)); if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++; Double_t theVariable_Best = -1.; // variable defining the best candidate Bool_t bestCandidateExists = kFALSE; Int_t nGoodClustersBestCandidate = 0; Int_t idBestCandidate = 0; Double_t chi2HistoryForBestCandidate[fNMaxPlanes] = {0}; // chi2 on each plane, for the best candidate Double_t nClustersPerPlane[fNMaxPlanes] = {0}; for (Int_t iPlane=0; iPlane GetEntries(); } fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks)); Int_t nClustersMC = 0; for (Int_t iPlane=0; iPlaneUncheckedAt(iTrack); Double_t chi2AtPlane[fNMaxPlanes] = {0}; Int_t nGoodClusters = 0; Int_t nMFTClusters = finalTrack->GetNMFTClusters(); // Int_t nMUONClusters = finalTrack->GetNMUONClusters(); Int_t plane = 0; for (Int_t iCluster=0; iClusterPlaneExists(plane)) plane++; AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster); chi2AtPlane[plane] = localCluster->GetLocalChi2(); if (IsCorrectMatch(localCluster)) nGoodClusters++; // Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope // Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack); // chi2AtPlane[plane] /= Double_t(ndfGlobalTrack); plane++; } if (fIsCurrentMuonTrackable) { // fOutputQAFile->cd(); fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters); } // fOutputQAFile->cd(); Float_t finalCandidatesInfo[] = {Double_t(fRun), Double_t(fEv), Double_t(fCountRealTracksAnalyzedOfEvent), Double_t(nFinalTracks), Double_t(fLabelMC>=0), xVtx, yVtx, zVtx, motherPdg, Double_t(fMuonTrackReco->GetMatchTrigger()), Double_t(nClustersMC), Double_t(nGoodClusters), pt, theta, eta, chi2AtPlane[0], chi2AtPlane[1], chi2AtPlane[2], chi2AtPlane[3], chi2AtPlane[4], chi2AtPlane[5], chi2AtPlane[6], chi2AtPlane[7], chi2AtPlane[8]}; fNtuFinalCandidates -> Fill(finalCandidatesInfo); // now comparing the tracks with various criteria, in order to find the best one Double_t theVariable = 0.; // theVariable = chi2AtPlane[0]; for (Int_t iCluster=0; iClusterUncheckedAt(idBestCandidate); AliInfo(Form("fFinalBestCandidate->GetNMFTClusters() = %d\n", fFinalBestCandidate->GetNMFTClusters())); PrintParticleHistory(); FillPlanesWithTrackHistory(); AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate); newTrack -> SetNWrongClustersMC(newTrack->GetNMFTClusters() - nGoodClustersBestCandidate); newTrack -> SetTrackMCId(fRun*100000+fEv*1000+fCountRealTracksAnalyzedOfEvent); new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack); } // fOutputQAFile->cd(); Float_t finalBestCandidatesInfo[] = {Double_t(fRun), Double_t(fEv), Double_t(fCountRealTracksAnalyzedOfEvent), Double_t(nFinalTracks), Double_t(fLabelMC>=0), xVtx, yVtx, zVtx, motherPdg, Double_t(fMuonTrackReco->GetMatchTrigger()), Double_t(nClustersMC), Double_t(nGoodClustersBestCandidate), pt, theta, eta, chi2HistoryForBestCandidate[0], chi2HistoryForBestCandidate[1], chi2HistoryForBestCandidate[2], chi2HistoryForBestCandidate[3], chi2HistoryForBestCandidate[4], chi2HistoryForBestCandidate[5], chi2HistoryForBestCandidate[6], chi2HistoryForBestCandidate[7], chi2HistoryForBestCandidate[8], nClustersPerPlane[0], nClustersPerPlane[1], nClustersPerPlane[2], nClustersPerPlane[3], nClustersPerPlane[4], nClustersPerPlane[5], nClustersPerPlane[6], nClustersPerPlane[7], nClustersPerPlane[8]}; fNtuFinalBestCandidates -> Fill(finalBestCandidatesInfo); if (fDrawOption && bestCandidateExists) { fTxtTrackGoodClusters = new TLatex(0.20, 0.51, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate)); DrawPlanes(); } // ------------------------------------------------------------------------------------------- fCandidateTracks->Delete(); fFinalBestCandidate = NULL; return 5; } //=========================================================================================================================================== Int_t AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) { AliDebug(2, Form(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId)); // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available // propagate track to plane #planeId (both to front and back active sensors) // look for compatible clusters // update TrackParam at found cluster (if any) using Kalman Filter AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack; if (planeId == fNPlanesMFT-1) { // last plane of the telecope currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First()))); currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First()))); currentParamForResearchFront = currentParamFront; currentParamForResearchBack = currentParamBack; Double_t xExtrap = gRandom->Gaus(0,fVertexErrorX); Double_t yExtrap = gRandom->Gaus(0,fVertexErrorY); Double_t zExtrap = gRandom->Gaus(0,fVertexErrorZ); if (fBransonCorrection) { AliMUONTrackExtrap::ExtrapToVertex(¤tParamFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); AliMUONTrackExtrap::ExtrapToVertex(¤tParamBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); } else { AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamFront, zExtrap); AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamBack, zExtrap); } AliMUONTrackExtrap::ExtrapToVertex(¤tParamForResearchFront, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); AliMUONTrackExtrap::ExtrapToVertex(¤tParamForResearchBack, xExtrap, yExtrap, zExtrap, fVertexErrorX, fVertexErrorY); } else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First()))); currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First()))); currentParamForResearchFront = currentParamFront; currentParamForResearchBack = currentParamBack; AliMUONTrackExtrap::AddMCSEffect(¤tParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+ fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.); AliMUONTrackExtrap::AddMCSEffect(¤tParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+ fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.); AliMUONTrackExtrap::AddMCSEffect(¤tParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+ fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.); AliMUONTrackExtrap::AddMCSEffect(¤tParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+ fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.); } // for all planes: extrapolation to the Z of the plane AliMUONTrackExtrap::ExtrapToZCov(¤tParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront()); AliMUONTrackExtrap::ExtrapToZCov(¤tParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront()); AliMUONTrackExtrap::ExtrapToZCov(¤tParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack()); AliMUONTrackExtrap::ExtrapToZCov(¤tParamForResearchBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack()); //--------------------------------------------------------------------------------------- TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances(); TMatrixD covBack(5,5); covBack = currentParamForResearchBack.GetCovariances(); Double_t squaredError_X_Front = covFront(0,0); Double_t squaredError_Y_Front = covFront(2,2); Double_t squaredError_X_Back = covBack(0,0); Double_t squaredError_Y_Back = covBack(2,2); Double_t corrFact = 1.0; Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front); Double_t researchRadiusBack = TMath::Sqrt(squaredError_X_Back + squaredError_Y_Back); if (0.5*(researchRadiusFront+researchRadiusBack)cd(); fHistResearchRadius[planeId] -> Fill(0.5*(researchRadiusFront+researchRadiusBack)); } Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor(); Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor(); Double_t position_X_Back = currentParamForResearchBack.GetNonBendingCoor(); Double_t position_Y_Back = currentParamForResearchBack.GetBendingCoor(); Double_t radialPositionOfTrackFront = TMath::Sqrt(position_X_Front*position_X_Front + position_Y_Front*position_Y_Front); Double_t radialPositionOfTrackBack = TMath::Sqrt(position_X_Back*position_X_Back + position_Y_Back*position_Y_Back); //--------------------------------------------------------------------------------------- Double_t chi2cut = 2.*fScaleSigmaClusterCut*fScaleSigmaClusterCut*fSigmaClusterCut*fSigmaClusterCut; // depends on the number of variables (here, 2) // Analyizing the clusters: FRONT ACTIVE ELEMENTS Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries(); AliDebug(2, Form("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId)); for (Int_t iCluster=0; iClusterAt(iCluster); Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster); // describes the compatibility between the track and the cluster if (chi2GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY()); if (planeId == fNPlanesMFT-1) { if (TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront)cd(); if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf } if (isGoodChi2) { AliDebug(3, Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut)); AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack); if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return 2; newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters())); newTrack->SetPlaneExists(planeId); AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable)); if (fIsCurrentMuonTrackable) { Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2(); AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2())); Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack); // fOutputQAFile->cd(); if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack)); else fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack)); } fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY()); } else AliDebug(3, Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut)); } // Analyizing the clusters: BACK ACTIVE ELEMENTS Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries(); AliDebug(2, Form("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId)); for (Int_t iCluster=0; iClusterAt(iCluster); Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster); // describes the compatibility between the track and the cluster if (chi2GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY()); if (planeId == fNPlanesMFT-1) { if (TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack)cd(); if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf } if (isGoodChi2) { AliDebug(3,Form("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut)); AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack); if (fCandidateTracks->GetEntriesFast() > fMaxNCandidates) return 2; newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster AliDebug(2, Form("After plane %02d: newTrack->GetNMFTClusters() = %d (fCurrentTrack->GetNMFTClusters() = %d)", planeId, newTrack->GetNMFTClusters(), fCurrentTrack->GetNMFTClusters())); newTrack->SetPlaneExists(planeId); AliDebug(2, Form("current muon is trackable: %d\n", fIsCurrentMuonTrackable)); if (fIsCurrentMuonTrackable) { Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2(); AliDebug(2, Form("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2())); Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack); // fOutputQAFile->cd(); if (IsCorrectMatch(cluster)) fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack)); else fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack)); } fGrMFTPlane[kClustersGoodChi2][planeId] -> SetPoint(fGrMFTPlane[kClustersGoodChi2][planeId]->GetN(), cluster->GetX(), cluster->GetY()); } else AliDebug(3,Form("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut)); } //--------------------------------------------------------------------------------------------- if (planeId == fNPlanesMFT-1) { if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) { // fOutputQAFile->cd(); fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane- fDistanceFromGoodClusterAndTrackAtLastPlane)); fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane); } } return 0; } //========================================================================================================================================== void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) { AliDebug(1, Form(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId)); AliMUONTrackParam currentParamFront, currentParamBack; if (planeId == fNPlanesMFT-1) { // last plane of the telecope currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First()))); currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First()))); AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamFront, 0.); AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(¤tParamBack, 0.); } else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed AliDebug(2, Form("fCurrentTrack = %p\n", fCurrentTrack)); currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First()))); currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First()))); AliMUONTrackExtrap::AddMCSEffect(¤tParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+ fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/fRadLengthSi,-1.); AliMUONTrackExtrap::AddMCSEffect(¤tParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+ fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/fRadLengthSi,-1.); } // for all planes: linear extrapolation to the Z of the plane AliMUONTrackExtrap::ExtrapToZCov(¤tParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront()); AliMUONTrackExtrap::ExtrapToZCov(¤tParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack()); Bool_t goodClusterFound = kFALSE; // Analyizing the clusters: FRONT ACTIVE ELEMENTS Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries(); AliDebug(1, Form("nClustersFront = %d\n", nClustersFront)); for (Int_t iCluster=0; iClusterUncheckedAt(iCluster); AliDebug(2, Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack)); if (IsCorrectMatch(cluster)) { fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster fCurrentTrack->SetPlaneExists(planeId); goodClusterFound = kTRUE; break; } } if (goodClusterFound) return; // Analyizing the clusters: BACK ACTIVE ELEMENTS Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries(); AliDebug(1, Form("nClustersBack = %d\n", nClustersBack)); for (Int_t iCluster=0; iClusterUncheckedAt(iCluster); AliDebug(2,Form("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack)); if (IsCorrectMatch(cluster)) { fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster fCurrentTrack->SetPlaneExists(planeId); goodClusterFound = kTRUE; break; } } } //========================================================================================================================================== void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() { for (Int_t iPlane=0; iPlaneGetEntriesFast(); for (Int_t iCluster=0; iClusterAt(iCluster); for (Int_t iTrack=0; iTrackGetNMCTracks(); iTrack++) { if (cluster->GetMCLabel(iTrack)==fLabelMC) { fIsGoodClusterInPlane[iPlane] = kTRUE; break; } } } } fIsCurrentMuonTrackable = kTRUE; for (Int_t iPlane=0; iPlaneGetNMFTClusters() = %d\n", fFinalBestCandidate->GetNMFTClusters())); for (Int_t iPlane=0; iPlanePlaneExists(iPlane)) { AliMFTCluster *trackCluster = fFinalBestCandidate->GetMFTCluster(cluster++); fGrMFTPlane[kClusterOfTrack][iPlane] -> SetPoint(fGrMFTPlane[kClusterOfTrack][iPlane]->GetN(), trackCluster->GetX(), trackCluster->GetY()); } Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast(); for (Int_t iCluster=0; iClusterUncheckedAt(iCluster); fGrMFTPlane[kAllClusters][iPlane] -> SetPoint(fGrMFTPlane[kAllClusters][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY()); if (IsCorrectMatch(myCluster)) { fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetPoint(fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN(), myCluster->GetX(), myCluster->GetY()); } } } } //====================================================================================================================================== Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) { Bool_t result = kFALSE; // check if the cluster belongs to the correct MC track for (Int_t iTrack=0; iTrackGetNMCTracks(); iTrack++) { if (cluster->GetMCLabel(iTrack)==fLabelMC) { result = kTRUE; break; } } AliDebug(2,Form("returning %d\n", result)); return result; } //====================================================================================================================================== Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) { // Test the compatibility between the track and the cluster (using trackParam's covariance matrix): // return the corresponding Chi2 // assume the track parameters are given at the Z of the cluster // Set differences between trackParam and cluster in the bending and non bending directions Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor(); Double_t dY = cluster->GetY() - trackParam.GetBendingCoor(); AliDebug(3,Form("dX = %f, dY = %f\n", dX, dY)); // Calculate errors and covariances const TMatrixD& kParamCov = trackParam.GetCovariances(); Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2(); Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2(); AliDebug(3, Form("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2)); Double_t covXY = kParamCov(0,2); Double_t det = sigmaX2 * sigmaY2 - covXY * covXY; // Compute chi2 if (det==0.) return 1.e10; return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det; } //========================================================================================================================================= void AliMuonForwardTrackFinder::SeparateFrontBackClusters() { for (Int_t iPlane=0; iPlaneDelete(); fMFTClusterArrayBack[iPlane] ->Delete(); for (Int_t iCluster=0; iClusterGetEntries(); iCluster++) { AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster); if (TMath::Abs(cluster->GetZ())GetPlane(iPlane)->GetZCenter())) { new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster); } else { new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster); } } } } //========================================================================================================================================= Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) { // the same definition as in AliMUONTrack is implemented, since here we just add more clusters to the Muon track Int_t ndf = 2 * nClusters - 5; return (ndf > 0) ? ndf : 0; } //============================================================================================================================================ void AliMuonForwardTrackFinder::BookHistos() { const Int_t nMaxNewTracks[] = {150, 200, 250, 600, 1000}; const Double_t radiusPlane[] = {0.010, 0.010, 0.050, 0.5, 1.5}; fHistRadiusEndOfAbsorber = new TH1D("hRadiusEndOfAbsorber", "Track radial distance at the end of the absorber", 1000, 0, 100.); fHistNGoodClustersForFinalTracks = new TH1D("hNGoodClustersForFinalTracks", "Number of Good Clusters per Final Track", 20, -0.25, 9.75); fHistDistanceGoodClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackAtLastPlane", "Distance of MC Good Cluster from Track in last MFT plane", 200, 0., 2.); fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane", "Good Cluster distance from track - Best Cluster distance from track in last MFT plane", 200, 0., 2.); for (Int_t iPlane=0; iPlane Sumw2(); fHistNGoodClustersForFinalTracks -> Sumw2(); fHistDistanceGoodClusterFromTrackAtLastPlane -> Sumw2(); fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Sumw2(); for (Int_t iPlane=0; iPlane Sumw2(); fHistResearchRadius[iPlane] -> Sumw2(); fHistChi2Cluster_GoodCluster[iPlane] -> Sumw2(); fHistChi2Cluster_BadCluster[iPlane] -> Sumw2(); fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2(); fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Sumw2(); } 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"); 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"); } //============================================================================================================================================ void AliMuonForwardTrackFinder::SetTitleHistos() { fHistRadiusEndOfAbsorber -> SetXTitle("R_{abs} [cm]"); fHistNGoodClustersForFinalTracks -> SetXTitle("N_{GoodClusters}"); fHistDistanceGoodClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]"); fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]"); for (Int_t iPlane=0; iPlane SetXTitle("N_{tracks}"); fHistResearchRadius[iPlane] -> SetXTitle("Research Radius [cm]"); fHistChi2Cluster_GoodCluster[iPlane] -> SetXTitle("#chi^{2}/ndf"); fHistChi2Cluster_BadCluster[iPlane] -> SetXTitle("#chi^{2}/ndf"); fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf"); fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf"); } } //=========================================================================================================================================== void AliMuonForwardTrackFinder::BookPlanes() { for (Int_t iPlane=0; iPlane SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane)); fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerStyle(20); // fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.5); // fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.3); fGrMFTPlane[kAllClusters][iPlane] -> SetMarkerSize(0.2); } for (Int_t iPlane=0; iPlane SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane)); fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerStyle(20); // fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.8); // fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.4); fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerSize(0.3); fGrMFTPlane[kClustersGoodChi2][iPlane] -> SetMarkerColor(kBlue); } for (Int_t iPlane=0; iPlane SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane)); fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerStyle(25); // fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(1.2); fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerSize(0.9); fGrMFTPlane[kClusterOfTrack][iPlane] -> SetMarkerColor(kRed); fGrMFTPlane[kClusterOfTrack][iPlane] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane])); } for (Int_t iPlane=0; iPlane SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane)); fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerStyle(20); // fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.8); fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerSize(0.5); fGrMFTPlane[kClusterCorrectMC][iPlane] -> SetMarkerColor(kGreen); } for (Int_t iPlane=0; iPlane SetMarkerSize(0.5); fMrkClustGoodChi2 = new TMarker(0.10, 0.26, 20); fMrkClustGoodChi2 -> SetMarkerSize(0.8); fMrkClustGoodChi2 -> SetMarkerColor(kBlue); fMrkClustMC = new TMarker(0.10, 0.20, 20); fMrkClustMC -> SetMarkerSize(0.8); fMrkClustMC -> SetMarkerColor(kGreen); fMrkClustOfTrack = new TMarker(0.10, 0.14, 25); fMrkClustOfTrack -> SetMarkerSize(1.2); fMrkClustOfTrack -> SetMarkerColor(kRed); fTxtAllClust = new TLatex(0.15, 0.30, "All Clusters"); fTxtAllClust -> SetTextSize(0.040); fTxtClustGoodChi2 = new TLatex(0.15, 0.24, "Clusters involved in the research"); fTxtClustGoodChi2 -> SetTextSize(0.040); fTxtClustMC = new TLatex(0.15, 0.18, "MC good clusters"); fTxtClustMC -> SetTextSize(0.040); fTxtClustOfTrack = new TLatex(0.15, 0.12, "Clusters of the best candidate"); fTxtClustOfTrack -> SetTextSize(0.040); } //=========================================================================================================================================== void AliMuonForwardTrackFinder::ResetPlanes() { for (Int_t iPlane=0; iPlaneGetN(); for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iGr][iPlane]->RemovePoint(iPoint); } } } //=========================================================================================================================================== void AliMuonForwardTrackFinder::PrintParticleHistory() { AliDebug(1, "Entering"); TString history = ""; TParticle *part = 0; if (fLabelMC>=0) part = fStack->Particle(fLabelMC); AliDebug(1, Form("fStack->Particle(%d) = %p", fLabelMC, part)); AliDebug(1, Form("fStack->Particle(%d)->GetPdgCode() = %d", fLabelMC, part->GetPdgCode())); if (part) { if (part->GetFirstMother() != -1) { TParticle *partMother = fStack->Particle(part->GetFirstMother()); AliDebug(1, Form("fStack->Particle(%d) = %p", part->GetFirstMother(), partMother)); if (partMother) { Char_t newName[100]; if (partMother->GetFirstMother() != -1) history += "... #rightarrow "; PDGNameConverter(partMother->GetName(), newName); history += Form("%s #rightarrow ", newName); } } Char_t newName[100]; AliDebug(1, Form("fStack->Particle(%d)->GetPdgCode() = %d", fLabelMC, part->GetPdgCode())); AliDebug(1, Form("fStack->Particle(%d)->GetName() = %s", fLabelMC, part->GetName())); PDGNameConverter(part->GetName(), newName); history += Form("%s at z = %5.1f cm", newName, part->Vz()); // printf("%s", history.Data()); } else history += "NO AVAILABLE HISTORY"; fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data()); // Filling particle history in the fFinalBestCandidate if (part) { for (Int_t iParent=0; iParentGetFirstMother() == -1) break; if (!(fStack->Particle(part->GetFirstMother()))) break; AliDebug(1, Form("fStack->Particle(part->GetFirstMother() = %p", fStack->Particle(part->GetFirstMother()))); fFinalBestCandidate->SetParentMCLabel(iParent, part->GetFirstMother()); fFinalBestCandidate->SetParentPDGCode(iParent, fStack->Particle(part->GetFirstMother())->GetPdgCode()); part = fStack->Particle(part->GetFirstMother()); } } } //=========================================================================================================================================== Bool_t AliMuonForwardTrackFinder::IsMother(const Char_t *nameMother) { Bool_t result = kFALSE; TParticle *part = 0; if (fLabelMC>=0) part = fStack->Particle(fLabelMC); if (part) { if (part->GetFirstMother() != -1) { TParticle *partMother = fStack->Particle(part->GetFirstMother()); if (partMother) { if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE; } } } return result; } //=========================================================================================================================================== void AliMuonForwardTrackFinder::DrawPlanes() { fCanvas -> Clear(); if (fNPlanesMFT <= 5) fCanvas -> Divide(3,2); else if (fNPlanesMFT <= 11) fCanvas -> Divide(4,3); else if (fNPlanesMFT <= 19) fCanvas -> Divide(5,4); for (Int_t iPlane=0; iPlanecd(fNPlanesMFT-iPlane+1); fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]); fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]); fGrMFTPlane[kClusterOfTrack][iPlane] -> GetXaxis() -> SetTitle("X [cm]"); fGrMFTPlane[kClusterOfTrack][iPlane] -> GetYaxis() -> SetTitle("Y [cm]"); fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("ap"); fCircleExt[iPlane] -> Draw("same"); fCircleInt[iPlane] -> Draw("same"); if (fGrMFTPlane[kAllClusters][iPlane]->GetN()) fGrMFTPlane[kAllClusters][iPlane] -> Draw("psame"); if (fGrMFTPlane[kClustersGoodChi2][iPlane]->GetN()) fGrMFTPlane[kClustersGoodChi2][iPlane] -> Draw("psame"); if (fGrMFTPlane[kClusterOfTrack][iPlane]->GetN()) fGrMFTPlane[kClusterOfTrack][iPlane] -> Draw("psame"); if (fGrMFTPlane[kClusterCorrectMC][iPlane]->GetN()) fGrMFTPlane[kClusterCorrectMC][iPlane] -> Draw("psame"); fTxtTrackChi2[iPlane] -> Draw("same"); } fCanvas -> cd(1); fTxtMuonHistory -> Draw(); fTxtDummy -> Draw("same"); if (fMatchingMode==kRealMatching) fTxtTrackGoodClusters -> Draw("same"); fTxtTrackFinalChi2 -> Draw("same"); fTxtTrackMomentum -> Draw("same"); if (fMatchingMode==kRealMatching) fTxtFinalCandidates -> Draw("same"); fMrkAllClust -> Draw("same"); fMrkClustGoodChi2 -> Draw("same"); fMrkClustMC -> Draw("same"); fMrkClustOfTrack -> Draw("same"); fTxtAllClust -> Draw("same"); fTxtClustGoodChi2 -> Draw("same"); fTxtClustMC -> Draw("same"); fTxtClustOfTrack -> Draw("same"); // fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent)); fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent)); if (IsMother("phi")) { fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent)); fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent)); } if (IsMother("J/psi")) { fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent)); fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent)); } } //=========================================================================================================================================== void AliMuonForwardTrackFinder::Terminate() { AliInfo(""); AliInfo("---------------------------------------------------------------------------------------------------------------"); AliInfo(Form("%8d tracks analyzed", fCountRealTracksAnalyzed)); AliInfo(Form("%8d tracks with MC ref", fCountRealTracksWithRefMC)); AliInfo(Form("%8d tracks with MC ref & trigger match", fCountRealTracksWithRefMC_andTrigger)); if (fMatchingMode==kRealMatching) { AliInfo(Form("%8d tracks analyzed with final candidates", fCountRealTracksAnalyzedWithFinalCandidates)); } else { AliInfo(Form("%8d tracks matched with their MC clusters", fCountRealTracksAnalyzedWithFinalCandidates)); } // printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c", fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut); // printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut); AliInfo("---------------------------------------------------------------------------------------------------------------"); WriteOutputTree(); WriteHistos(); } //========================================================================================================================================== void AliMuonForwardTrackFinder::FillOutputTree() { if (!fMuonForwardTracks || !fOutputEventTree) return; AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p", fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0))); // fOutputTreeFile->cd(); fOutputEventTree->Fill(); AliDebug(1, Form("\nFilled Tree: nEvents = %d!!!!\n", Int_t(fOutputEventTree->GetEntries()))); } //========================================================================================================================================== void AliMuonForwardTrackFinder::WriteOutputTree() { if (!fOutputEventTree || !fOutputTreeFile) return; fOutputTreeFile -> cd(); fOutputEventTree -> Write(); fOutputTreeFile -> Close(); } //========================================================================================================================================== void AliMuonForwardTrackFinder::WriteHistos() { fOutputQAFile = new TFile(Form("MuonGlobalTracking.QA.run%d.root", fRun), "recreate"); fOutputQAFile -> cd(); fHistRadiusEndOfAbsorber -> Write(); fHistNGoodClustersForFinalTracks -> Write(); fHistDistanceGoodClusterFromTrackAtLastPlane -> Write(); fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Write(); for (Int_t iPlane=0; iPlane Write(); fHistResearchRadius[iPlane] -> Write(); fHistChi2Cluster_GoodCluster[iPlane] -> Write(); fHistChi2Cluster_BadCluster[iPlane] -> Write(); fHistGlobalChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write(); fHistGlobalChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Write(); } fNtuFinalCandidates -> Write(); fNtuFinalBestCandidates -> Write(); fOutputQAFile -> Close(); } //=========================================================================================================================================== void AliMuonForwardTrackFinder::PDGNameConverter(const Char_t *nameIn, Char_t *nameOut) { if (!strcmp(nameIn, "mu+")) snprintf(nameOut, 50, "#mu^{+}"); else if (!strcmp(nameIn, "mu-")) snprintf(nameOut, 50, "#mu^{-}"); else if (!strcmp(nameIn, "pi+")) snprintf(nameOut, 50, "#pi^{+}"); else if (!strcmp(nameIn, "pi-")) snprintf(nameOut, 50, "#pi^{-}"); else if (!strcmp(nameIn, "K+")) snprintf(nameOut, 50, "K^{+}"); else if (!strcmp(nameIn, "K-")) snprintf(nameOut, 50, "K^{-}"); else if (!strcmp(nameIn, "K*+")) snprintf(nameOut, 50, "K^{*+}"); else if (!strcmp(nameIn, "K*-")) snprintf(nameOut, 50, "K^{*-}"); else if (!strcmp(nameIn, "K_S0")) snprintf(nameOut, 50, "K_{S}^{0}"); else if (!strcmp(nameIn, "K_L0")) snprintf(nameOut, 50, "K_{L}^{0}"); else if (!strcmp(nameIn, "K0")) snprintf(nameOut, 50, "K^{0}"); else if (!strcmp(nameIn, "K0_bar")) snprintf(nameOut, 50, "#bar{K}^{0}"); else if (!strcmp(nameIn, "K*0")) snprintf(nameOut, 50, "K^{*0}"); else if (!strcmp(nameIn, "K*0_bar")) snprintf(nameOut, 50, "#bar{K}^{*0}"); else if (!strcmp(nameIn, "rho0")) snprintf(nameOut, 50, "#rho^{0}"); else if (!strcmp(nameIn, "rho+")) snprintf(nameOut, 50, "#rho^{+}"); else if (!strcmp(nameIn, "rho-")) snprintf(nameOut, 50, "#rho^{-}"); else if (!strcmp(nameIn, "omega")) snprintf(nameOut, 50, "#omega"); else if (!strcmp(nameIn, "eta'")) snprintf(nameOut, 50, "#eta'"); else if (!strcmp(nameIn, "phi")) snprintf(nameOut, 50, "#phi"); else if (!strcmp(nameIn, "D-")) snprintf(nameOut, 50, "D^{-}"); else if (!strcmp(nameIn, "D+")) snprintf(nameOut, 50, "D^{+}"); else if (!strcmp(nameIn, "D0")) snprintf(nameOut, 50, "D^{0}"); else if (!strcmp(nameIn, "D0_bar")) snprintf(nameOut, 50, "#bar{D}^{0}"); else if (!strcmp(nameIn, "D*-")) snprintf(nameOut, 50, "D^{*-}"); else if (!strcmp(nameIn, "D*+")) snprintf(nameOut, 50, "D^{*+}"); else if (!strcmp(nameIn, "D_s+")) snprintf(nameOut, 50, "D_{s}^{+}"); else if (!strcmp(nameIn, "D*_s+")) snprintf(nameOut, 50, "D_{s}^{*+}"); else if (!strcmp(nameIn, "B-")) snprintf(nameOut, 50, "B^{-}"); else if (!strcmp(nameIn, "B+")) snprintf(nameOut, 50, "B^{+}"); else if (!strcmp(nameIn, "B_s0_bar")) snprintf(nameOut, 50, "#bar{B}_{s}^{0}"); else if (!strcmp(nameIn, "antiproton")) snprintf(nameOut, 50, "#bar{p}"); else if (!strcmp(nameIn, "proton")) snprintf(nameOut, 50, "p"); else if (!strcmp(nameIn, "neutron")) snprintf(nameOut, 50, "n"); else if (!strcmp(nameIn, "Sigma+")) snprintf(nameOut, 50, "#Sigma^{+}"); else if (!strcmp(nameIn, "Delta+")) snprintf(nameOut, 50, "#Delta{+}"); else if (!strcmp(nameIn, "Delta--")) snprintf(nameOut, 50, "#Delta{--}"); else if (!strcmp(nameIn, "Lambda0")) snprintf(nameOut, 50, "#Lambda_0"); else if (!strcmp(nameIn, "Lambda0_bar")) snprintf(nameOut, 50, "#bar{Lambda}_0"); else snprintf(nameOut, 50, "%s", nameIn); } //=========================================================================================================================================== void AliMuonForwardTrackFinder::SetDraw(Bool_t drawOption) { fDrawOption = drawOption; if (!fCanvas) { fCanvas = new TCanvas("tracking", "tracking", 1200, 800); fCanvas -> Divide(3,2); } } //=========================================================================================================================================== Bool_t AliMuonForwardTrackFinder::InitGRP() { //------------------------------------ // Initialization of the GRP entry //------------------------------------ AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data"); if (entry) { TMap* m = dynamic_cast(entry->GetObject()); // old GRP entry if (m) { AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject"); m->Print(); fGRPData = new AliGRPObject(); fGRPData->ReadValuesFromMap(m); } else { AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it"); fGRPData = dynamic_cast(entry->GetObject()); // new GRP entry entry->SetOwner(0); } // FIX ME: The unloading of GRP entry is temporarily disabled // because ZDC and VZERO are using it in order to initialize // their reconstructor objects. In the future one has to think // of propagating AliRunInfo to the reconstructors. // AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data"); } if (!fGRPData) { AliError("No GRP entry found in OCDB!"); return kFALSE; } TString lhcState = fGRPData->GetLHCState(); if (lhcState==AliGRPObject::GetInvalidString()) { AliError("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN"); lhcState = "UNKNOWN"; } TString beamType = fGRPData->GetBeamType(); if (beamType==AliGRPObject::GetInvalidString()) { AliError("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN"); beamType = "UNKNOWN"; } Float_t beamEnergy = fGRPData->GetBeamEnergy(); if (beamEnergy==AliGRPObject::GetInvalidFloat()) { AliError("GRP/GRP/Data entry: missing value for the beam energy ! Using 0"); beamEnergy = 0; } TString runType = fGRPData->GetRunType(); if (runType==AliGRPObject::GetInvalidString()) { AliError("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN"); runType = "UNKNOWN"; } Int_t activeDetectors = fGRPData->GetDetectorMask(); if (activeDetectors==AliGRPObject::GetInvalidUInt()) { AliError("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399"); activeDetectors = 1074790399; } AliDebug(1, Form("activeDetectors = %d", activeDetectors)); fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors); fRunInfo->Dump(); // *** Dealing with the magnetic field map if ( TGeoGlobalMagField::Instance()->IsLocked() ) { if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) { AliInfo("ExpertMode!!! GRP information will be ignored !"); AliInfo("ExpertMode!!! Running with the externally locked B field !"); } else { AliInfo("Destroying existing B field instance!"); delete TGeoGlobalMagField::Instance(); } } if ( !TGeoGlobalMagField::Instance()->IsLocked() ) { // Construct the field map out of the information retrieved from GRP. Bool_t ok = kTRUE; // L3 Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0); if (l3Current == AliGRPObject::GetInvalidFloat()) { AliError("GRP/GRP/Data entry: missing value for the L3 current !"); ok = kFALSE; } Char_t l3Polarity = fGRPData->GetL3Polarity(); if (l3Polarity == AliGRPObject::GetInvalidChar()) { AliError("GRP/GRP/Data entry: missing value for the L3 polarity !"); ok = kFALSE; } // Dipole Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0); if (diCurrent == AliGRPObject::GetInvalidFloat()) { AliError("GRP/GRP/Data entry: missing value for the dipole current !"); ok = kFALSE; } Char_t diPolarity = fGRPData->GetDipolePolarity(); if (diPolarity == AliGRPObject::GetInvalidChar()) { AliError("GRP/GRP/Data entry: missing value for the dipole polarity !"); ok = kFALSE; } // read special bits for the polarity convention and map type Int_t polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008; Bool_t uniformB = fGRPData->IsUniformBMap(); if (ok) { AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1), TMath::Abs(diCurrent) * (diPolarity ? -1:1), polConvention,uniformB,beamEnergy, beamType.Data()); if (fld) { TGeoGlobalMagField::Instance()->SetField( fld ); TGeoGlobalMagField::Instance()->Lock(); AliInfo("Running with the B field constructed out of GRP !"); } else AliFatal("Failed to create a B field map !"); } else AliFatal("B field is neither set nor constructed from GRP ! Exitig..."); } return kTRUE; } //==================================================================================================================================================== Bool_t AliMuonForwardTrackFinder::SetRunNumber() { AliCDBManager *man = AliCDBManager::Instance(); if (!fRunLoader) { AliError("No run loader found!"); return kFALSE; } else { fRunLoader->LoadHeader(); // read run number from gAlice if (fRunLoader->GetHeader()) { man->SetRun(fRunLoader->GetHeader()->GetRun()); fRunLoader->UnloadHeader(); } else { AliError("No run-loader header found!"); return kFALSE; } } return kTRUE; } //====================================================================================================================================================