1) New class "AliMUONRefitter" to:
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 20 Mar 2008 13:34:13 +0000 (13:34 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 20 Mar 2008 13:34:13 +0000 (13:34 +0000)
   - re-clusterize the ESD clusters from the attached ESD pad
   - refit the ESD tracks using the attached ESD clusters or
     the re-clusterized ones
2) The macro MUONRefitter.C provides a complete example using
   this new class (including modification of the digit charge
   before refitting, comparison of results with the original data
   and storage of these results in a new AliESDs.root file).
(Philippe P.)

18 files changed:
MUON/AliMUONClusterFinderMLEM.cxx
MUON/AliMUONClusterFinderPeakCOG.cxx
MUON/AliMUONClusterFinderPeakFit.cxx
MUON/AliMUONRefitter.cxx [new file with mode: 0644]
MUON/AliMUONRefitter.h [new file with mode: 0644]
MUON/AliMUONSimpleClusterServer.cxx
MUON/AliMUONSimpleClusterServer.h
MUON/AliMUONTrackReconstructor.cxx
MUON/AliMUONTrackReconstructor.h
MUON/AliMUONTrackReconstructorK.cxx
MUON/AliMUONTrackReconstructorK.h
MUON/AliMUONTracker.cxx
MUON/AliMUONVClusterServer.h
MUON/AliMUONVTrackReconstructor.h
MUON/MUONRefit.C [new file with mode: 0644]
MUON/MUONrecLinkDef.h
MUON/TestRecPoints.C
MUON/libMUONrec.pkg

