From f202486b0582208867a62e459a8b11fa19cfaef2 Mon Sep 17 00:00:00 2001 From: hristov Date: Mon, 2 Nov 2009 12:51:07 +0000 Subject: [PATCH] Fixes in reconstruction: 1) Define in a single place (AliMUONTrack) the algorithm to be used to associate 2 tracks by comparing the position of their clusters. Update all the classes/macro using that algorithm to find the simulated track corresponding to a given reconstructed one: AliMUONTrack: - Change the function IsValid(...) to be able to check the reconstructibility of the track depending on the tracking algorithm - Change the name of the function CompatibleTrack to FindCompatibleCluster that is actually what the function does - New function Match(...) to associate, or not, 2 tracks by comparing the position of their clusters (can be used to assign MC labels to reconstructed tracks). AliMUONRecoCheck: - Check the track reconstructibility according to given informations about the tracking algorithm (those info can be extracted from RecoParam). - New function FindCompatibleTrack(...) to find in the given store a track that match the given track. Matching is done by using MC labels or by comparing clusters' position. All other classes/macros/tasks: - use these new functionnalities 2) Additional correction: DecodeRecoCocktail.C: - correct the algorithm to make it working even if there are ghost tracks - avoid accessing OCDB when not needed MUONRecoCheck.C: - new histograms/graphs to study the momentum resolution as a function of the momentum - new histograms/graphs to study the cluster resolution at each chamber - avoid loading magnetic field when not needed AliMUONTrackLight: - correct constructor - avoid using magnetic field when not needed (Philippe P.) --- MUON/AliMUONRecoCheck.cxx | 80 ++++-- MUON/AliMUONRecoCheck.h | 9 +- MUON/AliMUONTrack.cxx | 56 +++- MUON/AliMUONTrack.h | 5 +- MUON/AliMUONTrackLight.cxx | 110 ++------ MUON/AliMUONTrackLight.h | 2 - MUON/DIMUONFakes.C | 74 +----- MUON/DecodeRecoCocktail.C | 83 +++--- MUON/MUONFakes.C | 132 ++-------- MUON/MUONRecoCheck.C | 241 ++++++++++-------- .../AliAnalysisTaskESDMCLabelAddition.cxx | 114 +-------- PWG3/muondep/AliAnalysisTaskRecoCheck.cxx | 68 +++-- 12 files changed, 368 insertions(+), 606 deletions(-) diff --git a/MUON/AliMUONRecoCheck.cxx b/MUON/AliMUONRecoCheck.cxx index 0daff813097..aaef0b0e9cc 100644 --- a/MUON/AliMUONRecoCheck.cxx +++ b/MUON/AliMUONRecoCheck.cxx @@ -220,14 +220,15 @@ AliMUONVTrackStore* AliMUONRecoCheck::TrackRefs(Int_t event) } //_____________________________________________________________________________ -AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event) +AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event, UInt_t requestedStationMask, Bool_t request2ChInSameSt45) { - /// Return a track store containing the reconstructible tracks for a given event + /// Return a track store containing the reconstructible tracks for a given event, + /// according to the mask of requested stations and the minimum number of chambers hit in stations 4 & 5. if (!fESDEventOwner) { if (fRecoTrackRefStore == 0x0) { if (TrackRefs(event) == 0x0) return 0x0; - MakeReconstructibleTracks(); + MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45); } return fRecoTrackRefStore; } @@ -240,7 +241,7 @@ AliMUONVTrackStore* AliMUONRecoCheck::ReconstructibleTracks(Int_t event) if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore; else { if (TrackRefs(event) == 0x0) return 0x0; - MakeReconstructibleTracks(); + MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45); return fRecoTrackRefStore; } } @@ -502,7 +503,7 @@ void AliMUONRecoCheck::CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefSt } //_____________________________________________________________________________ -void AliMUONRecoCheck::MakeReconstructibleTracks() +void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45) { /// Isolate the reconstructible tracks if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return; @@ -510,36 +511,59 @@ void AliMUONRecoCheck::MakeReconstructibleTracks() // create iterator on trackRef TIter next(fTrackRefStore->CreateIterator()); - // loop over trackRef + // loop over trackRef and add reconstructible tracks to fRecoTrackRefStore AliMUONTrack* track; while ( ( track = static_cast(next()) ) ) { + if (track->IsValid(requestedStationMask, request2ChInSameSt45)) fRecoTrackRefStore->Add(*track); + } + +} + +//_____________________________________________________________________________ +AliMUONTrack* AliMUONRecoCheck::FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore, + Int_t &nMatchClusters, Bool_t useLabel, Double_t sigmaCut) +{ + /// Return the track from the store matched with the given track (or 0x0) and the number of matched clusters. + /// Matching is done by using the MC label of by comparing cluster/TrackRef positions according to the flag "useLabel". + /// WARNING: Who match who matters since the matching algorithm uses the *fraction* of matched clusters of the given track + + AliMUONTrack *matchedTrack = 0x0; + nMatchClusters = 0; + + if (useLabel) { // by using the MC label - Bool_t* chamberInTrack = new Bool_t(AliMUONConstants::NTrackingCh()); - for (Int_t iCh = 0; iCh < AliMUONConstants::NTrackingCh(); iCh++) chamberInTrack[iCh] = kFALSE; + // get the corresponding simulated track if any + Int_t label = track.GetMCLabel(); + matchedTrack = (AliMUONTrack*) trackStore.FindObject(label); - // loop over trackRef's hits to get hit chambers - Int_t nTrackHits = track->GetNClusters(); - for (Int_t iHit = 0; iHit < nTrackHits; iHit++) { - AliMUONVCluster* hit = ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->UncheckedAt(iHit))->GetClusterPtr(); - chamberInTrack[hit->GetChamberId()] = kTRUE; - } + // get the fraction of matched clusters + if (matchedTrack) { + Int_t nClusters = track.GetNClusters(); + for (Int_t iCl = 0; iCl < nClusters; iCl++) + if (((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label) + nMatchClusters++; + } - // track is reconstructible if the particle is depositing a hit - // in the following chamber combinations: - Bool_t trackOK = kTRUE; - if (!chamberInTrack[0] && !chamberInTrack[1]) trackOK = kFALSE; - if (!chamberInTrack[2] && !chamberInTrack[3]) trackOK = kFALSE; - if (!chamberInTrack[4] && !chamberInTrack[5]) trackOK = kFALSE; - Int_t nHitsInLastStations = 0; - for (Int_t iCh = 6; iCh < AliMUONConstants::NTrackingCh(); iCh++) - if (chamberInTrack[iCh]) nHitsInLastStations++; - if(nHitsInLastStations < 3) trackOK = kFALSE; + } else { // by comparing cluster/TrackRef positions - // Add reconstructible tracks to fRecoTrackRefStore - if (trackOK) fRecoTrackRefStore->Add(*track); + // look for the corresponding simulated track if any + TIter next(trackStore.CreateIterator()); + AliMUONTrack* track2; + while ( ( track2 = static_cast(next()) ) ) { + + // check compatibility + Int_t n = 0; + if (track.Match(*track2, sigmaCut, n)) { + matchedTrack = track2; + nMatchClusters = n; + break; + } + + } - delete [] chamberInTrack; } - + + return matchedTrack; + } diff --git a/MUON/AliMUONRecoCheck.h b/MUON/AliMUONRecoCheck.h index 0897c1b2851..1d0aa412fde 100644 --- a/MUON/AliMUONRecoCheck.h +++ b/MUON/AliMUONRecoCheck.h @@ -18,6 +18,7 @@ class TTree; class AliESDEvent; class AliMCEventHandler; class AliMUONVTrackStore; +class AliMUONTrack; class AliMUONRecoCheck : public TObject { @@ -33,7 +34,7 @@ public: AliMUONVTrackStore* TrackRefs(Int_t event); /// Return reconstructible reference tracks - AliMUONVTrackStore* ReconstructibleTracks(Int_t event); + AliMUONVTrackStore* ReconstructibleTracks(Int_t event, UInt_t requestedStationMask = 0x1F, Bool_t request2ChInSameSt45 = kTRUE); /// Return the run number of the current ESD event Int_t GetRunNumber(); @@ -47,6 +48,10 @@ public: /// Return the interface to the Monte Carlo data of current event const AliMCEventHandler* GetMCEventHandler() { return fMCEventHandler; } + /// Return the track from the store matched with the given track (or 0x0) and the fraction of matched clusters. + static AliMUONTrack* FindCompatibleTrack(AliMUONTrack &track, AliMUONVTrackStore &trackStore, + Int_t &nMatchClusters, Bool_t useLabel = kFALSE, Double_t sigmaCut = 10.); + private: /// Not implemented AliMUONRecoCheck(const AliMUONRecoCheck& rhs); @@ -61,7 +66,7 @@ private: void CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore); - void MakeReconstructibleTracks(); + void MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45 = kTRUE); private: AliMCEventHandler* fMCEventHandler; ///< to access MC truth information diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index 8bcc1c9bb43..d402a983df3 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -472,15 +472,16 @@ Bool_t AliMUONTrack::UpdateCovTrackParamAtCluster() } //__________________________________________________________________________ -Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask) +Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask, Bool_t request2ChInSameSt45) { /// check the validity of the current track: /// at least one cluster per requested station /// and at least 2 chambers in stations 4 & 5 that contain cluster(s) + /// + if request2ChInSameSt45 = kTRUE: 2 chambers hit in the same station (4 or 5) Int_t nClusters = GetNClusters(); AliMUONTrackParam *trackParam; - Int_t currentCh, currentSt, previousCh = -1, nChHitInSt45 = 0; + Int_t currentCh, currentSt, previousCh = -1, nChHitInSt4 = 0, nChHitInSt5 = 0; UInt_t presentStationMask(0); // first loop over clusters @@ -493,15 +494,28 @@ Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask) // build present station mask presentStationMask |= ( 1 << currentSt ); - // count the number of chambers in station 4 & 5 that contain cluster(s) - if (currentCh > 5 && currentCh != previousCh) { - nChHitInSt45++; + // count the number of chambers hit in station 4 that contain cluster(s) + if (currentSt == 3 && currentCh != previousCh) { + nChHitInSt4++; + previousCh = currentCh; + } + + // count the number of chambers hit in station 5 that contain cluster(s) + if (currentSt == 4 && currentCh != previousCh) { + nChHitInSt5++; previousCh = currentCh; } } - return (((requestedStationMask & presentStationMask) == requestedStationMask) && (nChHitInSt45 >= 2)); + // at least one cluster per requested station + if ((requestedStationMask & presentStationMask) != requestedStationMask) return kFALSE; + + // 2 chambers hit in the same station (4 or 5) + if (request2ChInSameSt45) return (nChHitInSt4 == 2 || nChHitInSt5 == 2); + // or 2 chambers hit in station 4 & 5 together + else return (nChHitInSt4+nChHitInSt5 >= 2); + } //__________________________________________________________________________ @@ -1066,10 +1080,11 @@ Double_t AliMUONTrack::GetNormalizedChi2() const } //__________________________________________________________________________ -Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const +Int_t AliMUONTrack::FindCompatibleClusters(AliMUONTrack &track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const { - /// for each chamber: return kTRUE (kFALSE) if clusters are compatible (not compatible). - /// nMatchClusters = number of clusters of "this" track matched with one cluster of track "track" + /// Try to match clusters from this track with clusters from the given track within the provided sigma cut: + /// - Fill the array compatibleCluster[iCh] with kTRUE if a compatible cluster has been found in chamber iCh. + /// - Return the number of clusters of "this" track matched with one cluster of the given track. AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2; AliMUONVCluster *cluster1, *cluster2; Double_t chi2, dX, dY; @@ -1079,7 +1094,7 @@ Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool Int_t nMatchClusters = 0; for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE; - if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters; + if (!track.fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters; // Loop over clusters of first track trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First(); @@ -1088,13 +1103,13 @@ Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool cluster1 = trackParamAtCluster1->GetClusterPtr(); // Loop over clusters of second track - trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First(); + trackParamAtCluster2 = (AliMUONTrackParam*) track.fTrackParamAtCluster->First(); while (trackParamAtCluster2) { cluster2 = trackParamAtCluster2->GetClusterPtr(); //prepare next step - trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2); + trackParamAtCluster2 = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster2); // check DE Id if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue; @@ -1116,6 +1131,23 @@ Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool return nMatchClusters; } +//__________________________________________________________________________ +Bool_t AliMUONTrack::Match(AliMUONTrack &track, Double_t sigmaCut, Int_t &nMatchClusters) const +{ + /// Try to match this track with the given track. Matching conditions: + /// - more than 50% of clusters from this track matched with clusters from the given track + /// - at least 1 cluster matched before and 1 cluster matched after the dipole + + Bool_t compTrack[10]; + nMatchClusters = FindCompatibleClusters(track, sigmaCut, compTrack); + + if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2 + (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5 + 2 * nMatchClusters > GetNClusters()) return kTRUE; // more than 50% of clusters matched + else return kFALSE; + +} + //__________________________________________________________________________ void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam) { diff --git a/MUON/AliMUONTrack.h b/MUON/AliMUONTrack.h index 943f52a73eb..e6f5a5dac62 100644 --- a/MUON/AliMUONTrack.h +++ b/MUON/AliMUONTrack.h @@ -38,7 +38,7 @@ class AliMUONTrack : public TObject Bool_t UpdateTrackParamAtCluster(); Bool_t UpdateCovTrackParamAtCluster(); - Bool_t IsValid(UInt_t requestedStationMask); + Bool_t IsValid(UInt_t requestedStationMask, Bool_t request2ChInSameSt45 = kFALSE); void TagRemovableClusters(UInt_t requestedStationMask); @@ -94,7 +94,8 @@ class AliMUONTrack : public TObject Int_t GetNDF() const; Double_t GetNormalizedChi2() const; - Int_t CompatibleTrack(AliMUONTrack* track, Double_t sigma2Cut, Bool_t compatibleCluster[10]) const; + Int_t FindCompatibleClusters(AliMUONTrack &track, Double_t sigma2Cut, Bool_t compatibleCluster[10]) const; + Bool_t Match(AliMUONTrack &track, Double_t sigma2Cut, Int_t &nMatchClusters) const; /// return pointer to track parameters at vertex (can be 0x0) AliMUONTrackParam* GetTrackParamAtVertex() const {return fTrackParamAtVertex;} diff --git a/MUON/AliMUONTrackLight.cxx b/MUON/AliMUONTrackLight.cxx index 2ac944899b6..ff842c1aa79 100644 --- a/MUON/AliMUONTrackLight.cxx +++ b/MUON/AliMUONTrackLight.cxx @@ -40,13 +40,11 @@ #include "AliMUONTrack.h" #include "AliMUONConstants.h" #include "AliMUONVTrackStore.h" -#include "AliMUONTrackExtrap.h" #include "AliMUONTrackParam.h" #include "AliESDMuonTrack.h" -#include "AliRunLoader.h" #include "AliStack.h" -#include "AliHeader.h" +#include "AliLog.h" #include "TDatabasePDG.h" #include "TParticle.h" @@ -127,8 +125,16 @@ AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack) fWeight(1) { /// constructor - //AliMUONTrackLight(); - FillFromESD(muonTrack); + fPgen.SetPxPyPzE(0.,0.,0.,0.); + for (Int_t npar = 0; npar < fgkNParentsMax; npar++){ + fParentPDGCode[npar] = -1; + fParentPythiaLine[npar] = -1; + } + for (Int_t i = 0; i < 4; i++){ + fQuarkPDGCode[i] = -1; + fQuarkPythiaLine[i] = -1; + } + FillFromESD(muonTrack); } //============================================ @@ -141,19 +147,18 @@ AliMUONTrackLight::~AliMUONTrackLight() void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){ /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack - AliMUONTrackParam trPar; - if (trackReco->GetTrackParamAtVertex()) trPar = *(trackReco->GetTrackParamAtVertex()); - else { - trPar = *((AliMUONTrackParam*) trackReco->GetTrackParamAtCluster()->First()); - AliMUONTrackExtrap::ExtrapToVertex(&trPar,0.,0.,0.,0.,0.); + AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex(); + if (!trPar) { + AliError("The track must contain the parameters at vertex"); + return; } - this->SetCharge(Int_t(TMath::Sign(1.,trPar.GetInverseBendingMomentum()))); - this->SetPxPyPz(trPar.Px(),trPar.Py(), trPar.Pz()); + this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum()))); + this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz()); this->SetTriggered(trackReco->GetMatchTrigger()); - Double_t xyz[3] = { trPar.GetNonBendingCoor(), - trPar.GetBendingCoor(), - trPar.GetZ()}; + Double_t xyz[3] = { trPar->GetNonBendingCoor(), + trPar->GetBendingCoor(), + trPar->GetZ()}; if (zvert!=-9999) xyz[2] = zvert; this->SetVertex(xyz); } @@ -191,81 +196,6 @@ void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){ fPrec.SetPxPyPzE(px,py,pz,energy); } -//============================================ -TParticle* AliMUONTrackLight::FindRefTrack( - AliMUONTrack* trackReco, AliMUONVTrackStore* trackRefArray, - AliStack* stack) -{ - /// Find the Monte Carlo (MC) particle that corresponds to a given reconstructed track. - /// @param trackReco This is the reconstructed track for which we want to find a - /// corresponding MC particle. - /// @param trackRefArray The list of MC reference tracks as generated by - /// AliMUONRecoCheck::ReconstructedTracks for example. - /// @param stack The MC stack of simulated particles. - - TParticle *part = 0; - const Double_t kSigmaCut = 4; // 4 sigmas cut - Int_t compPart = 0; - Bool_t compTrack; - - TIter next(trackRefArray->CreateIterator()); - AliMUONTrack* trackRef; - while ( (trackRef = static_cast(next())) ) { - // check if trackRef is compatible with trackReco: - compTrack = kFALSE; - - if (trackReco->GetNClusters() > 1) { - - // check cluster by cluster if trackReco contain info at each cluster - Bool_t compTrackArray[10]; - trackRef->CompatibleTrack(trackReco,kSigmaCut,compTrackArray); - if (TrackCheck(compTrackArray) == 4) compTrack = kTRUE; - - } else { - - // otherwise check only parameters at the z position of the first trackRef - AliMUONTrackParam *refParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First(); - AliMUONTrackParam recoParam(*((AliMUONTrackParam*) trackReco->GetTrackParamAtCluster()->First())); - Double_t chi2; - if (AliMUONTrackExtrap::ExtrapToZCov(&recoParam, refParam->GetZ()) && - refParam->CompatibleTrackParam(recoParam, kSigmaCut, chi2)) compTrack = kTRUE; - - } - - if (compTrack) { - compPart++; - Int_t trackID = trackRef->GetUniqueID(); - this->SetTrackPythiaLine(trackID); - part = stack->Particle(trackID); - fTrackPDGCode = part->GetPdgCode(); - } - } - if (compPart>1) { - AliFatal("More than one particle compatible with the reconstructed track."); - } - return part; -} - -//============================================ -Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){ - /// Apply reconstruction requirements - /// Return number of validated conditions - /// If all the tests are verified then TrackCheck = 4 (good track) - - Int_t iTrack = 0; - Int_t hitsInLastStations = 0; - - // apply reconstruction requirements - if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0 - if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1 - if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2 - for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) { - if (compTrack[ch]) hitsInLastStations++; - } - if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4 - return iTrack; -} - //============================================ void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){ /// scans the muon history to determine parents pdg code and pythia line diff --git a/MUON/AliMUONTrackLight.h b/MUON/AliMUONTrackLight.h index fe7a49b2288..73c982cb018 100644 --- a/MUON/AliMUONTrackLight.h +++ b/MUON/AliMUONTrackLight.h @@ -95,8 +95,6 @@ class AliMUONTrackLight : public TObject { void SetTriggered(Bool_t isTriggered) { fIsTriggered = isTriggered; } /// Return flag for trigger Bool_t IsTriggered() const { return fIsTriggered; } - TParticle* FindRefTrack(AliMUONTrack* trackReco, AliMUONVTrackStore* trackRefArray, AliStack* stack); - Int_t TrackCheck(Bool_t *compTrack); /// Return acually filled no. of *fragmented* parents Int_t GetNParents() const {return fNParents;} void FillMuonHistory(AliStack *stack, TParticle *part); diff --git a/MUON/DIMUONFakes.C b/MUON/DIMUONFakes.C index 95b7ace37e8..1c15a65d3cf 100644 --- a/MUON/DIMUONFakes.C +++ b/MUON/DIMUONFakes.C @@ -33,9 +33,7 @@ /// Results are saved in the root file DiFakes.root /// Results are relevent provided that you use the same recoParams as for the reconstruction -Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut); -AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore, - Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut); +Double_t sigmaCut = -1.; //----------------------------------------------------------------------- void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1, @@ -78,7 +76,7 @@ void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent if (!recoParam) return; // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used - Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); + sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); TLorentzVector vMu1, vMu2, vDiMu; @@ -106,8 +104,8 @@ void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent AliMUONTrack* muonTrack = (AliMUONTrack*) muonTrackStore->FindObject(esdTrack->GetUniqueID()); // try to match the reconstructed track with a simulated one - Float_t fractionOfMatchCluster = 0.; - AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut); + Int_t nMatchClusters = 0; + AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut); // take actions according to matching result if (matchedTrackRef) { @@ -270,67 +268,3 @@ void DIMUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent } -//----------------------------------------------------------------------- -Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut) -{ - /// Try to match 2 tracks - - Bool_t compTrack[10]; - Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack); - fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters()); - - if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2 - (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5 - fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched - else return kFALSE; - -} - -//----------------------------------------------------------------------- -AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore, - Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut) -{ - /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters - - AliMUONTrack *matchedTrackRef = 0x0; - fractionOfMatchCluster = 0.; - - if (useLabel) { // by using the MC label - - // get the corresponding simulated track if any - Int_t label = muonTrack.GetMCLabel(); - matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label); - - // get the fraction of matched clusters - if (matchedTrackRef) { - Int_t nMatchClusters = 0; - Int_t nClusters = muonTrack.GetNClusters(); - for (Int_t iCl = 0; iCl < nClusters; iCl++) - if (((AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label) - nMatchClusters++; - fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)nClusters); - } - - } else { // by comparing cluster/TrackRef positions - - // look for the corresponding simulated track if any - TIter next(trackRefStore.CreateIterator()); - AliMUONTrack* trackRef; - while ( ( trackRef = static_cast(next()) ) ) { - - // check compatibility - Float_t f = 0.; - if (TrackMatched(muonTrack, *trackRef, f, sigmaCut)) { - matchedTrackRef = trackRef; - fractionOfMatchCluster = f; - break; - } - - } - - } - - return matchedTrackRef; - -} - diff --git a/MUON/DecodeRecoCocktail.C b/MUON/DecodeRecoCocktail.C index 6b2b7df8a81..118f715dfb5 100644 --- a/MUON/DecodeRecoCocktail.C +++ b/MUON/DecodeRecoCocktail.C @@ -45,20 +45,17 @@ #include #include #include -#include -#include "AliMUONCDB.h" #include "AliMUONRecoCheck.h" #include "AliMUONTrack.h" #include "AliMUONTrackLight.h" #include "AliMUONPairLight.h" #include "AliMUONVTrackStore.h" -#include "AliMUONTrackExtrap.h" #include "AliESDEvent.h" #include "AliESDVertex.h" +#include "AliESDMuonTrack.h" #include "AliMCEventHandler.h" #include "AliMCEvent.h" #include "AliStack.h" -#include "AliCDBManager.h" /*TODO: need to update this with changes made to ITS #include "AliITSVertexerPPZ.h" #include "AliITSLoader.h" @@ -69,8 +66,7 @@ void DecodeRecoCocktail( char* recodir=".", // The directory containing galice.root for reconstructed data. char* simdir="generated/", // The directory containing galice.root for simulated data. - char* outFileName = "MuonLight.root", // The output filename containing AliMUONTrackLight and AliMUONPairLight objects. - char* geoFilename = "geometry.root" // The filename containing the geometry. + char* outFileName = "MuonLight.root" // The output filename containing AliMUONTrackLight and AliMUONPairLight objects. ) { /// \param recodir The directory containing galice.root for reconstructed data. @@ -92,24 +88,7 @@ void DecodeRecoCocktail( treeOut->Branch("muons",&muonArray); treeOut->Branch("dimuons",&dimuonArray); - // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex) - if (!gGeoManager) { - TGeoManager::Import(geoFilename); - if (!gGeoManager) { - Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename); - return; - } - } - - // load necessary data from OCDB - AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"); - AliCDBManager::Instance()->SetRun(0); - if (!AliMUONCDB::LoadField()) return; - // set the magnetic field for track extrapolations - AliMUONTrackExtrap::SetField(); - AliMUONRecoCheck *rc = new AliMUONRecoCheck("AliESDs.root", simdir); - Int_t nev = rc->NumberOfEvents(); /*TODO: need to update this with changes made to ITS AliITSLoader* ITSloader = (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader"); @@ -125,6 +104,7 @@ void DecodeRecoCocktail( TLorentzVector v; + Int_t nev = rc->NumberOfEvents(); for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events /*TODO: need to update this with changes made to ITS @@ -143,36 +123,45 @@ void DecodeRecoCocktail( AliMUONVTrackStore* trackRefs = rc->TrackRefs(ievent); AliStack* pstack = (const_cast(rc->GetMCEventHandler()))->MCEvent()->Stack(); - TIter next(recoTracks->CreateIterator()); - AliMUONTrack* trackReco = NULL; - - Int_t nTrackReco = recoTracks->GetSize(); - Int_t nTracksESD = rc->GetESDEvent()->GetNumberOfMuonTracks(); - if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD); + // loop over ESD tracks Int_t nreftracks = 0; - Int_t itrRec = 0; - while ( (trackReco = static_cast(next())) != NULL ) - { - // assign parameters concerning the reconstructed tracks - AliMUONTrackLight muLight; + const AliESDEvent* esd = rc->GetESDEvent(); + Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ; + for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { + + AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack); + + // skip ghosts + if (!esdTrack->ContainTrackerData()) continue; - muLight.FillFromESD(rc->GetESDEvent()->GetMuonTrack(itrRec)); - // muLight.FillFromAliMUONTrack(trackReco); + // find the corresponding MUON track + AliMUONTrack* trackReco = (AliMUONTrack*) recoTracks->FindObject(esdTrack->GetUniqueID()); - // find the reference track and store further information - TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack); - if (part) { + // try to match the reconstructed track with a simulated one + Int_t nMatchClusters = 0; + AliMUONTrack* matchedTrackRef = rc->FindCompatibleTrack(*trackReco, *trackRefs, nMatchClusters, kFALSE, 10.); + + if (matchedTrackRef) { + + //store new referenced track in the muonArray + AliMUONTrackLight* muLight = new ((*muonArray)[nreftracks++]) AliMUONTrackLight(); + + // assign parameters concerning the reconstructed tracks + muLight->FillFromESD(esdTrack); + // muLight->FillFromAliMUONTrack(trackReco); + + // store further information related to the simulated track + muLight->SetTrackPythiaLine(matchedTrackRef->GetUniqueID()); + TParticle *part = pstack->Particle(matchedTrackRef->GetUniqueID()); + muLight->SetTrackPDGCode(part->GetPdgCode()); v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy()); - muLight.SetPGen(v); - muLight.FillMuonHistory(pstack, part); + muLight->SetPGen(v); + muLight->FillMuonHistory(pstack, part); // muLight.PrintInfo("A"); - //store the referenced track in the muonArray: - TClonesArray &muons = *muonArray; - new (muons[nreftracks++]) AliMUONTrackLight(muLight); - } + + } - itrRec++; - } // end reco track + } // end esd track // now loop over muon pairs to build dimuons Int_t nmuons = muonArray->GetEntriesFast(); diff --git a/MUON/MUONFakes.C b/MUON/MUONFakes.C index 60ff78f119b..271938d0fab 100644 --- a/MUON/MUONFakes.C +++ b/MUON/MUONFakes.C @@ -34,12 +34,12 @@ /// Results are saved in the root file Fakes.root /// Results are relevent provided that you use the same recoParams as for the reconstruction -Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut); -Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam); -AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore, - Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut); -Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam, - Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters); +UInt_t requestedStationMask = 0; +Bool_t request2ChInSameSt45 = kFALSE; +Double_t sigmaCut = -1.; + +Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, + Bool_t useLabel, TH1F &hFractionOfConnectedClusters); //----------------------------------------------------------------------- void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = -1, @@ -89,7 +89,11 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = if (!recoParam) return; // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used - Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); + sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); + // compute the mask of requested stations from recoParam + for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i ); + // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible + request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates(); // initialize global counters Int_t nReconstructibleTracks = 0; @@ -121,7 +125,7 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = TIter next(trackRefStore->CreateIterator()); AliMUONTrack* trackRef; while ( ( trackRef = static_cast(next()) ) ) { - if (IsRecontructible(*trackRef,*recoParam)) nReconstructibleTracks++; + if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++; } // loop over ESD tracks @@ -158,21 +162,21 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = hPhi->Fill(phi); // try to match the reconstructed track with a simulated one - Float_t fractionOfMatchCluster = 0.; - AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*muonTrack, *trackRefStore, fractionOfMatchCluster, useLabel, sigmaCut); + Int_t nMatchClusters = 0; + AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(*muonTrack, *trackRefStore, nMatchClusters, useLabel, sigmaCut); // take actions according to matching result if (matchedTrackRef) { // global counter nTotMatchedTracks++; - if (!IsRecontructible(*matchedTrackRef,*recoParam)) { + if (!matchedTrackRef->IsValid(requestedStationMask, request2ChInSameSt45)) { trackReconstructedYet = kTRUE; nTotTracksReconstructedYet++; } // fill histograms - hFractionOfMatchedClusters->Fill(fractionOfMatchCluster); + hFractionOfMatchedClusters->Fill(((Float_t) nMatchClusters) / ((Float_t) nClusters)); hNumberOfClustersMC->Fill(matchedTrackRef->GetNClusters()); hNumberOfClustersM->Fill(nClusters); hChi2PerDofM->Fill(normalizedChi2); @@ -213,8 +217,7 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = if (fakeTrackStore->GetSize() > 0) { // remove the most connected fake tracks - Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, *recoParam, - useLabel, sigmaCut, *hFractionOfConnectedClusters); + Int_t nFreeMissingTracks = RemoveConnectedFakes(*fakeTrackStore, *trackRefStore, useLabel, *hFractionOfConnectedClusters); // remove the remaining free reconstructible tracks Int_t nAdditionalTracks = fakeTrackStore->GetSize() - nFreeMissingTracks; @@ -347,101 +350,8 @@ void MUONFakes(Bool_t useLabel = kFALSE, Int_t FirstEvent = 0, Int_t LastEvent = } //----------------------------------------------------------------------- -Bool_t TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef, Float_t &fractionOfMatchCluster, Double_t sigmaCut) -{ - /// Try to match 2 tracks - - Bool_t compTrack[10]; - Int_t nMatchClusters = track.CompatibleTrack(&trackRef, sigmaCut, compTrack); - fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)track.GetNClusters()); - - if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2 - (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5 - fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched - else return kFALSE; - -} - -//----------------------------------------------------------------------- -Bool_t IsRecontructible(AliMUONTrack &track, AliMUONRecoParam &recoParam) -{ - /// Check il the track is reconstructible - Int_t nMinChHitInSt45 = (recoParam.MakeMoreTrackCandidates()) ? 2 : 3; - Int_t currentCh, previousCh = -1, nChHitInSt45 = 0; - Bool_t clusterInSt[5]; - for (Int_t iSt = 0; iSt < 5; iSt++) clusterInSt[iSt] = !recoParam.RequestStation(iSt); - - AliMUONTrackParam* trackParam = static_cast(track.GetTrackParamAtCluster()->First()); - while (trackParam) { - - currentCh = trackParam->GetClusterPtr()->GetChamberId(); - - clusterInSt[currentCh/2] = kTRUE; - - if (currentCh > 5 && currentCh != previousCh) { - nChHitInSt45++; - previousCh = currentCh; - } - - trackParam = static_cast(track.GetTrackParamAtCluster()->After(trackParam)); - } - - return (clusterInSt[0] && clusterInSt[1] && clusterInSt[2] && - clusterInSt[3] && clusterInSt[4] && nChHitInSt45 >= nMinChHitInSt45); - -} - -//----------------------------------------------------------------------- -AliMUONTrack* MatchWithTrackRef(AliMUONTrack &muonTrack, AliMUONVTrackStore &trackRefStore, - Float_t &fractionOfMatchCluster, Bool_t useLabel, Double_t sigmaCut) -{ - /// Return if the trackRef matched with the reconstructed track and the fraction of matched clusters - - AliMUONTrack *matchedTrackRef = 0x0; - fractionOfMatchCluster = 0.; - - if (useLabel) { // by using the MC label - - // get the corresponding simulated track if any - Int_t label = muonTrack.GetMCLabel(); - matchedTrackRef = (AliMUONTrack*) trackRefStore.FindObject(label); - - // get the fraction of matched clusters - if (matchedTrackRef) { - Int_t nMatchClusters = 0; - Int_t nClusters = muonTrack.GetNClusters(); - for (Int_t iCl = 0; iCl < nClusters; iCl++) - if (((AliMUONTrackParam*) muonTrack.GetTrackParamAtCluster()->UncheckedAt(iCl))->GetClusterPtr()->GetMCLabel() == label) - nMatchClusters++; - fractionOfMatchCluster = ((Float_t)nMatchClusters) / ((Float_t)nClusters); - } - - } else { // by comparing cluster/TrackRef positions - - // look for the corresponding simulated track if any - TIter next(trackRefStore.CreateIterator()); - AliMUONTrack* trackRef; - while ( ( trackRef = static_cast(next()) ) ) { - - // check compatibility - Float_t f = 0.; - if (TrackMatched(muonTrack, *trackRef, f, sigmaCut)) { - matchedTrackRef = trackRef; - fractionOfMatchCluster = f; - break; - } - - } - - } - - return matchedTrackRef; - -} - -//----------------------------------------------------------------------- -Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, AliMUONRecoParam &recoParam, - Bool_t useLabel, Double_t sigmaCut, TH1F &hFractionOfConnectedClusters) +Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStore &trackRefStore, + Bool_t useLabel, TH1F &hFractionOfConnectedClusters) { /// loop over reconstructible TrackRef not associated with reconstructed track: /// for each of them, find and remove the most connected the fake track, if any, @@ -456,7 +366,7 @@ Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStor while ( ( trackRef = static_cast(next()) ) ) { // skip not reconstructible trackRefs - if (!IsRecontructible(*trackRef,recoParam)) continue; + if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue; Int_t label = trackRef->GetUniqueID(); @@ -475,7 +385,7 @@ Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStor nConnectedClusters++; } else { // by comparing cluster/TrackRef positions Bool_t compTrack[10]; - nConnectedClusters = fakeTrack->CompatibleTrack(trackRef, sigmaCut, compTrack); + nConnectedClusters = fakeTrack->FindCompatibleClusters(*trackRef, sigmaCut, compTrack); } // skip non-connected fake tracks diff --git a/MUON/MUONRecoCheck.C b/MUON/MUONRecoCheck.C index cc62800fb17..c15e3302fba 100644 --- a/MUON/MUONRecoCheck.C +++ b/MUON/MUONRecoCheck.C @@ -26,8 +26,12 @@ /// \author Jean-Pierre Cussonneau, Subatech // ROOT includes +#include "TMath.h" #include "TClonesArray.h" #include "TH1.h" +#include "TH2.h" +#include "TGraphErrors.h" +#include "TF1.h" #include "TFile.h" // STEER includes @@ -39,54 +43,76 @@ #include "AliMUONTrack.h" #include "AliMUONRecoCheck.h" #include "AliMUONTrackParam.h" -#include "AliMUONTrackExtrap.h" #include "AliMUONRecoParam.h" #include "AliMUONVTrackStore.h" - -Int_t TrackCheck( Bool_t *compTrack); +#include "AliMUONVCluster.h" void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esdFileName="AliESDs.root", - const char* ocdbPath = "local://$ALICE_ROOT/OCDB") + const char* ocdbPath = "local://$ALICE_ROOT/OCDB") { - Bool_t compTrack[10]; - Bool_t compTrackOK[10]; - Int_t nClusterOk = 0; - Int_t testTrack = 0; - Int_t iTrack = 0; - AliMUONTrack* trackOK(0x0); - Int_t trackID = 0; - Double_t maxChi2 = 999.; - AliMUONTrackParam *trackParam; - Double_t x1,y1,z1,pX1,pY1,pZ1,p1; - Double_t x2,y2,z2,pX2,pY2,pZ2,p2; - // File for histograms and histogram booking TFile *histoFile = new TFile("MUONRecoCheck.root", "RECREATE"); TH1F *hReconstructible = new TH1F("hReconstructible"," Nb of reconstructible tracks ",15,-0.5,14.5); TH1F *hReco = new TH1F("hReco"," Nb of reconstructed tracks / evt",15,-0.5,14.5); TH1F *hNClusterComp = new TH1F("hNClusterComp"," Nb of compatible clusters / track ",15,-0.5,14.5); - TH1F *hTestTrack = new TH1F("hTestTrack"," Reconstruction requirement / track",15,-0.5,14.5); TH1F *hTrackRefID = new TH1F("hTrackRefID"," track reference ID ",100,-0.5,99.5); + // momentum resolution TH1F *hResMomVertex = new TH1F("hResMomVertex"," delta P vertex (GeV/c)",100,-10.,10); TH1F *hResMomFirstCluster = new TH1F("hResMomFirstCluster"," delta P first cluster (GeV/c)",100,-10.,10); + TH2D *hResMomVertexVsMom = new TH2D("hResMomVertexVsMom","#Delta_{p} at vertex versus p (GeV/c)",30,0.,300.,800,-20.,20.); + TH2D *hResMomFirstClusterVsMom = new TH2D("hResMomFirstClusterVsMom","#Delta_{p} at first cluster versus p (GeV/c)",30,0.,300.,800,-20.,20.); + TGraphErrors* gResMomVertexVsMom = new TGraphErrors(30); + gResMomVertexVsMom->SetName("gResMomVertexVsMom"); + gResMomVertexVsMom->SetTitle("#Delta_{p}/p at vertex versus p;p (GeV/c);#sigma_{p}/p (%)"); + TGraphErrors* gResMomFirstClusterVsMom = new TGraphErrors(30); + gResMomFirstClusterVsMom->SetName("gResMomFirstClusterVsMom"); + gResMomFirstClusterVsMom->SetTitle("#Delta_{p}/p at first cluster versus p;p (GeV/c);#sigma_{p}/p (%)"); + TF1* f = new TF1("f","gaus"); + + // residuals at clusters + histoFile->mkdir("residuals","residuals"); + histoFile->cd("residuals"); + TH1F* hResidualXInCh[AliMUONConstants::NTrackingCh()]; + TH1F* hResidualYInCh[AliMUONConstants::NTrackingCh()]; + for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { + hResidualXInCh[i] = new TH1F(Form("hResidualXInCh%d",i+1), Form("cluster-track residual-X distribution in chamber %d;#Delta_{X} (cm)",i+1), 1000, -1., 1.); + hResidualYInCh[i] = new TH1F(Form("hResidualYInCh%d",i+1), Form("cluster-track residual-Y distribution in chamber %d;#Delta_{Y} (cm)",i+1), 1000, -0.5, 0.5); + } + TGraphErrors* gResidualXPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh()); + gResidualXPerChMean->SetName("gResidualXPerChMean"); + gResidualXPerChMean->SetTitle("cluster-track residual-X per Ch: mean;chamber ID;<#Delta_{Y}> (cm)"); + gResidualXPerChMean->SetMarkerStyle(kFullDotLarge); + TGraphErrors* gResidualYPerChMean = new TGraphErrors(AliMUONConstants::NTrackingCh()); + gResidualYPerChMean->SetName("gResidualYPerChMean"); + gResidualYPerChMean->SetTitle("cluster-track residual-Y per Ch: mean;chamber ID;<#Delta_{Y}> (cm)"); + gResidualYPerChMean->SetMarkerStyle(kFullDotLarge); + TGraphErrors* gResidualXPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh()); + gResidualXPerChSigma->SetName("gResidualXPerChSigma"); + gResidualXPerChSigma->SetTitle("cluster-track residual-X per Ch: sigma;chamber ID;#sigma_{X} (cm)"); + gResidualXPerChSigma->SetMarkerStyle(kFullDotLarge); + TGraphErrors* gResidualYPerChSigma = new TGraphErrors(AliMUONConstants::NTrackingCh()); + gResidualYPerChSigma->SetName("gResidualYPerChSigma"); + gResidualYPerChSigma->SetTitle("cluster-track residual-Y per Ch: sigma;chamber ID;#sigma_{X} (cm)"); + gResidualYPerChSigma->SetMarkerStyle(kFullDotLarge); AliMUONRecoCheck rc(esdFileName, pathSim); // load necessary data from OCDB AliCDBManager::Instance()->SetDefaultStorage(ocdbPath); AliCDBManager::Instance()->SetRun(rc.GetRunNumber()); - AliMUONCDB::LoadField(); AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam(); if (!recoParam) return; - // set the magnetic field for track extrapolations - AliMUONTrackExtrap::SetField(); - - // get sigma cut from recoParam to associate clusters with TrackRefs + // get sigma cut from recoParam to associate clusters with TrackRefs in case the label are not used Double_t sigmaCut = (recoParam->ImproveTracks()) ? recoParam->GetSigmaCutForImprovement() : recoParam->GetSigmaCutForTracking(); + // compute the mask of requested stations from recoParam + UInt_t requestedStationMask = 0; + for (Int_t i = 0; i < 5; i++) if (recoParam->RequestStation(i)) requestedStationMask |= ( 1 << i ); + // get from recoParam whether a track need 2 chambers hit in the same station (4 or 5) or not to be reconstructible + Bool_t request2ChInSameSt45 = !recoParam->MakeMoreTrackCandidates(); Int_t nevents = rc.NumberOfEvents(); @@ -96,13 +122,16 @@ void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esd Int_t nReconstructibleTracks = 0; Int_t nReconstructedTracks = 0; Int_t nReconstructibleTracksCheck = 0; + AliMUONTrackParam *trackParam; + Double_t x1,y1,z1,pX1,pY1,pZ1,p1,pT1; + Double_t x2,y2,z2,pX2,pY2,pZ2,p2,pT2; for (ievent=0; ieventFill(trackRefStore->GetSize()); hReco->Fill(trackStore->GetSize()); @@ -110,68 +139,36 @@ void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esd nReconstructibleTracks += trackRefStore->GetSize(); nReconstructedTracks += trackStore->GetSize(); + // loop over trackRef TIter next(trackRefStore->CreateIterator()); AliMUONTrack* trackRef; - while ( ( trackRef = static_cast(next()) ) ) { - maxChi2 = 999.; - testTrack = 0; - trackOK = 0x0; - for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) - { - compTrackOK[ch] = kFALSE; - } + hTrackRefID->Fill(trackRef->GetMCLabel()); + + AliMUONTrack* trackMatched = 0x0; + Int_t nMatchClusters = 0; + + // loop over trackReco and look for compatible track TIter next2(trackStore->CreateIterator()); AliMUONTrack* trackReco; - while ( ( trackReco = static_cast(next2()) ) ) { - // check if trackRef is compatible with trackReco - if (trackReco->GetNClusters() > 1) { - - // check cluster by cluster if trackReco contain info at each cluster - trackRef->CompatibleTrack(trackReco,sigmaCut,compTrack); - - iTrack = TrackCheck(compTrack); - - if (iTrack > testTrack) - { - nClusterOk = 0; - for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) - { - if (compTrack[ch]) nClusterOk++; - compTrackOK[ch] = compTrack[ch]; - } - testTrack = iTrack; - trackOK = trackReco; - } - - } else { - - // check only parameters at the z position of the first trackRef - AliMUONTrackParam *refParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First(); - AliMUONTrackParam recoParam(*((AliMUONTrackParam*) trackReco->GetTrackParamAtCluster()->First())); - Double_t chi2; - if (AliMUONTrackExtrap::ExtrapToZCov(&recoParam, refParam->GetZ()) && - refParam->CompatibleTrackParam(recoParam, sigmaCut, chi2) && chi2 < maxChi2) - { - maxChi2 = chi2; - trackOK = trackReco; - } - + + // check if trackReco is compatible with trackRef + Int_t n = 0; + if (trackReco->Match(*trackRef, sigmaCut, nMatchClusters)) { + trackMatched = trackReco; + nMatchClusters = n; + break; } - + } - hTestTrack->Fill(testTrack); - trackID = trackRef->GetMCLabel(); - hTrackRefID->Fill(trackID); - - if (testTrack == 4 || maxChi2 < 5.*sigmaCut*sigmaCut) { // tracking requirements verified, track is found + if (trackMatched) { // tracking requirements verified, track is found nReconstructibleTracksCheck++; - hNClusterComp->Fill(nClusterOk); + hNClusterComp->Fill(nMatchClusters); trackParam = trackRef->GetTrackParamAtVertex(); x1 = trackParam->GetNonBendingCoor(); y1 = trackParam->GetBendingCoor(); @@ -180,10 +177,10 @@ void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esd pY1 = trackParam->Py(); pZ1 = trackParam->Pz(); p1 = trackParam->P(); + pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1); // printf(" Ref. track at vertex: x,y,z: %f %f %f px,py,pz,p: %f %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1,p1); - trackReco = trackOK; - trackParam = new AliMUONTrackParam(*((AliMUONTrackParam*)(trackReco->GetTrackParamAtVertex()))); + trackParam = trackMatched->GetTrackParamAtVertex(); x2 = trackParam->GetNonBendingCoor(); y2 = trackParam->GetBendingCoor(); z2 = trackParam->GetZ(); @@ -191,11 +188,12 @@ void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esd pY2 = trackParam->Py(); pZ2 = trackParam->Pz(); p2 = trackParam->P(); - delete trackParam; + pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2); // printf(" Reconst. track at vertex: x,y,z: %f %f %f px,py,pz: %f %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2,p2); hResMomVertex->Fill(p2-p1); - + hResMomVertexVsMom->Fill(p1,p2-p1); + trackParam = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First(); x1 = trackParam->GetNonBendingCoor(); y1 = trackParam->GetBendingCoor(); @@ -204,8 +202,10 @@ void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esd pY1 = trackParam->Py(); pZ1 = trackParam->Pz(); p1 = trackParam->P(); + pT1 = TMath::Sqrt(pX1*pX1 + pY1*pY1); + // printf(" Ref. track at 1st cluster: x,y,z: %f %f %f px,py,pz: %f %f %f \n",x1,y1,z1,pX1,pY1,pZ1); - trackParam = (AliMUONTrackParam*) trackOK->GetTrackParamAtCluster()->First(); + trackParam = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First(); x2 = trackParam->GetNonBendingCoor(); y2 = trackParam->GetBendingCoor(); z2 = trackParam->GetZ(); @@ -213,42 +213,83 @@ void MUONRecoCheck (Int_t nEvent = -1, char * pathSim="./generated/", char * esd pY2 = trackParam->Py(); pZ2 = trackParam->Pz(); p2 = trackParam->P(); + pT2 = TMath::Sqrt(pX2*pX2 + pY2*pY2); // printf(" Reconst. track at 1st cluster: x,y,z: %f %f %f px,py,pz: %f %f %f \n",x2,y2,z2,pX2,pY2,pZ2); hResMomFirstCluster->Fill(p2-p1); - + hResMomFirstClusterVsMom->Fill(p1,p2-p1); + + // Fill residuals + // Loop over clusters of first track + AliMUONTrackParam* trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->First(); + while (trackParamAtCluster1) { + AliMUONVCluster* cluster1 = trackParamAtCluster1->GetClusterPtr(); + AliMUONTrackParam* trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->First(); + while (trackParamAtCluster2) { + AliMUONVCluster* cluster2 = trackParamAtCluster2->GetClusterPtr(); + if (cluster1->GetDetElemId() == cluster2->GetDetElemId()) { + hResidualXInCh[cluster1->GetChamberId()]->Fill(cluster1->GetX() - cluster2->GetX()); + hResidualYInCh[cluster1->GetChamberId()]->Fill(cluster1->GetY() - cluster2->GetY()); + break; + } + trackParamAtCluster2 = (AliMUONTrackParam*) trackRef->GetTrackParamAtCluster()->After(trackParamAtCluster2); + } + trackParamAtCluster1 = (AliMUONTrackParam*) trackMatched->GetTrackParamAtCluster()->After(trackParamAtCluster1); + } + } + } // end loop track ref. } // end loop on event + // compute momentum resolution versus p + for (Int_t i = 1; i <= hResMomVertexVsMom->GetNbinsX(); i++) { + TH1D *tmp = hResMomVertexVsMom->ProjectionY("tmp",i,i,"e"); + Double_t p = hResMomVertexVsMom->GetBinCenter(i); + gResMomVertexVsMom->SetPoint(i-1,p,100.*tmp->GetRMS()/p); + gResMomVertexVsMom->SetPointError(i-1,hResMomVertexVsMom->GetBinWidth(i)/2.,100.*tmp->GetRMSError()/p); + delete tmp; + } + for (Int_t i = 1; i <= hResMomFirstClusterVsMom->GetNbinsX(); i++) { + TH1D *tmp = hResMomFirstClusterVsMom->ProjectionY("tmp",i,i,"e"); + Double_t p = hResMomFirstClusterVsMom->GetBinCenter(i); + f->SetParameters(1.,0.,1.); + f->SetParLimits(0,0.,1.e3); + tmp->Fit("f","WWN"); + gResMomFirstClusterVsMom->SetPoint(i-1,p,100.*f->GetParameter(2)/p); + gResMomFirstClusterVsMom->SetPointError(i-1,hResMomFirstClusterVsMom->GetBinWidth(i)/2.,100.*f->GetParError(2)/p); + delete tmp; + } + + // compute residual mean and dispersion + for (Int_t i = 0; i < AliMUONConstants::NTrackingCh(); i++) { + hResidualXInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualXInCh[i]->GetRMS(), 3.*hResidualXInCh[i]->GetRMS()); + gResidualXPerChMean->SetPoint(i, i+1, hResidualXInCh[i]->GetMean()); + gResidualXPerChMean->SetPointError(i, 0., hResidualXInCh[i]->GetMeanError()); + gResidualXPerChSigma->SetPoint(i, i+1, hResidualXInCh[i]->GetRMS()); + gResidualXPerChSigma->SetPointError(i, 0., hResidualXInCh[i]->GetRMSError()); + hResidualXInCh[i]->GetXaxis()->SetRange(0,0); + hResidualYInCh[i]->GetXaxis()->SetRangeUser(-3.*hResidualYInCh[i]->GetRMS(), 3.*hResidualYInCh[i]->GetRMS()); + gResidualYPerChMean->SetPoint(i, i+1, hResidualYInCh[i]->GetMean()); + gResidualYPerChMean->SetPointError(i, 0., hResidualYInCh[i]->GetMeanError()); + gResidualYPerChSigma->SetPoint(i, i+1, hResidualYInCh[i]->GetRMS()); + gResidualYPerChSigma->SetPointError(i, 0., hResidualYInCh[i]->GetRMSError()); + hResidualYInCh[i]->GetXaxis()->SetRange(0,0); + } + printf(" nb of reconstructible tracks: %d \n", nReconstructibleTracks); printf(" nb of reconstructed tracks: %d \n", nReconstructedTracks); printf(" nb of reconstructible tracks which are reconstructed: %d \n", nReconstructibleTracksCheck); histoFile->Write(); + histoFile->cd(); + gResMomVertexVsMom->Write(); + gResMomFirstClusterVsMom->Write(); + gResidualXPerChMean->Write(); + gResidualXPerChSigma->Write(); + gResidualYPerChMean->Write(); + gResidualYPerChSigma->Write(); histoFile->Close(); } - -Int_t TrackCheck( Bool_t *compTrack) -{ - // Apply reconstruction requirements - // Return number of validated conditions - // If all the tests are verified then TrackCheck = 4 (good track) - Int_t iTrack = 0; - Int_t nCompClustersInLastStations = 0; - - // apply reconstruction requirements - if (compTrack[0] || compTrack[1]) iTrack++; // at least one compatible cluster in st. 0 - if (compTrack[2] || compTrack[3]) iTrack++; // at least one compatible cluster in st. 1 - if (compTrack[4] || compTrack[5]) iTrack++; // at least one compatible cluster in st. 2 - for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) { - if (compTrack[ch]) nCompClustersInLastStations++; - } - if (nCompClustersInLastStations > 2) iTrack++; // at least 3 compatible clusters in st. 3 & 4 - - return iTrack; - -} - diff --git a/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx b/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx index 26b0d8d6c49..1d291c461ea 100644 --- a/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx +++ b/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx @@ -107,23 +107,21 @@ void AliAnalysisTaskESDMCLabelAddition::AddMCLabel() for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) { esdTrack = esd->GetMuonTrack(nMuTrack); - esdTrack->SetLabel(-1); // skip ghosts if (!esdTrack->ContainTrackerData()) continue; + // convert ESD track to MUON track (without recomputing track parameters at each clusters) + AliMUONTrack muonTrack; + AliMUONESDInterface::ESDToMUON(*esdTrack, muonTrack, kFALSE); + // try to match the reconstructed track with a simulated one - AliMUONTrack* matchedTrackRef = MatchWithTrackRef(*esdTrack, *trackRefStore); + Int_t nMatchClusters = 0; + AliMUONTrack* matchedTrackRef = rc.FindCompatibleTrack(muonTrack, *trackRefStore, nMatchClusters, kFALSE, fgkSigmaCut); - if (matchedTrackRef) { - - // set the MC label - esdTrack->SetLabel(matchedTrackRef->GetUniqueID()); - - // remove already matched trackRefs - trackRefStore->Remove(*matchedTrackRef); - - } + // set the MC label + if (matchedTrackRef) esdTrack->SetLabel(matchedTrackRef->GetUniqueID()); + else esdTrack->SetLabel(-1); } @@ -138,97 +136,3 @@ void AliAnalysisTaskESDMCLabelAddition::Terminate(Option_t */*option*/) AliDebug(2, "Terminate()"); } - -//---------------------------------------------------------------------- -AliMUONTrack* AliAnalysisTaskESDMCLabelAddition::ESDToMUON(AliESDMuonTrack &esdTrack) -{ - /// Convert an ESD track into a MUON track with dummy parameters (just to fill the clusters). - /// It is the responsability of the user to delete the track afterward - - AliMUONTrack *track = new AliMUONTrack(); - - // ckeck whether the ESD track contains clusters - if(!esdTrack.ClustersStored()) return track; - - // track parameters at first cluster - AliMUONTrackParam param; - AliMUONESDInterface::GetParamAtFirstCluster(esdTrack, param); - AliMUONESDInterface::GetParamCov(esdTrack, param); - - // create empty cluster - AliMUONVClusterStore* cStore = AliMUONESDInterface::NewClusterStore(); - AliMUONVCluster* cluster = cStore->CreateCluster(0,0,0); - - // loop over ESD clusters - AliESDMuonCluster *esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().First(); - while (esdCluster) { - - // copy cluster information - AliMUONESDInterface::ESDToMUON(*esdCluster, *cluster); - - // only set the Z parameter to avoid error in the AddTrackParamAtCluster(...) method - param.SetZ(cluster->GetZ()); - - // add common track parameters at current cluster - track->AddTrackParamAtCluster(param, *cluster, kTRUE); - - esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().After(esdCluster); - } - - // clean memory - delete cluster; - delete cStore; - - return track; - -} - - -//---------------------------------------------------------------------- -AliMUONTrack* AliAnalysisTaskESDMCLabelAddition::MatchWithTrackRef(AliESDMuonTrack &esdTrack, - AliMUONVTrackStore &trackRefStore) -{ - /// Return the trackRef matched with the reconstructed track and the fraction of matched clusters - - AliMUONTrack *matchedTrackRef = 0x0; - - // convert ESD track to MUON track - AliMUONTrack *track = ESDToMUON(esdTrack); - - // look for the corresponding simulated track if any - TIter next(trackRefStore.CreateIterator()); - AliMUONTrack* trackRef; - while ( ( trackRef = static_cast(next()) ) ) { - - // check compatibility - if (TrackMatched(*track, *trackRef)) { - matchedTrackRef = trackRef; - break; - } - - } - - // clean memory - delete track; - - return matchedTrackRef; - -} - - -//---------------------------------------------------------------------- -Bool_t AliAnalysisTaskESDMCLabelAddition::TrackMatched(AliMUONTrack &track, AliMUONTrack &trackRef) -{ - /// Try to match 2 tracks - - Bool_t compTrack[10]; - Int_t nMatchClusters = track.CompatibleTrack(&trackRef, fgkSigmaCut, compTrack); - Double_t fractionOfMatchCluster = ((Double_t)nMatchClusters) / ((Double_t)track.GetNClusters()); - - if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2 - (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5 - fractionOfMatchCluster > 0.5) return kTRUE; // more than 50% of clusters matched - else return kFALSE; - -} - diff --git a/PWG3/muondep/AliAnalysisTaskRecoCheck.cxx b/PWG3/muondep/AliAnalysisTaskRecoCheck.cxx index 7555104c52c..031f21e7b74 100644 --- a/PWG3/muondep/AliAnalysisTaskRecoCheck.cxx +++ b/PWG3/muondep/AliAnalysisTaskRecoCheck.cxx @@ -4,8 +4,6 @@ #include "TIterator.h" #include "AliStack.h" -#include "AliMagF.h" -#include "AliTracker.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" @@ -23,8 +21,8 @@ #include "AliMUONPairLight.h" #include "AliMUONVTrackStore.h" #include "AliMUONTrack.h" +#include "AliMUONESDInterface.h" #include "AliMUONRecoCheck.h" -#include "AliMUONTrackExtrap.h" // analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck) // Authors: Bogdan Vulpescu @@ -67,9 +65,6 @@ void AliAnalysisTaskRecoCheck::ConnectInputData(Option_t *) fESDEvent = esdH->GetEvent(); } - // calculate the filed map in the L3 magnet using the current value - AliMUONTrackExtrap::SetField(); - } //________________________________________________________________________ @@ -111,8 +106,6 @@ void AliAnalysisTaskRecoCheck::Exec(Option_t *) return; } - Int_t eventNumber = fESDEvent->GetEventNumberInFile(); - fArray1Mu->Clear(); fArray2Mu->Clear(); @@ -120,42 +113,43 @@ void AliAnalysisTaskRecoCheck::Exec(Option_t *) AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler); - AliMUONVTrackStore* trackRefs = rc->TrackRefs(eventNumber); - AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(eventNumber); + AliMUONVTrackStore* trackRefs = rc->TrackRefs(-1); AliStack* pstack = mcEvent->Stack(); - TIter next(recoTracks->CreateIterator()); - AliMUONTrack* trackReco = NULL; - - Int_t nTrackReco = recoTracks->GetSize(); - Int_t nTracksESD = fESDEvent->GetNumberOfMuonTracks(); - //printf("Tracks in recoTracks (%d) and in ESD (%d).\n", nTrackReco, nTracksESD); - - if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD); - + // loop over ESD tracks Int_t nreftracks = 0; - Int_t itrRec = 0; - while ( (trackReco = static_cast(next())) != NULL ) { - // assign parameters concerning the reconstructed tracks - AliMUONTrackLight muLight; + for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { + + AliESDMuonTrack* esdTrack = fESDEvent->GetMuonTrack(iTrack); - muLight.FillFromESD(fESDEvent->GetMuonTrack(itrRec)); + // skip ghosts + if (!esdTrack->ContainTrackerData()) continue; - // find the reference track and store further information - TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack); - if (part) { + // convert ESD track to MUON track (without recomputing track parameters at each clusters) + AliMUONTrack muonTrack; + AliMUONESDInterface::ESDToMUON(*esdTrack, muonTrack, kFALSE); + + // try to match the reconstructed track with a simulated one + Int_t nMatchClusters = 0; + AliMUONTrack* matchedTrackRef = rc->FindCompatibleTrack(muonTrack, *trackRefs, nMatchClusters, kFALSE, 10.); + + if (matchedTrackRef) { + + //store new referenced track in the muonArray with parameters from the reconstructed tracks + AliMUONTrackLight* muLight = new ((*fArray1Mu)[nreftracks++]) AliMUONTrackLight(esdTrack); + + // store further information related to the simulated track + muLight->SetTrackPythiaLine(matchedTrackRef->GetUniqueID()); + TParticle *part = pstack->Particle(matchedTrackRef->GetUniqueID()); + muLight->SetTrackPDGCode(part->GetPdgCode()); v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy()); - muLight.SetPGen(v); - muLight.FillMuonHistory(pstack, part); - //muLight.PrintInfo("A"); - //store the referenced track in the muonArray: - TClonesArray &muons = *fArray1Mu; - new (muons[nreftracks++]) AliMUONTrackLight(muLight); - } + muLight->SetPGen(v); + muLight->FillMuonHistory(pstack, part); + + } - itrRec++; - } // end reco track - + } // end esd track + // now loop over muon pairs to build dimuons Int_t nmuons = fArray1Mu->GetEntriesFast(); Int_t ndimuons = 0; -- 2.43.0