MUON/AliMUONRefitter.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // $Id$
  17
  18 #include <cstdlib>
  19 #include "AliMUONRefitter.h"
  20 #include "AliMUONGeometryTransformer.h"
  21 #include "AliMUONClusterFinderCOG.h"
  22 #include "AliMUONClusterFinderMLEM.h"
  23 #include "AliMUONClusterFinderSimpleFit.h"
  24 #include "AliMUONPreClusterFinder.h"
  25 #include "AliMUONPreClusterFinderV2.h"
  26 #include "AliMUONPreClusterFinderV3.h"
  27 #include "AliMUONSimpleClusterServer.h"
  28 #include "AliMUONReconstructor.h"
  29 #include "AliMUONTrackReconstructor.h"
  30 #include "AliMUONTrackReconstructorK.h"
  31 #include "AliMUONRecoParam.h"
  32 #include "AliMUONESDInterface.h"
  33 #include "AliMUONVClusterStore.h"
  34 #include "AliMUONVTrackStore.h"
  35 #include "AliMUONTrack.h"
  36 #include "AliMUONTracker.h"
  37 #include "AliMUONTrackParam.h"
  38 #include "AliLog.h"
  39
  40 //-----------------------------------------------------------------------------
  41 /// \class AliMUONRefitter
  42 ///
  43 /// This class has been developped to simplify the re-reconstruction of the MUON tracks
  44 /// stored into ESD with different recoParams and/or after having re-calibrated the digits.
  45 /// It creates new MUON object from ESD objects given as input (through the ESDInterface) then:
  46 ///
  47 /// - re-clusterize the ESD clusters using the attached ESD pads
  48 ///   (several new clusters can be reconstructed per ESD cluster)
  49 /// - re-fit the ESD tracks using the attached ESD clusters
  50 /// - reconstruct the ESD tracks from ESD pads (i.e. re-clusterized the attached clusters)
  51 ///
  52 /// note:
  53 /// - connexion between an ESD cluster and corresponding MUON clusters from re-clustering
  54 ///   can be made through the detection element ID
  55 /// - connexion between an ESD track and the corresponding refitted MUON track
  56 ///   can be made through their unique ID
  57 ///
  58 /// \author Philippe Pillot
  59 //-----------------------------------------------------------------------------
  60
  61 /// \cond CLASSIMP
  62 ClassImp(AliMUONRefitter)
  63 /// \endcond
  64
  65 //_____________________________________________________________________________
  66 AliMUONRefitter::AliMUONRefitter(const AliMUONRecoParam* recoParam)
  67 : TObject(),
  68   fkRecoParam(recoParam),
  69   fkESDInterface(0x0),
  70   fGeometryTransformer(0x0),
  71   fClusterServer(0x0),
  72   fTracker(0x0)
  73 {
  74   /// Default constructor
  75   CreateGeometryTransformer();
  76   CreateClusterServer(*fGeometryTransformer);
  77   if (fClusterServer) fTracker = AliMUONTracker::CreateTrackReconstructor(recoParam,fClusterServer);
  78   if (!fClusterServer || !fTracker) {
  79     AliFatal("refitter initialization failed");
  80     exit(-1);
  81   }
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliMUONRefitter::~AliMUONRefitter()
  86 {
  87   /// Destructor
  88   delete fGeometryTransformer;
  89   delete fClusterServer;
  90   delete fTracker;
  91 }
  92
  93 //_____________________________________________________________________________
  94 AliMUONVTrackStore* AliMUONRefitter::ReconstructFromDigits()
  95 {
  96   /// re-reconstruct all tracks and attached clusters from the digits
  97   /// it is the responsability of the user to delete the returned store
  98
  99   if (!fkESDInterface) {
 100     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 101     return 0x0;
 102   }
 103
 104   // prepare new track(s)
 105   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
 106   if (!newTrackStore) return 0x0;
 107
 108   // loop over tracks and refit them (create new tracks)
 109   AliMUONTrack *track;
 110   TIter next(fkESDInterface->CreateTrackIterator());
 111   while ((track = static_cast<AliMUONTrack*>(next()))) {
 112     AliMUONTrack *newTrack = RetrackFromDigits(*track);
 113     if (newTrack) newTrackStore->Add(newTrack);
 114     delete newTrack;
 115   }
 116
 117   return newTrackStore;
 118 }
 119
 120 //_____________________________________________________________________________
 121 AliMUONVTrackStore* AliMUONRefitter::ReconstructFromClusters()
 122 {
 123   /// refit all tracks from the attached clusters
 124   /// it is the responsability of the user to delete the returned store
 125
 126   if (!fkESDInterface) {
 127     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 128     return 0x0;
 129   }
 130
 131   // prepare new track(s)
 132   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
 133   if (!newTrackStore) return 0x0;
 134
 135   // loop over tracks and refit them (create new tracks)
 136   AliMUONTrack *track;
 137   TIter next(fkESDInterface->CreateTrackIterator());
 138   while ((track = static_cast<AliMUONTrack*>(next()))) {
 139     AliMUONTrack* newTrack = newTrackStore->Add(*track);
 140     if (!fTracker->RefitTrack(*newTrack)) newTrackStore->Remove(*newTrack);
 141   }
 142
 143   return newTrackStore;
 144 }
 145
 146 //_____________________________________________________________________________
 147 AliMUONTrack* AliMUONRefitter::RetrackFromDigits(UInt_t trackId)
 148 {
 149   /// refit track "trackId" from the digits (i.e. re-clusterized the attached clusters)
 150   /// it is the responsability of the user to delete the returned track
 151
 152   if (!fkESDInterface) {
 153     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 154     return 0x0;
 155   }
 156
 157   // get the track to refit
 158   AliMUONTrack* track = fkESDInterface->FindTrack(trackId);
 159
 160   return track ? RetrackFromDigits(*track) : 0x0;
 161 }
 162
 163 //_____________________________________________________________________________
 164 AliMUONTrack* AliMUONRefitter::RetrackFromClusters(UInt_t trackId)
 165 {
 166   /// refit track "trackId" from the clusters (i.e. do not re-clusterize)
 167   /// it is the responsability of the user to delete the returned track
 168
 169   if (!fkESDInterface) {
 170     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 171     return 0x0;
 172   }
 173
 174   // get the track to refit
 175   AliMUONTrack* track = fkESDInterface->FindTrack(trackId);
 176   if (!track) return 0x0;
 177
 178   // refit the track (create a new one)
 179   AliMUONTrack* newTrack = new AliMUONTrack(*track);
 180   if (!fTracker->RefitTrack(*newTrack)) {
 181     delete newTrack;
 182     return 0x0;
 183   }
 184
 185   return newTrack;
 186 }
 187
 188 //_____________________________________________________________________________
 189 AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t trackId, UInt_t clusterId)
 190 {
 191   /// re-clusterize cluster numbered "clusterId" in track "trackId"
 192   /// several new clusters may be reconstructed
 193   /// it is the responsability of the user to delete the returned store
 194
 195   if (!fkESDInterface) {
 196     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 197     return 0x0;
 198   }
 199
 200   // get the cluster to re-clusterize
 201   AliMUONVCluster* cluster = fkESDInterface->FindCluster(trackId,clusterId);
 202   if (!cluster) return 0x0;
 203
 204   // check if digits exist
 205   if (cluster->GetNDigits() == 0) {
 206     AliError(Form("no digit attached to cluster #%d in track %d",clusterId,trackId));
 207     return 0x0;
 208   }
 209
 210   // create the cluster store
 211   AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
 212   if (!clusterStore) return 0x0;
 213
 214   // re-clusterize
 215   TIter next(fkESDInterface->CreateDigitIterator(trackId, clusterId));
 216   fClusterServer->UseDigits(next,fkESDInterface->GetDigits());
 217   fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea(),fkRecoParam);
 218
 219   return clusterStore;
 220 }
 221
 222 //_____________________________________________________________________________
 223 AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t clusterId)
 224 {
 225   /// re-clusterize cluster "clusterId"
 226   /// several new clusters may be reconstructed
 227   /// it is the responsability of the user to delete the returned store
 228
 229   if (!fkESDInterface) {
 230     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 231     return 0x0;
 232   }
 233
 234   // get the cluster to re-clusterize
 235   AliMUONVCluster* cluster = fkESDInterface->FindCluster(clusterId);
 236   if (!cluster) return 0x0;
 237
 238   // check if digits exist
 239   if (cluster->GetNDigits() == 0) {
 240     AliError(Form("no digit attached to cluster %d",clusterId));
 241     return 0x0;
 242   }
 243
 244   // create the cluster store
 245   AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
 246   if (!clusterStore) return 0x0;
 247
 248   // re-clusterize
 249   TIter next(fkESDInterface->CreateDigitIteratorInCluster(clusterId));
 250   fClusterServer->UseDigits(next,fkESDInterface->GetDigits());
 251   fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea(),fkRecoParam);
 252
 253   return clusterStore;
 254 }
 255
 256 //_____________________________________________________________________________
 257 void AliMUONRefitter::CreateGeometryTransformer()
 258 {
 259   /// Create geometry transformer (local<->global)
 260   /// and load geometry data
 261   fGeometryTransformer = new AliMUONGeometryTransformer();
 262   fGeometryTransformer->LoadGeometryData();
 263 }
 264
 265 //_____________________________________________________________________________
 266 void AliMUONRefitter::CreateClusterServer(AliMUONGeometryTransformer& transformer)
 267 {
 268   /// Create cluster server
 269   AliMUONVClusterFinder* clusterFinder = AliMUONReconstructor::CreateClusterFinder(fkRecoParam->GetClusteringMode());
 270   fClusterServer = clusterFinder ? new AliMUONSimpleClusterServer(clusterFinder,transformer) : 0x0;
 271 }
 272
 273 //_____________________________________________________________________________
 274 AliMUONTrack* AliMUONRefitter::RetrackFromDigits(const AliMUONTrack& track)
 275 {
 276   /// refit the given track from the digits (i.e. re-clusterized the attached clusters):
 277   /// several new clusters may be reconstructed per initial ESD cluster:
 278   /// -> all the combinations of clusters are considered to build the new tracks
 279   /// -> return the best track (largest number of clusters or best chi2 in case of equality)
 280
 281   // check if digits exist
 282   UInt_t trackId = track.GetUniqueID();
 283   if (fkESDInterface->GetNDigits(trackId) == 0) {
 284     AliError(Form("no digit attached to track #%d",trackId));
 285     return 0x0;
 286   }
 287
 288   // prepare new track(s)
 289   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
 290   if (!newTrackStore) return 0x0;
 291   newTrackStore->Add(track)->Clear("C");
 292
 293   // prepare new cluster store
 294   AliMUONVClusterStore* newClusterStore = AliMUONESDInterface::NewClusterStore();
 295   if (!newClusterStore) {
 296     delete newTrackStore;
 297     return 0x0;
 298   }
 299
 300   // loop over clusters, re-clusterize and build new tracks
 301   AliMUONVCluster* cluster;
 302   TIter nextCluster(fkESDInterface->CreateClusterIterator(trackId));
 303   while ((cluster = static_cast<AliMUONVCluster*>(nextCluster()))) {
 304
 305     // reset the new cluster store
 306     newClusterStore->Clear();
 307
 308     // re-clusterize current cluster
 309     TIter nextDigit(fkESDInterface->CreateDigitIterator(trackId, cluster->GetUniqueID()));
 310     fClusterServer->UseDigits(nextDigit,fkESDInterface->GetDigits());
 311     Int_t nNewClusters = fClusterServer->Clusterize(cluster->GetChamberId(),*newClusterStore,AliMpArea(),fkRecoParam);
 312
 313     // check that re-clusterizing succeeded
 314     if (nNewClusters == 0) {
 315       AliWarning(Form("refit gave no cluster (chamber %d)",cluster->GetChamberId()));
 316       AliInfo("initial ESD cluster:");
 317       cluster->Print("FULL");
 318       continue;
 319     }
 320
 321     // add the new cluster(s) to the tracks
 322     if (!AddClusterToTracks(*newClusterStore, *newTrackStore)) {
 323       delete newClusterStore;
 324       delete newTrackStore;
 325       return 0x0;
 326     }
 327
 328   }
 329
 330   if (newTrackStore->GetSize() > 1000) AliInfo(Form("%d tracks to refit... be patient!!",newTrackStore->GetSize()));
 331
 332   // refit the tracks and pick up the best one
 333   AliMUONTrack *currentTrack, *bestTrack = 0x0;
 334   Double_t currentChi2, bestChi2 = AliMUONTrack::MaxChi2();
 335   Int_t currentNCluster, bestNClusters = 0;
 336   TIter next(newTrackStore->CreateIterator());
 337   while ((currentTrack = static_cast<AliMUONTrack*>(next()))) {
 338
 339     // set the track parameters at first cluster if any (used as seed in original tracking)
 340     AliMUONTrackParam* param = (AliMUONTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
 341     if (param) *param = *((AliMUONTrackParam*) track.GetTrackParamAtCluster()->First());
 342
 343     // refit the track
 344     if (!fTracker->RefitTrack(*currentTrack)) continue;
 345
 346     // find best track (the one with the higher number of cluster or the best chi2 in case of equality)
 347     currentNCluster = currentTrack->GetNClusters();
 348     currentChi2 = currentTrack->GetGlobalChi2();
 349     if (currentNCluster > bestNClusters || (currentNCluster == bestNClusters && currentChi2 < bestChi2)) {
 350       bestTrack = currentTrack;
 351       bestNClusters = currentNCluster;
 352       bestChi2 = currentChi2;
 353     }
 354
 355   }
 356
 357   // copy best track and free memory
 358   AliMUONTrack* newTrack = bestTrack ? new AliMUONTrack(*bestTrack) : 0x0;
 359   delete newClusterStore;
 360   delete newTrackStore;
 361
 362   return newTrack;
 363 }
 364
 365 //_____________________________________________________________________________
 366 Bool_t AliMUONRefitter::AddClusterToTracks(const AliMUONVClusterStore &clusterStore, AliMUONVTrackStore &trackStore)
 367 {
 368   /// add clusters to each of the given tracks
 369   /// duplicate the tracks if there are several clusters and add one cluster per copy
 370
 371   // create new track store if there are more than 1 cluster to add per track
 372   Int_t nClusters = clusterStore.GetSize();
 373   if (nClusters < 1) return kTRUE;
 374
 375   // check if we will exceed the maximum allowed number of tracks
 376   if (nClusters * trackStore.GetSize() > fkRecoParam->GetMaxTrackCandidates()) {
 377     AliError(Form("Too many track candidates (%d tracks). Stop refitting.", nClusters * trackStore.GetSize()));
 378     return kFALSE;
 379   }
 380
 381   AliMUONTrackParam dummyParam;
 382   AliMUONTrack *currentTrack, *track;
 383   AliMUONVCluster *newCluster;
 384   Int_t nTracks = trackStore.GetSize();
 385   Int_t iTrack = 0;
 386   Int_t iCluster = 0;
 387
 388   // loop over existing tracks to add the cluster(s)
 389   TIter nextTrack(trackStore.CreateIterator());
 390   while ((currentTrack = static_cast<AliMUONTrack*>(nextTrack())) && (iTrack < nTracks)) {
 391
 392     iTrack++;
 393
 394     // add the new cluster(s) to the tracks
 395     // duplicate the tracks if there are several clusters
 396     // stop the loop after loading the last cluster which is added to the current track
 397     iCluster = 0;
 398     TIter nextCluster(clusterStore.CreateIterator());
 399     while ((newCluster = static_cast<AliMUONVCluster*>(nextCluster())) && (iCluster < nClusters - 1)) {
 400
 401       iCluster++;
 402
 403       // add a copy of the current track to the store
 404       track = trackStore.Add(AliMUONTrack(*currentTrack));
 405
 406       // only set Z parameter to avoid error in AddTrackParamAtCluster()
 407       // the rest will be recomputed during refit
 408       dummyParam.SetZ(newCluster->GetZ());
 409
 410       // add new cluster to the new track
 411       track->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
 412
 413     }
 414
 415     // only set Z parameter to avoid error in AddTrackParamAtCluster()
 416     // the rest will be recomputed during refit
 417     dummyParam.SetZ(newCluster->GetZ());
 418
 419     // add new cluster to the current track
 420     currentTrack->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
 421
 422   }
 423
 424   return kTRUE;
 425
 426 }
 427