index 67884aa..c4142b0 100644 (file)
@@ -125,7 +125,8 @@ AliMUONClusterFinderMLEM::Prepare(Int_t detElemId,
   fSplitter->SetDebug(fDebug);
     
   // find out current event number, and reset the cluster number
-  fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
+  AliRunLoader *runLoader = AliRunLoader::GetRunLoader();
+  fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
   fClusterNumber = -1;
   fClusterList.Delete();
   
@@ -198,7 +199,7 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster()
   //  AliCodeTimerAuto("")     
 
   if (fDebug) {
-    cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber() 
+    cout << " *** Event # " << fEventNumber 
         << " det. elem.: " << fDetElemId << endl;
     for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
       AliMUONPad* pad = fPreCluster->Pad(j);
index 94f245b..0790944 100644 (file)
@@ -112,7 +112,8 @@ AliMUONClusterFinderPeakCOG::Prepare(Int_t detElemId, TClonesArray* pads[2],
   fDetElemId = detElemId;
   
   // find out current event number, and reset the cluster number
-  fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
+  AliRunLoader *runLoader = AliRunLoader::GetRunLoader();
+  fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
   fClusterNumber = -1;
   fClusterList.Delete();
   
@@ -185,7 +186,7 @@ AliMUONClusterFinderPeakCOG::WorkOnPreCluster()
   //  AliCodeTimerAuto("")     
 
   if (fDebug) {
-    cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber() 
+    cout << " *** Event # " << fEventNumber 
         << " det. elem.: " << fDetElemId << endl;
     for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
       AliMUONPad* pad = fPreCluster->Pad(j);
index c44b02b..7a370a1 100644 (file)
@@ -185,7 +185,8 @@ AliMUONClusterFinderPeakFit::Prepare(Int_t detElemId, TClonesArray* pads[2],
   fDetElemId = detElemId;
   
   // find out current event number, and reset the cluster number
-  fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
+  AliRunLoader *runLoader = AliRunLoader::GetRunLoader();
+  fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
   fClusterNumber = -1;
   fClusterList.Delete();
   
@@ -278,7 +279,7 @@ AliMUONClusterFinderPeakFit::WorkOnPreCluster()
   //  AliCodeTimerAuto("")     
 
   if (fDebug) {
-    cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber() 
+    cout << " *** Event # " << fEventNumber
         << " det. elem.: " << fDetElemId << endl;
     for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
       AliMUONPad* pad = fPreCluster->Pad(j);
diff --git a/MUON/AliMUONRefitter.cxx b/MUON/AliMUONRefitter.cxx
new file mode 100644 (file)
index 0000000..8764784
--- /dev/null
@@ -0,0 +1,445 @@
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  *
+**************************************************************************/
+
+// $Id$
+
+#include "AliMUONRefitter.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONClusterFinderCOG.h"
+#include "AliMUONClusterFinderMLEM.h"
+#include "AliMUONClusterFinderSimpleFit.h"
+#include "AliMUONPreClusterFinder.h"
+#include "AliMUONPreClusterFinderV2.h"
+#include "AliMUONPreClusterFinderV3.h"
+#include "AliMUONSimpleClusterServer.h"
+#include "AliMUONTrackReconstructor.h"
+#include "AliMUONTrackReconstructorK.h"
+#include "AliMUONRecoParam.h"
+#include "AliMUONESDInterface.h"
+#include "AliMUONVClusterStore.h"
+#include "AliMUONVTrackStore.h"
+#include "AliMUONTrack.h"
+
+#include "AliLog.h"
+
+//-----------------------------------------------------------------------------
+/// \class AliMUONRefitter
+///
+/// create new MUON object from ESD objects given as input (through the ESDInterface):
+///
+/// - re-clusterize the ESD clusters using the attached ESD pads
+///   (several new clusters can be reconstructed per ESD cluster)
+/// - re-fit the ESD tracks using the attached ESD clusters
+/// - reconstruct the ESD tracks from ESD pads (i.e. re-clusterized the attached clusters)
+///
+/// note:
+/// - connexion between an ESD cluster and corresponding MUON clusters from re-clustering
+///   can be made through the detection element ID
+/// - connexion between an ESD track and the corresponding refitted MUON track
+///   can be made through their unique ID
+///
+/// \author Philippe Pillot
+//-----------------------------------------------------------------------------
+
+/// \cond CLASSIMP
+ClassImp(AliMUONRefitter)
+/// \endcond
+
+//_____________________________________________________________________________
+AliMUONRefitter::AliMUONRefitter()
+: TObject(),
+  fGeometryTransformer(0x0),
+  fClusterServer(0x0),
+  fTracker(0x0),
+  fESDInterface(0x0)
+{
+  /// Default constructor
+  CreateGeometryTransformer();
+  CreateClusterServer(*fGeometryTransformer);
+  if (fClusterServer) CreateTrackReconstructor(*fClusterServer);
+  if (!fClusterServer || !fTracker) {
+    AliFatal("refitter initialization failed");
+    exit(-1);
+  }
+}
+
+//_____________________________________________________________________________
+AliMUONRefitter::~AliMUONRefitter()
+{
+  /// Destructor
+  delete fGeometryTransformer;
+  delete fClusterServer;
+  delete fTracker;
+}
+
+//_____________________________________________________________________________
+AliMUONVTrackStore* AliMUONRefitter::ReconstructFromDigits()
+{
+  /// re-reconstruct all tracks and attached clusters from the digits
+  /// it is the responsability of the user to delete the returned store
+  
+  if (!fESDInterface) {
+    AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
+    return 0x0;
+  }
+  
+  // prepare new track(s)
+  AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
+  if (!newTrackStore) return 0x0;
+  
+  // loop over tracks and refit them (create new tracks)
+  Int_t nTracks = fESDInterface->GetNTracks();
+  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
+    AliMUONTrack *track = RetrackFromDigits(iTrack);
+    newTrackStore->Add(track);
+    delete track;
+  }
+  
+  return newTrackStore;
+}
+
+//_____________________________________________________________________________
+AliMUONVTrackStore* AliMUONRefitter::ReconstructFromClusters()
+{
+  /// refit all tracks from the attached clusters
+  /// it is the responsability of the user to delete the returned store
+  
+  if (!fESDInterface) {
+    AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
+    return 0x0;
+  }
+  
+  // prepare new track(s)
+  AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
+  if (!newTrackStore) return 0x0;
+  
+  // loop over tracks and refit them (create new tracks)
+  AliMUONTrack *track;
+  TIter next(fESDInterface->CreateTrackIterator());
+  while ((track = static_cast<AliMUONTrack*>(next()))) {
+    AliMUONTrack* newTrack = newTrackStore->Add(*track);
+    if (!fTracker->RefitTrack(*newTrack)) newTrackStore->Remove(*newTrack);
+  }
+  
+  return newTrackStore;
+}
+
+//_____________________________________________________________________________
+AliMUONTrack* AliMUONRefitter::RetrackFromDigits(Int_t iTrack)
+{
+  /// refit track "iTrack" from the digits (i.e. re-clusterized the attached clusters):
+  /// several new clusters may be reconstructed per initial ESD cluster:
+  /// -> all the combinations of clusters are considered to build the new tracks
+  /// -> return the best track (largest number of clusters or best chi2 in case of equality)
+  /// it is the responsability of the user to delete the returned track
+  
+  if (!fESDInterface) {
+    AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
+    return 0x0;
+  }
+  
+  // get the track to refit
+  AliMUONTrack* track = fESDInterface->GetTrack(iTrack);
+  if (!track) return 0x0;
+  
+  // check if digits exist
+  if (fESDInterface->GetNDigits(iTrack) == 0) {
+    AliError(Form("no digit attached to track #%d",iTrack));
+    return 0x0;
+  }
+  
+  // prepare new track(s)
+  AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
+  if (!newTrackStore) return 0x0;
+  newTrackStore->Add(*track)->Clear("C");
+  
+  // prepare new cluster store
+  AliMUONVClusterStore* newClusterStore = AliMUONESDInterface::NewClusterStore();
+  if (!newClusterStore) {
+    delete newTrackStore;
+    return 0x0;
+  }
+  
+  // loop over clusters, re-clusterize and build new tracks
+  Int_t nClusters = track->GetNClusters();
+  for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
+    
+    // reset the new cluster store
+    newClusterStore->Clear();
+    
+    // get the current cluster
+    AliMUONVCluster* cluster = fESDInterface->GetClusterFast(iTrack,iCluster);
+    
+    // re-clusterize current cluster
+    TIter next(fESDInterface->CreateDigitIterator(iTrack, iCluster));
+    fClusterServer->UseDigits(next);
+    Int_t nNewClusters = fClusterServer->Clusterize(cluster->GetChamberId(),*newClusterStore,AliMpArea());
+    
+    // check that re-clusterizing succeeded
+    if (nNewClusters == 0) {
+      AliWarning(Form("refit gave no cluster (chamber %d)",cluster->GetChamberId()));
+      AliInfo("initial ESD cluster:");
+      cluster->Print("FULL");
+      continue;
+    }
+    
+    // add the new cluster(s) to the tracks
+    AddClusterToTracks(*newClusterStore, *newTrackStore);
+    
+  }
+  
+  // refit the tracks and pick up the best one
+  AliMUONTrack *currentTrack, *bestTrack = 0x0;
+  Double_t currentChi2, bestChi2 = 1.e10;
+  Int_t currentNCluster, bestNClusters = 0;
+  TIter next(newTrackStore->CreateIterator());
+  while ((currentTrack = static_cast<AliMUONTrack*>(next()))) {
+    
+    // set the track parameters at first cluster if any (used as seed in original tracking)
+    AliMUONTrackParam* param = (AliMUONTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
+    if (param) *param = *((AliMUONTrackParam*) track->GetTrackParamAtCluster()->First());
+    
+    // refit the track
+    if (!fTracker->RefitTrack(*currentTrack)) break;
+    
+    // find best track (the one with the higher number of cluster or the best chi2 in case of equality)
+    currentNCluster = currentTrack->GetNClusters();
+    currentChi2 = currentTrack->GetGlobalChi2();
+    if (currentNCluster > bestNClusters || (currentNCluster == bestNClusters && currentChi2 < bestChi2)) {
+      bestTrack = currentTrack;
+      bestNClusters = currentNCluster;
+      bestChi2 = currentChi2;
+    }
+    
+  }
+  
+  // copy best track and free memory
+  AliMUONTrack* newTrack = bestTrack ? new AliMUONTrack(*bestTrack) : 0x0;
+  delete newClusterStore;
+  delete newTrackStore;
+  
+  return newTrack;
+}
+
+//_____________________________________________________________________________
+AliMUONTrack* AliMUONRefitter::RetrackFromClusters(Int_t iTrack)
+{
+  /// refit track "iTrack" form the clusters (i.e. do not re-clusterize)
+  /// it is the responsability of the user to delete the returned track
+  
+  if (!fESDInterface) {
+    AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
+    return 0x0;
+  }
+  
+  // get the track to refit
+  AliMUONTrack* track = fESDInterface->GetTrack(iTrack);
+  if (!track) return 0x0;
+  
+  // refit the track (create a new one)
+  AliMUONTrack* newTrack = new AliMUONTrack(*track);
+  if (!fTracker->RefitTrack(*newTrack)) {
+    delete newTrack;
+    return 0x0;
+  }
+  
+  return newTrack;
+}
+
+//_____________________________________________________________________________
+AliMUONVClusterStore* AliMUONRefitter::ReClusterize(Int_t iTrack, Int_t iCluster)
+{
+  /// re-clusterize cluster numbered "iCluster" in track "iTrack"
+  /// several new clusters may be reconstructed
+  /// it is the responsability of the user to delete the returned store
+  
+  if (!fESDInterface) {
+    AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
+    return 0x0;
+  }
+  
+  // get the cluster to re-clusterize
+  AliMUONVCluster* cluster = fESDInterface->GetCluster(iTrack,iCluster);
+  if (!cluster) return 0x0;
+  
+  // check if digits exist
+  if (cluster->GetNDigits() == 0) {
+    AliError(Form("no digit attached to cluster #%d in track %d",iCluster,iTrack));
+    return 0x0;
+  }
+  
+  // create the cluster store
+  AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
+  if (!clusterStore) return 0x0;
+  
+  // re-clusterize
+  TIter next(fESDInterface->CreateDigitIterator(iTrack, iCluster));
+  fClusterServer->UseDigits(next);
+  fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea());
+  
+  return clusterStore;
+}
+
+//_____________________________________________________________________________
+AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t clusterId)
+{
+  /// re-clusterize cluster "clusterId"
+  /// several new clusters may be reconstructed
+  /// it is the responsability of the user to delete the returned store
+  
+  if (!fESDInterface) {
+    AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
+    return 0x0;
+  }
+  
+  // get the cluster to re-clusterize
+  AliMUONVCluster* cluster = fESDInterface->FindCluster(clusterId);
+  if (!cluster) return 0x0;
+  
+  // check if digits exist
+  if (cluster->GetNDigits() == 0) {
+    AliError(Form("no digit attached to cluster %d",clusterId));
+    return 0x0;
+  }
+  
+  // create the cluster store
+  AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
+  if (!clusterStore) return 0x0;
+  
+  // re-clusterize
+  TIter next(fESDInterface->CreateDigitIteratorInCluster(clusterId));
+  fClusterServer->UseDigits(next);
+  fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea());
+  
+  return clusterStore;
+}
+
+//_____________________________________________________________________________
+void AliMUONRefitter::CreateGeometryTransformer()
+{
+  /// Create geometry transformer (local<->global)
+  /// and load geometry data
+  fGeometryTransformer = new AliMUONGeometryTransformer();
+  fGeometryTransformer->LoadGeometryData();
+}
+
+//_____________________________________________________________________________
+void AliMUONRefitter::CreateClusterServer(AliMUONGeometryTransformer& transformer)
+{
+  /// Create cluster server
+  AliMUONVClusterFinder* clusterFinder = CreateClusterFinder();
+  fClusterServer = clusterFinder ? new AliMUONSimpleClusterServer(*clusterFinder,transformer) : 0x0;
+}
+
+//_____________________________________________________________________________
+AliMUONVClusterFinder* AliMUONRefitter::CreateClusterFinder()
+{
+  /// Create a given cluster finder instance
+  AliMUONVClusterFinder* clusterFinder;
+  Option_t *opt = AliMUONReconstructor::GetRecoParam()->GetClusteringMode();
+  
+  if      (strstr(opt,"PRECLUSTERV2")) clusterFinder = new AliMUONPreClusterFinderV2;
+  else if (strstr(opt,"PRECLUSTERV3")) clusterFinder = new AliMUONPreClusterFinderV3;
+  else if (strstr(opt,"PRECLUSTER"))   clusterFinder = new AliMUONPreClusterFinder;
+  else if (strstr(opt,"COG"))          clusterFinder = new AliMUONClusterFinderCOG(new AliMUONPreClusterFinder);
+  else if (strstr(opt,"SIMPLEFITV3"))  clusterFinder = new AliMUONClusterFinderSimpleFit(new AliMUONClusterFinderCOG(new AliMUONPreClusterFinderV3));
+  else if (strstr(opt,"SIMPLEFIT"))    clusterFinder = new AliMUONClusterFinderSimpleFit(new AliMUONClusterFinderCOG(new AliMUONPreClusterFinder));
+  else if (strstr(opt,"MLEM:DRAW"))    clusterFinder = new AliMUONClusterFinderMLEM(kTRUE,new AliMUONPreClusterFinder);
+  else if (strstr(opt,"MLEMV3"))       clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinderV3);
+  else if (strstr(opt,"MLEMV2"))       clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinderV2);
+  else if (strstr(opt,"MLEM"))         clusterFinder = new AliMUONClusterFinderMLEM(kFALSE,new AliMUONPreClusterFinder);
+  else clusterFinder = 0x0;
+  
+  if (clusterFinder) {
+    AliInfo(Form("Will use %s for clusterizing",clusterFinder->ClassName()));
+  } else {
+    AliError(Form("clustering mode \"%s\" does not exist",opt));
+  }
+  
+  return clusterFinder;
+}
+
+//_____________________________________________________________________________
+void AliMUONRefitter::CreateTrackReconstructor(AliMUONVClusterServer& clusterServer)
+{
+  /// Create track reconstructor, depending on tracking mode set in RecoParam
+  Option_t *opt = AliMUONReconstructor::GetRecoParam()->GetTrackingMode();
+  
+  if      (strstr(opt,"ORIGINAL")) fTracker = new AliMUONTrackReconstructor(clusterServer);
+  else if (strstr(opt,"KALMAN"))   fTracker = new AliMUONTrackReconstructorK(clusterServer);
+  else fTracker = 0x0;
+  
+  if (fTracker) {
+    AliInfo(Form("Will use %s for tracking",fTracker->ClassName()));
+  } else {
+    AliError(Form("tracking mode \"%s\" does not exist",opt));
+  }
+}
+
+//_____________________________________________________________________________
+void AliMUONRefitter::AddClusterToTracks(const AliMUONVClusterStore &clusterStore, AliMUONVTrackStore &trackStore)
+{
+  /// add clusters to each of the given tracks
+  /// duplicate the tracks if there are several clusters and add one cluster per copy
+  
+  // create new track store if there are more than 1 cluster to add per track
+  Int_t nClusters = clusterStore.GetSize();
+  if (nClusters < 1) return;
+  
+  AliMUONTrackParam dummyParam;
+  AliMUONTrack *currentTrack, *track;
+  AliMUONVCluster *newCluster;
+  Int_t nTracks = trackStore.GetSize();
+  Int_t iTrack = 0;
+  Int_t iCluster = 0;
+  
+  // loop over existing tracks to add the cluster(s)
+  TIter nextTrack(trackStore.CreateIterator());
+  while ((currentTrack = static_cast<AliMUONTrack*>(nextTrack())) && (iTrack < nTracks)) {
+    
+    iTrack++;
+    
+    // add the new cluster(s) to the tracks
+    // duplicate the tracks if there are several clusters
+    // the loop after loading the last cluster which is added to the current track
+    iCluster = 0;
+    TIter nextCluster(clusterStore.CreateIterator());
+    while ((newCluster = static_cast<AliMUONVCluster*>(nextCluster())) && (iCluster < nClusters - 1)) {
+      
+      iCluster++;
+      
+      // add a copy of the current track to the store
+      track = trackStore.Add(AliMUONTrack(*currentTrack));
+      
+      // only set Z parameter to avoid error in AddTrackParamAtCluster()
+      // the rest will be recomputed during refit
+      dummyParam.SetZ(newCluster->GetZ());
+      
+      // add new cluster to the new track
+      track->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
+      
+    }
+    
+    // only set Z parameter to avoid error in AddTrackParamAtCluster()
+    // the rest will be recomputed during refit
+    dummyParam.SetZ(newCluster->GetZ());
+    
+    // add new cluster to the current track
+    currentTrack->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
+    
+  }
+  
+}
+
diff --git a/MUON/AliMUONRefitter.h b/MUON/AliMUONRefitter.h
new file mode 100644 (file)
index 0000000..7e98234
--- /dev/null
@@ -0,0 +1,76 @@
+#ifndef ALIMUONREFITTER_H
+#define ALIMUONREFITTER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice                               */
+
+// $Id$
+
+/// \ingroup rec
+/// \class AliMUONRefitter
+/// \brief class to refit the ESD clusters/tracks
+/// 
+//  Author Philippe Pillot
+
+#include <TObject.h>
+
+class AliMUONGeometryTransformer;
+class AliMUONVClusterFinder;
+class AliMUONVClusterServer;
+class AliMUONVTrackReconstructor;
+class AliMUONESDInterface;
+class AliMUONVClusterStore;
+class AliMUONVTrackStore;
+class AliMUONTrack;
+
+class AliMUONRefitter : public TObject
+{
+public:
+  
+  AliMUONRefitter();
+  virtual ~AliMUONRefitter();
+  
+  void Connect(AliMUONESDInterface* esdInterface) {fESDInterface = esdInterface;}
+  
+  // re-reconstruct all tracks (clusters) in the ESD event
+  AliMUONVTrackStore* ReconstructFromDigits();
+  AliMUONVTrackStore* ReconstructFromClusters();
+  
+  // refit a particular track in the ESD event
+  AliMUONTrack* RetrackFromDigits(Int_t iTrack);
+  AliMUONTrack* RetrackFromClusters(Int_t iTrack);
+  
+  // re-clusterize a particular cluster in the ESD event
+  AliMUONVClusterStore* ReClusterize(Int_t iTrack, Int_t iCluster);
+  AliMUONVClusterStore* ReClusterize(UInt_t clusterId);
+  
+  
+protected:
+  
+  AliMUONRefitter (const AliMUONRefitter&); ///< copy constructor
+  AliMUONRefitter& operator=(const AliMUONRefitter&); ///< assignment operator
+  
+  
+private:
+  
+  void                   CreateGeometryTransformer();
+  void                   CreateClusterServer(AliMUONGeometryTransformer& transformer);
+  AliMUONVClusterFinder* CreateClusterFinder();
+  void                   CreateTrackReconstructor(AliMUONVClusterServer& clusterServer);
+  
+  void AddClusterToTracks(const AliMUONVClusterStore &localClusterStore, AliMUONVTrackStore &trackStore);
+  
+  
+private:
+    
+  AliMUONGeometryTransformer* fGeometryTransformer; ///< geometry transformer (owner)
+  AliMUONVClusterServer*      fClusterServer;       ///< clusterizer (owner)
+  AliMUONVTrackReconstructor* fTracker;             ///< tracker (owner)
+  AliMUONESDInterface*        fESDInterface;        ///< container of MUON tracks/clusters/digits (not owner)
+  
+  
+  ClassDef(AliMUONRefitter,0)
+};
+
+#endif
+
index 75be873..853df2b 100644 (file)
 #include "AliMpSegmentation.h"
 #include "AliMpVPadIterator.h"
 #include "AliMpVSegmentation.h"
