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.)
}
//_____________________________________________________________________________
-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;
}
if (fRecoTrackRefStore != 0x0) return fRecoTrackRefStore;
else {
if (TrackRefs(event) == 0x0) return 0x0;
- MakeReconstructibleTracks();
+ MakeReconstructibleTracks(requestedStationMask, request2ChInSameSt45);
return fRecoTrackRefStore;
}
}
}
//_____________________________________________________________________________
-void AliMUONRecoCheck::MakeReconstructibleTracks()
+void AliMUONRecoCheck::MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
{
/// Isolate the reconstructible tracks
if (!(fRecoTrackRefStore = AliMUONESDInterface::NewTrackStore())) return;
// 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<AliMUONTrack*>(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<AliMUONTrack*>(next()) ) ) {
+
+ // check compatibility
+ Int_t n = 0;
+ if (track.Match(*track2, sigmaCut, n)) {
+ matchedTrack = track2;
+ nMatchClusters = n;
+ break;
+ }
+
+ }
- delete [] chamberInTrack;
}
-
+
+ return matchedTrack;
+
}
class AliESDEvent;
class AliMCEventHandler;
class AliMUONVTrackStore;
+class AliMUONTrack;
class AliMUONRecoCheck : public TObject
{
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();
/// 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);
void CleanMuonTrackRef(const AliMUONVTrackStore *tmpTrackRefStore);
- void MakeReconstructibleTracks();
+ void MakeReconstructibleTracks(UInt_t requestedStationMask, Bool_t request2ChInSameSt45 = kTRUE);
private:
AliMCEventHandler* fMCEventHandler; ///< to access MC truth information
}
//__________________________________________________________________________
-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
// 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);
+
}
//__________________________________________________________________________
}
//__________________________________________________________________________
-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;
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();
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;
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)
{
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);
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;}
#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"
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);
}
//============================================
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);
}
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<AliMUONTrack*>(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
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);
/// 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,
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;
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) {
}
-//-----------------------------------------------------------------------
-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<AliMUONTrack*>(next()) ) ) {
-
- // check compatibility
- Float_t f = 0.;
- if (TrackMatched(muonTrack, *trackRef, f, sigmaCut)) {
- matchedTrackRef = trackRef;
- fractionOfMatchCluster = f;
- break;
- }
-
- }
-
- }
-
- return matchedTrackRef;
-
-}
-
#include <TLorentzVector.h>
#include <TParticle.h>
#include <TSystem.h>
-#include <TGeoManager.h>
-#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"
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.
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");
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
AliMUONVTrackStore* trackRefs = rc->TrackRefs(ievent);
AliStack* pstack = (const_cast<AliMCEventHandler*>(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<AliMUONTrack*>(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();
/// 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,
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;
TIter next(trackRefStore->CreateIterator());
AliMUONTrack* trackRef;
while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
- if (IsRecontructible(*trackRef,*recoParam)) nReconstructibleTracks++;
+ if (trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) nReconstructibleTracks++;
}
// loop over ESD tracks
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);
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;
}
//-----------------------------------------------------------------------
-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<AliMUONTrackParam*>(track.GetTrackParamAtCluster()->First());
- while (trackParam) {
-
- currentCh = trackParam->GetClusterPtr()->GetChamberId();
-
- clusterInSt[currentCh/2] = kTRUE;
-
- if (currentCh > 5 && currentCh != previousCh) {
- nChHitInSt45++;
- previousCh = currentCh;
- }
-
- trackParam = static_cast<AliMUONTrackParam*>(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<AliMUONTrack*>(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,
while ( ( trackRef = static_cast<AliMUONTrack*>(next()) ) ) {
// skip not reconstructible trackRefs
- if (!IsRecontructible(*trackRef,recoParam)) continue;
+ if (!trackRef->IsValid(requestedStationMask, request2ChInSameSt45)) continue;
Int_t label = trackRef->GetUniqueID();
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
/// \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
#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();
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; ievent<nEvent; ievent++)
{
if (!(ievent%10)) printf(" **** event # %d \n",ievent);
AliMUONVTrackStore* trackStore = rc.ReconstructedTracks(ievent, kFALSE);
- AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent);
+ AliMUONVTrackStore* trackRefStore = rc.ReconstructibleTracks(ievent, requestedStationMask, request2ChInSameSt45);
hReconstructible->Fill(trackRefStore->GetSize());
hReco->Fill(trackStore->GetSize());
nReconstructibleTracks += trackRefStore->GetSize();
nReconstructedTracks += trackStore->GetSize();
+ // loop over trackRef
TIter next(trackRefStore->CreateIterator());
AliMUONTrack* trackRef;
-
while ( ( trackRef = static_cast<AliMUONTrack*>(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<AliMUONTrack*>(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();
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();
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();
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();
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;
-
-}
-
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);
}
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<AliMUONTrack*>(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;
-
-}
-
#include "TIterator.h"
#include "AliStack.h"
-#include "AliMagF.h"
-#include "AliTracker.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#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
fESDEvent = esdH->GetEvent();
}
- // calculate the filed map in the L3 magnet using the current value
- AliMUONTrackExtrap::SetField();
-
}
//________________________________________________________________________
return;
}
- Int_t eventNumber = fESDEvent->GetEventNumberInFile();
-
fArray1Mu->Clear();
fArray2Mu->Clear();
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<AliMUONTrack*>(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;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff