Fixes in reconstruction:
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Nov 2009 12:51:07 +0000 (12:51 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Nov 2009 12:51:07 +0000 (12:51 +0000)
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.)

12 files changed:
MUON/AliMUONRecoCheck.cxx
MUON/AliMUONRecoCheck.h
MUON/AliMUONTrack.cxx
MUON/AliMUONTrack.h
MUON/AliMUONTrackLight.cxx
MUON/AliMUONTrackLight.h
MUON/DIMUONFakes.C
MUON/DecodeRecoCocktail.C
MUON/MUONFakes.C
MUON/MUONRecoCheck.C
PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx
PWG3/muondep/AliAnalysisTaskRecoCheck.cxx

index 0daff81..aaef0b0 100644 (file)
@@ -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<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;
+  
 }
 
index 0897c1b..1d0aa41 100644 (file)
@@ -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
index 8bcc1c9..d402a98 100644 (file)
@@ -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;
@@ -1117,6 +1132,23 @@ Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool
 }
 
 //__________________________________________________________________________
+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)
 {
   /// set track parameters at vertex
index 943f52a..e6f5a5d 100644 (file)
@@ -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;}
index 2ac9448..ff842c1 100644 (file)
 #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); 
 }
@@ -192,81 +197,6 @@ void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
 }
 
 //============================================
-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
   Int_t countP = -1;
index fe7a49b..73c982c 100644 (file)
@@ -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);
index 95b7ace..1c15a65 100644 (file)
@@ -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<AliMUONTrack*>(next()) ) ) {
-      
-      // check compatibility
-      Float_t f = 0.;
-      if (TrackMatched(muonTrack, *trackRef, f, sigmaCut)) {
-       matchedTrackRef = trackRef;
-       fractionOfMatchCluster = f;
-       break;
-      }
-      
-    }
-    
-  }
-  
-  return matchedTrackRef;
-  
-}
-
index 6b2b7df..118f715 100644 (file)
 #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"
@@ -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<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(); 
index 60ff78f..271938d 100644 (file)
 /// 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<AliMUONTrack*>(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<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,
@@ -456,7 +366,7 @@ Int_t RemoveConnectedFakes(AliMUONVTrackStore &fakeTrackStore, AliMUONVTrackStor
   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();
     
@@ -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
index cc62800..c15e330 100644 (file)
 /// \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();
   
@@ -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; 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());
@@ -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<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();
@@ -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;
-  
-}
-
index 26b0d8d..1d291c4 100644 (file)
@@ -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<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;
-  
-}
-
index 7555104..031f21e 100644 (file)
@@ -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<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;