+#include "AliESDMuonPad.h"
 #include "AliLog.h"
 #include <float.h>
 #include <Riostream.h>
 #include <TClonesArray.h>
 #include <TString.h>
 
+
 /// \class AliMUONSimpleClusterServer
 ///
 /// Implementation of AliMUONVClusterServer interface
@@ -72,6 +74,9 @@ AliMUONSimpleClusterServer::AliMUONSimpleClusterServer(AliMUONVClusterFinder& cl
 {
     /// Ctor
     
+    fPads[0] = new AliMpExMap(true);
+    fPads[1] = new AliMpExMap(true);
+    
     AliMpDEIterator it;
     
     it.First();
@@ -144,6 +149,7 @@ AliMUONSimpleClusterServer::AliMUONSimpleClusterServer(AliMUONVClusterFinder& cl
 AliMUONSimpleClusterServer::~AliMUONSimpleClusterServer()
 {
   /// Dtor
+  delete &fClusterFinder;
   delete fPads[0];
   delete fPads[1];
   delete fDEAreas;
@@ -328,19 +334,14 @@ AliMUONSimpleClusterServer::PadArray(Int_t detElemId, Int_t cathode) const
 
 //_____________________________________________________________________________
 void 
-AliMUONSimpleClusterServer::UseDigitStore(const AliMUONVDigitStore& digitStore)
+AliMUONSimpleClusterServer::UseDigits(TIter& next)
 {
   /// Convert digitStore into two arrays of AliMUONPads
 
-  delete fPads[0];
-  delete fPads[1];
-  
-  fPads[0] = new AliMpExMap(true);
-  fPads[1] = new AliMpExMap(true);
+  fPads[0]->Clear();
+  fPads[1]->Clear();
   
-  TIter next(digitStore.CreateIterator());
   AliMUONVDigit* d;
-  
   while ( ( d = static_cast<AliMUONVDigit*>(next()) ) )
   {
     if ( ! d->Charge() > 0 ) continue; // skip void digits.
index 33be883..042da88 100644 (file)
@@ -18,8 +18,9 @@
 
 class AliMUONVClusterFinder;
 class AliMUONGeometryTransformer;
-class TClonesArray;
 class AliMpExMap;
+class AliESDMuonPad;
+class TClonesArray;
 
 class AliMUONSimpleClusterServer : public AliMUONVClusterServer
 {
@@ -33,7 +34,7 @@ public:
                    AliMUONVClusterStore& clusterStore,
                    const AliMpArea& area);
   
-  void UseDigitStore(const AliMUONVDigitStore& digitStore);
+  void UseDigits(TIter& next);
   
   void Print(Option_t* opt="") const;
   
@@ -50,8 +51,8 @@ private:
   TClonesArray* PadArray(Int_t detElemId, Int_t cathode) const;
   
 private:
-  AliMUONVClusterFinder& fClusterFinder; //!< the cluster finder
-  const AliMUONGeometryTransformer& fTransformer; //!< the geometry transformer
+  AliMUONVClusterFinder& fClusterFinder; //!< the cluster finder (owner)
+  const AliMUONGeometryTransformer& fTransformer; //!< the geometry transformer (not owner)
   AliMpExMap* fPads[2]; ///< map of TClonesArray of AliMUONPads
   AliMpExMap* fDEAreas; ///< map of detection element areas in global coordinates
   
index 4b03011..2570f74 100644 (file)
@@ -1156,3 +1156,40 @@ void AliMUONTrackReconstructor::FinalizeTrack(AliMUONTrack &track)
   if (!track.IsImproved()) track.UpdateCovTrackParamAtCluster();
 }
 
+//__________________________________________________________________________
+Bool_t AliMUONTrackReconstructor::RefitTrack(AliMUONTrack &track)
+{
+  /// re-fit the given track
+  
+  // check validity of the track
+  if (!track.IsValid()) {
+    AliWarning("the track does not contain enough clusters --> unable to refit");
+    return kFALSE;
+  }
+  
+  // reset the seed (i.e. parameters at first cluster) before fitting
+  AliMUONTrackParam* firstTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First();
+  if (firstTrackParam->GetInverseBendingMomentum() == 0.) {
+    AliWarning("track parameters at first chamber are not initialized --> unable to refit");
+    return kFALSE;
+  }
+  
+  // compute track parameters at each cluster from parameters at the first one
+  // necessary to compute multiple scattering effect during refitting
+  track.UpdateTrackParamAtCluster();
+  
+  // Re-fit the track:
+  // Take into account the multiple scattering
+  // Calculate the track parameter covariance matrix
+  Fit(track, kTRUE, kFALSE, kTRUE);
+  
+  // Improve the reconstructed tracks if required
+  if (AliMUONReconstructor::GetRecoParam()->ImproveTracks()) ImproveTrack(track);
+  
+  // Fill AliMUONTrack data members
+  FinalizeTrack(track);
+  
+  return kTRUE;
+  
+}
+
index beece1c..d4b63e8 100644 (file)
@@ -24,6 +24,8 @@ class AliMUONTrackReconstructor : public AliMUONVTrackReconstructor
   AliMUONTrackReconstructor(AliMUONVClusterServer& clusterServer); // default Constructor
   virtual ~AliMUONTrackReconstructor(); // Destructor
 
+  virtual Bool_t RefitTrack(AliMUONTrack &track);
+
 
  protected:
 
index ef2fd44..b996cc4 100644 (file)
@@ -1203,3 +1203,27 @@ void AliMUONTrackReconstructorK::FinalizeTrack(AliMUONTrack &track)
   
 }
 
+  //__________________________________________________________________________
+Bool_t AliMUONTrackReconstructorK::RefitTrack(AliMUONTrack &track)
+{
+  /// re-fit the given track
+  
+  // check validity of the track
+  if (!track.IsValid()) {
+    AliWarning("the track does not contain enough clusters --> unable to refit");
+    return kFALSE;
+  }
+  
+  // re-compute track parameters and covariances using Kalman filter
+  RetraceTrack(track,kTRUE);
+  
+  // Improve the reconstructed tracks if required
+  if (AliMUONReconstructor::GetRecoParam()->ImproveTracks()) ImproveTrack(track);
+  
+  // Fill AliMUONTrack data members
+  FinalizeTrack(track);
+  
+  return kTRUE;
+  
+}
+
index de8b4ee..93d71ab 100644 (file)
@@ -22,6 +22,8 @@ class AliMUONTrackReconstructorK : public AliMUONVTrackReconstructor
   
   AliMUONTrackReconstructorK(AliMUONVClusterServer& clusterServer); // default Constructor
   virtual ~AliMUONTrackReconstructorK(); // Destructor
+  
+  virtual Bool_t RefitTrack(AliMUONTrack &track);
 
 
  protected:
index 26d64d6..45c299b 100644 (file)
@@ -86,7 +86,8 @@ AliMUONTracker::AliMUONTracker(AliMUONVClusterServer& clusterServer,
   if (fTransformer && fDigitMaker)
     fTrackHitPatternMaker = new AliMUONTrackHitPattern(*fTransformer,*fDigitMaker);
   
-  fClusterServer.UseDigitStore(fDigitStore);
+  TIter next(fDigitStore.CreateIterator());
+  fClusterServer.UseDigits(next);
 }
 
 //_____________________________________________________________________________
index e79972d..3de623c 100644 (file)
 // Author Laurent Aphecetche, Subatech
 
 #ifndef ROOT_TObject
-#  include "TObject.h"
+#include "TObject.h"
 #endif
 
 class AliMUONVClusterStore;
-class AliMUONVDigitStore;
 class AliMpArea;
+class TIter;
 
 class AliMUONVClusterServer : public TObject
 {
@@ -32,7 +32,7 @@ public:
                            const AliMpArea& area) = 0;
   
   /// Use digits from the given digitstore to perform our job.
-  virtual void UseDigitStore(const AliMUONVDigitStore& digitStore) = 0;
+  virtual void UseDigits(TIter& next) = 0;
   
   ClassDef(AliMUONVClusterServer,1) // Cluster server interface
 };
index e2fad85..73264be 100644 (file)
@@ -49,6 +49,9 @@ class AliMUONVTrackReconstructor : public TObject {
                                  const AliMUONVTriggerStore& triggerStore,
                                  const AliMUONTrackHitPattern& trackHitPattern);
   
+  /// re-fit the given track
+  virtual Bool_t RefitTrack(AliMUONTrack &track) = 0;
+  
   
  protected:
 
diff --git a/MUON/MUONRefit.C b/MUON/MUONRefit.C
new file mode 100644 (file)
index 0000000..f3503bf
--- /dev/null
@@ -0,0 +1,244 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+/// \ingroup macros
+/// \file MUONRefit.C
+/// \brief Macro for refitting ESD tracks from ESD pads
+///
+/// \author Philippe Pillot, SUBATECH
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TStopwatch.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TString.h>
+#include <Riostream.h>
+#include <TGeoManager.h>
+#include <TRandom.h>
+#include <TROOT.h>
+
+// STEER includes
+#include "AliMagFMaps.h"
+#include "AliTracker.h"
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliRecoParam.h"
+#include "AliCDBManager.h"
+#include "AliGeomManager.h"
+
+// MUON includes
+#include "AliMpCDB.h"
+#include "AliMUONRecoParam.h"
+#include "AliMUONESDInterface.h"
+#include "AliMUONRefitter.h"
+#include "AliMUONVDigit.h"
+#include "AliMUONVCluster.h"
+#include "AliMUONVClusterStore.h"
+#include "AliMUONTrack.h"
+#include "AliMUONVTrackStore.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONTrackExtrap.h"
+#endif
+
+const Int_t printLevel = 1;
+
+void Prepare();
+TTree* GetESDTree(TFile *esdFile);
+
+//-----------------------------------------------------------------------
+void MUONRefit(Int_t nevents = -1, const char* esdFileNameIn = "AliESDs.root", const char* esdFileNameOut = "AliESDs_New.root")
+{
+  /// refit ESD tracks from ESD pads (i.e. re-clusterized the attached ESD clusters);
+  /// reset the charge of the digit using their raw charge before refitting;
+  /// compare results with original ESD tracks; 
+  /// write results in a new ESD file
+  
+  // prepare the refitting
+  gRandom->SetSeed(1);
+  Prepare();
+  AliMUONESDInterface esdInterface;
+  AliMUONRefitter refitter;
+  refitter.Connect(&esdInterface);
+  
+  // open the ESD file and tree
+  TFile* esdFile = TFile::Open(esdFileNameIn);
+  TTree* esdTree = GetESDTree(esdFile);
+  
+  // create the ESD output tree
+  gROOT->cd();
+  TTree* newESDTree = esdTree->CloneTree(0);
+  
+  // connect ESD event to the ESD tree
+  AliESDEvent* esd = new AliESDEvent();
+  esd->ReadFromTree(esdTree);
+
+  // timer start...
+  TStopwatch timer;
+  
+  // Loop over ESD events
+  if (nevents > 0) nevents = TMath::Min(nevents,(Int_t)esdTree->GetEntries());
+  else nevents = (Int_t)esdTree->GetEntries();
+  for (Int_t iEvent = 0; iEvent < nevents; iEvent++) {
+    if (printLevel>0) cout<<endl<<"            ****************event #"<<iEvent+1<<"****************"<<endl;
+    
+    //----------------------------------------------//
+    // -------------- process event --------------- //
+    //----------------------------------------------//
+    // get the ESD of current event
+    esdTree->GetEvent(iEvent);
+    if (!esd) {
+      Error("CheckESD", "no ESD object found for event %d", iEvent);
+      return;
+    }
+    Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
+    if (nTracks < 1) continue;
+    
+    // load the current event
+    esdInterface.LoadEvent(*esd);
+    
+    // loop over digit to modify their charge
+    AliMUONVDigit *digit;
+    TIter next(esdInterface.CreateDigitIterator());
+    while ((digit = static_cast<AliMUONVDigit*>(next()))) digit->SetCharge(digit->ADC());
+    
+    // refit the tracks from digits
+    AliMUONVTrackStore* newTrackStore = refitter.ReconstructFromDigits();
+    
+    //----------------------------------------------//
+    // ------ fill new ESD and print results ------ //
+    //----------------------------------------------//
+    // loop over the list of ESD tracks
+    TClonesArray *esdTracks = (TClonesArray*) esd->FindListObject("MuonTracks");
+    for (Int_t iTrack = 0; iTrack <  nTracks; iTrack++) {
+      
+      // get the ESD track
+      AliESDMuonTrack* esdTrack = (AliESDMuonTrack*) esdTracks->UncheckedAt(iTrack);
+      // get the corresponding MUON track
+      AliMUONTrack* track = esdInterface.GetTrackFast(iTrack);
+      // Find the corresponding re-fitted MUON track
+      AliMUONTrack* newTrack = (AliMUONTrack*) newTrackStore->FindObject(esdTrack->GetUniqueID());
+      
+      // replace the content of the current ESD track or remove it if the refitting has failed
+      if (newTrack) {
+       Double_t vertex[3] = {esdTrack->GetNonBendingCoor(), esdTrack->GetBendingCoor(), esdTrack->GetZ()};
+       AliMUONESDInterface::MUONToESD(*newTrack, *esdTrack, vertex, esdInterface.GetDigits());
+      } else {
+       esdTracks->Remove(esdTrack);
+      }
+      
+      // print initial and re-fitted track parameters at first cluster if any
+      if (printLevel>0) {
+       cout<<"            ----------------track #"<<iTrack+1<<"----------------"<<endl;
+       cout<<"before refit:"<<endl;
+       AliMUONTrackParam *param = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First();
+       param->Print("FULL");
+       if (printLevel>1) param->GetCovariances().Print();
+       if (!newTrack) continue;
+       cout<<"after refit:"<<endl;
+       param = (AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First();
+       param->Print("FULL");
+       if (printLevel>1) param->GetCovariances().Print();
+       cout<<"            ----------------------------------------"<<endl;
+      }
+      
+    }
+    
+    // free memory
+    delete newTrackStore;
+    
+    // fill new ESD tree with new tracks
+    esdTracks->Compress();
+    newESDTree->Fill();
+    
+    if (printLevel>0) cout<<"            ****************************************"<<endl;
+  }
+  
+  // ...timer stop
+  timer.Stop();
+  
+  // write output ESD tree
+  TFile* newESDFile = TFile::Open(esdFileNameOut, "RECREATE");
+  newESDFile->SetCompressionLevel(2);
+  newESDTree->Write();
+  newESDFile->Close();
+  
+  // free memory
+  esdFile->Close();
+  delete newESDTree;
+  delete esd;
+  
+  cout<<endl<<"time to refit: R:"<<timer.RealTime()<<" C:"<<timer.CpuTime()<<endl<<endl;
+}
+
+//-----------------------------------------------------------------------
+void Prepare()
+{
+  /// Set the geometry, the magnetic field, the mapping and the reconstruction parameters
+  
+  // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
+  if (!gGeoManager) {
+    AliGeomManager::LoadGeometry("geometry.root");
+    if (!gGeoManager) {
+      Error("MUONRefit", "getting geometry from file %s failed", "generated/galice.root");
+      return;
+    }
+  }
+  
+  // set mag field
+  printf("Loading field map...\n");
+  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
+  AliTracker::SetFieldMap(field, kFALSE);
+  AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
+  
+  // Load mapping
+  AliCDBManager* man = AliCDBManager::Instance();
+  man->SetDefaultStorage("local://$ALICE_ROOT");
+  man->SetRun(0);
+  if ( ! AliMpCDB::LoadDDLStore() ) {
+    Error("MUONRefit","Could not access mapping from OCDB !");
+    exit(-1);
+  }
+  
+  // set reconstruction parameters
+  AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLowFluxParam();
+  muonRecoParam->ImproveTracks(kFALSE);
+  muonRecoParam->Print("FULL");
+  AliRecoParam::Instance()->RegisterRecoParam(muonRecoParam);
+  
+}
+
+//-----------------------------------------------------------------------
+TTree* GetESDTree(TFile *esdFile)
+{
+  /// Check that the file is properly open
+  /// Return pointer to the ESD Tree
+  
+  if (!esdFile || !esdFile->IsOpen()) {
+    Error("GetESDTree", "opening ESD file failed");
+    exit(-1);
+  }
+  
+  TTree* tree = (TTree*) esdFile->Get("esdTree");
+  if (!tree) {
+    Error("GetESDTree", "no ESD tree found");
+    exit(-1);
+  }
+  
+  return tree;
+  
+}
+
index e7978a7..b7de853 100644 (file)
@@ -36,6 +36,7 @@
 #pragma link C++ class AliMUONClusterFinderMLEM+;
 #pragma link C++ class AliMUONClusterSplitterMLEM+;
 #pragma link C++ class AliMUONTrackHitPattern+;
+#pragma link C++ class AliMUONRefitter+;
 
 #pragma link C++ class AliMUONVClusterStore+;
 #pragma link C++ class AliMUONClusterStoreV1+;
index bc35c1a..1f95e36 100644 (file)
@@ -293,7 +293,8 @@ void TestRecPoints(TString baseDir=".", TString outDir=".", Float_t adcCut = 10.
         digitStoreTrackCut.Add(*mDigit, AliMUONVDigitStore::kDeny);
       } // loop on digits
       
-      clusterServer->UseDigitStore(digitStoreTrackCut);
+      TIter nextDigitTrackCut(digitStoreTrackCut->CreateIterator());
+      clusterServer->UseDigits(nextDigitTrackCut);
       
       for (Int_t ch = firstChamber; ch <= lastChamber; ++ch ){
         clusterServer->Clusterize(ch, clusterStore, AliMpArea());
index 9b42109..1359734 100644 (file)
@@ -41,6 +41,7 @@ SRCS:= AliMUONReconstructor.cxx \
        AliMUONQAChecker.cxx \
        AliMUONClusterFinderPeakCOG.cxx \
        AliMUONClusterFinderPeakFit.cxx \
+       AliMUONRefitter.cxx \
        AliMUONESDInterface.cxx
 
 HDRS:= $(SRCS:.cxx=.h)