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   nextClusterIndex(0)
  74 {
  75   /// Default constructor
  76   CreateGeometryTransformer();
  77   CreateClusterServer(*fGeometryTransformer);
  78   if (fClusterServer) fTracker = AliMUONTracker::CreateTrackReconstructor(recoParam,fClusterServer,fGeometryTransformer);
  79   if (!fClusterServer || !fTracker) {
  80     AliFatal("refitter initialization failed");
  81     exit(-1);
  82   }
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliMUONRefitter::~AliMUONRefitter()
  87 {
  88   /// Destructor
  89   delete fGeometryTransformer;
  90   delete fClusterServer;
  91   delete fTracker;
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliMUONVTrackStore* AliMUONRefitter::ReconstructFromDigits()
  96 {
  97   /// re-reconstruct all tracks and attached clusters from the digits
  98   /// it is the responsability of the user to delete the returned store
  99
 100   if (!fkESDInterface) {
 101     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 102     return 0x0;
 103   }
 104
 105   // prepare new track(s)
 106   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
 107   if (!newTrackStore) return 0x0;
 108
 109   // loop over tracks and refit them (create new tracks)
 110   AliMUONTrack *track;
 111   TIter next(fkESDInterface->CreateTrackIterator());
 112   while ((track = static_cast<AliMUONTrack*>(next()))) {
 113     AliMUONTrack *newTrack = RetrackFromDigits(*track);
 114     if (newTrack) newTrackStore->Add(newTrack);
 115     delete newTrack;
 116   }
 117
 118   return newTrackStore;
 119 }
 120
 121 //_____________________________________________________________________________
 122 AliMUONVTrackStore* AliMUONRefitter::ReconstructFromClusters()
 123 {
 124   /// refit all tracks from the attached clusters
 125   /// it is the responsability of the user to delete the returned store
 126
 127   if (!fkESDInterface) {
 128     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 129     return 0x0;
 130   }
 131
 132   // prepare new track(s)
 133   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
 134   if (!newTrackStore) return 0x0;
 135
 136   // loop over tracks and refit them (create new tracks)
 137   AliMUONTrack *track;
 138   TIter next(fkESDInterface->CreateTrackIterator());
 139   while ((track = static_cast<AliMUONTrack*>(next()))) {
 140     AliMUONTrack* newTrack = newTrackStore->Add(*track);
 141     if (!fTracker->RefitTrack(*newTrack)) newTrackStore->Remove(*newTrack);
 142   }
 143
 144   return newTrackStore;
 145 }
 146
 147 //_____________________________________________________________________________
 148 AliMUONTrack* AliMUONRefitter::RetrackFromDigits(UInt_t trackId)
 149 {
 150   /// refit track "trackId" from the digits (i.e. re-clusterized the attached clusters)
 151   /// it is the responsability of the user to delete the returned track
 152
 153   if (!fkESDInterface) {
 154     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 155     return 0x0;
 156   }
 157
 158   // get the track to refit
 159   AliMUONTrack* track = fkESDInterface->FindTrack(trackId);
 160
 161   return track ? RetrackFromDigits(*track) : 0x0;
 162 }
 163
 164 //_____________________________________________________________________________
 165 AliMUONTrack* AliMUONRefitter::RetrackFromClusters(UInt_t trackId)
 166 {
 167   /// refit track "trackId" from the clusters (i.e. do not re-clusterize)
 168   /// it is the responsability of the user to delete the returned track
 169
 170   if (!fkESDInterface) {
 171     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 172     return 0x0;
 173   }
 174
 175   // get the track to refit
 176   AliMUONTrack* track = fkESDInterface->FindTrack(trackId);
 177   if (!track) return 0x0;
 178
 179   // refit the track (create a new one)
 180   AliMUONTrack* newTrack = new AliMUONTrack(*track);
 181   if (!fTracker->RefitTrack(*newTrack)) {
 182     delete newTrack;
 183     return 0x0;
 184   }
 185
 186   return newTrack;
 187 }
 188
 189 //_____________________________________________________________________________
 190 AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t trackId, UInt_t clusterId)
 191 {
 192   /// re-clusterize cluster numbered "clusterId" in track "trackId"
 193   /// several new clusters may be reconstructed
 194   /// it is the responsability of the user to delete the returned store
 195
 196   if (!fkESDInterface) {
 197     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 198     return 0x0;
 199   }
 200
 201   // get the cluster to re-clusterize
 202   AliMUONVCluster* cluster = fkESDInterface->FindCluster(trackId,clusterId);
 203   if (!cluster) return 0x0;
 204
 205   // check if digits exist
 206   if (cluster->GetNDigits() == 0) {
 207     AliError(Form("no digit attached to cluster #%d in track %d",clusterId,trackId));
 208     return 0x0;
 209   }
 210
 211   // create the cluster store
 212   AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
 213   if (!clusterStore) return 0x0;
 214
 215   // re-clusterize
 216   TIter next(fkESDInterface->CreateDigitIterator(trackId, clusterId));
 217   fClusterServer->UseDigits(next,fkESDInterface->GetDigits());
 218   fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea(),fkRecoParam);
 219
 220   // set the uniqueID of the new clusters
 221   TIter nextCl(clusterStore->CreateIterator());
 222   AliMUONVCluster* newCluster = 0x0;
 223   while ((newCluster = static_cast<AliMUONVCluster*>(nextCl())))
 224     newCluster->SetUniqueID(AliMUONVCluster::BuildUniqueID(cluster->GetChamberId(), cluster->GetDetElemId(), nextClusterIndex++));
 225
 226   return clusterStore;
 227 }
 228
 229 //_____________________________________________________________________________
 230 AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t clusterId)
 231 {
 232   /// re-clusterize cluster "clusterId"
 233   /// several new clusters may be reconstructed
 234   /// it is the responsability of the user to delete the returned store
 235
 236   if (!fkESDInterface) {
 237     AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
 238     return 0x0;
 239   }
 240
 241   // get the cluster to re-clusterize
 242   AliMUONVCluster* cluster = fkESDInterface->FindCluster(clusterId);
 243   if (!cluster) return 0x0;
 244
 245   // check if digits exist
 246   if (cluster->GetNDigits() == 0) {
 247     AliError(Form("no digit attached to cluster %d",clusterId));
 248     return 0x0;
 249   }
 250
 251   // create the cluster store
 252   AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
 253   if (!clusterStore) return 0x0;
 254
 255   // re-clusterize
 256   TIter next(fkESDInterface->CreateDigitIteratorInCluster(clusterId));
 257   fClusterServer->UseDigits(next,fkESDInterface->GetDigits());
 258   fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea(),fkRecoParam);
 259
 260   // set the uniqueID of the new clusters
 261   TIter nextCl(clusterStore->CreateIterator());
 262   AliMUONVCluster* newCluster = 0x0;
 263   while ((newCluster = static_cast<AliMUONVCluster*>(nextCl())))
 264     newCluster->SetUniqueID(AliMUONVCluster::BuildUniqueID(cluster->GetChamberId(), cluster->GetDetElemId(), nextClusterIndex++));
 265
 266   return clusterStore;
 267 }
 268
 269 //_____________________________________________________________________________
 270 void AliMUONRefitter::CreateGeometryTransformer()
 271 {
 272   /// Create geometry transformer (local<->global)
 273   /// and load geometry data
 274   fGeometryTransformer = new AliMUONGeometryTransformer();
 275   fGeometryTransformer->LoadGeometryData();
 276 }
 277
 278 //_____________________________________________________________________________
 279 void AliMUONRefitter::CreateClusterServer(AliMUONGeometryTransformer& transformer)
 280 {
 281   /// Create cluster server
 282   AliMUONVClusterFinder* clusterFinder = AliMUONReconstructor::CreateClusterFinder(fkRecoParam->GetClusteringMode());
 283   fClusterServer = clusterFinder ? new AliMUONSimpleClusterServer(clusterFinder,transformer) : 0x0;
 284 }
 285
 286 //_____________________________________________________________________________
 287 AliMUONTrack* AliMUONRefitter::RetrackFromDigits(const AliMUONTrack& track)
 288 {
 289   /// refit the given track from the digits (i.e. re-clusterized the attached clusters):
 290   /// several new clusters may be reconstructed per initial ESD cluster:
 291   /// -> all the combinations of clusters are considered to build the new tracks
 292   /// -> return the best track (largest number of clusters or best chi2 in case of equality)
 293
 294   // check if digits exist
 295   UInt_t trackId = track.GetUniqueID();
 296   if (!fkESDInterface->DigitsStored(trackId)) {
 297     AliError(Form("no digit attached to track #%d",trackId));
 298     return 0x0;
 299   }
 300
 301   // prepare new track(s)
 302   AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
 303   if (!newTrackStore) return 0x0;
 304   newTrackStore->Add(track)->Clear("C");
 305
 306   // prepare new cluster store
 307   AliMUONVClusterStore* newClusterStore = AliMUONESDInterface::NewClusterStore();
 308   if (!newClusterStore) {
 309     delete newTrackStore;
 310     return 0x0;
 311   }
 312
 313   // loop over clusters, re-clusterize and build new tracks
 314   AliMUONVCluster* cluster;
 315   TIter nextCluster(fkESDInterface->CreateClusterIterator(trackId));
 316   while ((cluster = static_cast<AliMUONVCluster*>(nextCluster()))) {
 317
 318     // reset the new cluster store
 319     newClusterStore->Clear();
 320
 321     // re-clusterize current cluster
 322     TIter nextDigit(fkESDInterface->CreateDigitIterator(trackId, cluster->GetUniqueID()));
 323     fClusterServer->UseDigits(nextDigit,fkESDInterface->GetDigits());
 324     Int_t nNewClusters = fClusterServer->Clusterize(cluster->GetChamberId(),*newClusterStore,AliMpArea(),fkRecoParam);
 325
 326     // check that re-clusterizing succeeded
 327     if (nNewClusters == 0) {
 328       AliWarning(Form("refit gave no cluster (chamber %d)",cluster->GetChamberId()));
 329       AliInfo("initial ESD cluster:");
 330       cluster->Print("FULL");
 331       continue;
 332     }
 333
 334     // set the uniqueID of the new clusters
 335     TIter nextCl(newClusterStore->CreateIterator());
 336     AliMUONVCluster* newCluster = 0x0;
 337     while ((newCluster = static_cast<AliMUONVCluster*>(nextCl())))
 338       newCluster->SetUniqueID(AliMUONVCluster::BuildUniqueID(cluster->GetChamberId(), cluster->GetDetElemId(), nextClusterIndex++));
 339
 340     // add the new cluster(s) to the tracks
 341     if (!AddClusterToTracks(*newClusterStore, *newTrackStore)) {
 342       delete newClusterStore;
 343       delete newTrackStore;
 344       return 0x0;
 345     }
 346
 347   }
 348
 349   if (newTrackStore->GetSize() > 1000) AliInfo(Form("%d tracks to refit... be patient!!",newTrackStore->GetSize()));
 350
 351   // refit the tracks and pick up the best one
 352   AliMUONTrack *currentTrack, *bestTrack = 0x0;
 353   Double_t currentChi2, bestChi2 = AliMUONTrack::MaxChi2();
 354   Int_t currentNCluster, bestNClusters = 0;
 355   TIter next(newTrackStore->CreateIterator());
 356   while ((currentTrack = static_cast<AliMUONTrack*>(next()))) {
 357
 358     // set the track parameters at first cluster if any (used as seed in original tracking)
 359     AliMUONTrackParam* param = (AliMUONTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
 360     if (param) *param = *((AliMUONTrackParam*) track.GetTrackParamAtCluster()->First());
 361
 362     // refit the track
 363     if (!fTracker->RefitTrack(*currentTrack)) continue;
 364
 365     // check the status of the improvement if enabled
 366     if (fkRecoParam->ImproveTracks() && !currentTrack->IsImproved()) continue;
 367
 368     // find best track (the one with the highest number of cluster or the best chi2 in case of equality)
 369     currentNCluster = currentTrack->GetNClusters();
 370     currentChi2 = currentTrack->GetGlobalChi2();
 371     if (currentNCluster > bestNClusters || (currentNCluster == bestNClusters && currentChi2 < bestChi2)) {
 372       bestTrack = currentTrack;
 373       bestNClusters = currentNCluster;
 374       bestChi2 = currentChi2;
 375     }
 376
 377   }
 378
 379   // copy best track and free memory
 380   AliMUONTrack* newTrack = bestTrack ? new AliMUONTrack(*bestTrack) : 0x0;
 381   delete newClusterStore;
 382   delete newTrackStore;
 383
 384   return newTrack;
 385 }
 386
 387 //_____________________________________________________________________________
 388 Bool_t AliMUONRefitter::AddClusterToTracks(const AliMUONVClusterStore &clusterStore, AliMUONVTrackStore &trackStore)
 389 {
 390   /// add clusters to each of the given tracks
 391   /// duplicate the tracks if there are several clusters and add one cluster per copy
 392
 393   // create new track store if there are more than 1 cluster to add per track
 394   Int_t nClusters = clusterStore.GetSize();
 395   if (nClusters < 1) return kTRUE;
 396
 397   // check if we will exceed the maximum allowed number of tracks
 398   if (nClusters * trackStore.GetSize() > fkRecoParam->GetMaxTrackCandidates()) {
 399     AliError(Form("Too many track candidates (%d tracks). Stop refitting.", nClusters * trackStore.GetSize()));
 400     return kFALSE;
 401   }
 402
 403   AliMUONTrackParam dummyParam;
 404   AliMUONTrack *currentTrack, *track;
 405   AliMUONVCluster *newCluster;
 406   Int_t nTracks = trackStore.GetSize();
 407   Int_t iTrack = 0;
 408   Int_t iCluster = 0;
 409
 410   // loop over existing tracks to add the cluster(s)
 411   TIter nextTrack(trackStore.CreateIterator());
 412   while ((currentTrack = static_cast<AliMUONTrack*>(nextTrack())) && (iTrack < nTracks)) {
 413
 414     iTrack++;
 415
 416     // add the new cluster(s) to the tracks
 417     // duplicate the tracks if there are several clusters
 418     // stop the loop after loading the last cluster which is added to the current track
 419     iCluster = 0;
 420     TIter nextCluster(clusterStore.CreateIterator());
 421     while ((newCluster = static_cast<AliMUONVCluster*>(nextCluster())) && (iCluster < nClusters - 1)) {
 422
 423       iCluster++;
 424
 425       // add a copy of the current track to the store
 426       track = trackStore.Add(AliMUONTrack(*currentTrack));
 427
 428       // only set Z parameter to avoid error in AddTrackParamAtCluster()
 429       // the rest will be recomputed during refit
 430       dummyParam.SetZ(newCluster->GetZ());
 431
 432       // add new cluster to the new track
 433       track->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
 434
 435     }
 436
 437     // only set Z parameter to avoid error in AddTrackParamAtCluster()
 438     // the rest will be recomputed during refit
 439     dummyParam.SetZ(newCluster->GetZ());
 440
 441     // add new cluster to the current track
 442     currentTrack->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
 443
 444   }
 445
 446   return kTRUE;
 447
 448 }
 449