#include <cstdlib>
#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 "AliMUONReconstructor.h"
#include "AliMUONTrackReconstructor.h"
#include "AliMUONTrackReconstructorK.h"
#include "AliMUONRecoParam.h"
#include "AliMUONESDInterface.h"
#include "AliMUONVClusterStore.h"
#include "AliMUONVTrackStore.h"
#include "AliMUONTrack.h"
#include "AliMUONTracker.h"
#include "AliMUONTrackParam.h"
#include "AliLog.h"
ClassImp(AliMUONRefitter)
AliMUONRefitter::AliMUONRefitter(const AliMUONRecoParam* recoParam)
: TObject(),
fkRecoParam(recoParam),
fkESDInterface(0x0),
fGeometryTransformer(0x0),
fClusterServer(0x0),
fTracker(0x0),
nextClusterIndex(0)
{
CreateGeometryTransformer();
CreateClusterServer(*fGeometryTransformer);
if (fClusterServer) fTracker = AliMUONTracker::CreateTrackReconstructor(recoParam,fClusterServer,fGeometryTransformer);
if (!fClusterServer || !fTracker) {
AliFatal("refitter initialization failed");
exit(-1);
}
}
AliMUONRefitter::~AliMUONRefitter()
{
delete fGeometryTransformer;
delete fClusterServer;
delete fTracker;
}
AliMUONVTrackStore* AliMUONRefitter::ReconstructFromDigits()
{
if (!fkESDInterface) {
AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
return 0x0;
}
AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
if (!newTrackStore) return 0x0;
AliMUONTrack *track;
TIter next(fkESDInterface->CreateTrackIterator());
while ((track = static_cast<AliMUONTrack*>(next()))) {
AliMUONTrack *newTrack = RetrackFromDigits(*track);
if (newTrack) newTrackStore->Add(newTrack);
delete newTrack;
}
return newTrackStore;
}
AliMUONVTrackStore* AliMUONRefitter::ReconstructFromClusters()
{
if (!fkESDInterface) {
AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
return 0x0;
}
AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
if (!newTrackStore) return 0x0;
AliMUONTrack *track;
TIter next(fkESDInterface->CreateTrackIterator());
while ((track = static_cast<AliMUONTrack*>(next()))) {
AliMUONTrack* newTrack = newTrackStore->Add(*track);
if (!fTracker->RefitTrack(*newTrack)) newTrackStore->Remove(*newTrack);
}
return newTrackStore;
}
AliMUONTrack* AliMUONRefitter::RetrackFromDigits(UInt_t trackId)
{
if (!fkESDInterface) {
AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
return 0x0;
}
AliMUONTrack* track = fkESDInterface->FindTrack(trackId);
return track ? RetrackFromDigits(*track) : 0x0;
}
AliMUONTrack* AliMUONRefitter::RetrackFromClusters(UInt_t trackId)
{
if (!fkESDInterface) {
AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
return 0x0;
}
AliMUONTrack* track = fkESDInterface->FindTrack(trackId);
if (!track) return 0x0;
AliMUONTrack* newTrack = new AliMUONTrack(*track);
if (!fTracker->RefitTrack(*newTrack)) {
delete newTrack;
return 0x0;
}
return newTrack;
}
AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t trackId, UInt_t clusterId)
{
if (!fkESDInterface) {
AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
return 0x0;
}
AliMUONVCluster* cluster = fkESDInterface->FindCluster(trackId,clusterId);
if (!cluster) return 0x0;
if (cluster->GetNDigits() == 0) {
AliError(Form("no digit attached to cluster #%d in track %d",clusterId,trackId));
return 0x0;
}
AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
if (!clusterStore) return 0x0;
TIter next(fkESDInterface->CreateDigitIterator(trackId, clusterId));
fClusterServer->UseDigits(next,fkESDInterface->GetDigits());
fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea(),fkRecoParam);
TIter nextCl(clusterStore->CreateIterator());
AliMUONVCluster* newCluster = 0x0;
while ((newCluster = static_cast<AliMUONVCluster*>(nextCl())))
newCluster->SetUniqueID(AliMUONVCluster::BuildUniqueID(cluster->GetChamberId(), cluster->GetDetElemId(), nextClusterIndex++));
return clusterStore;
}
AliMUONVClusterStore* AliMUONRefitter::ReClusterize(UInt_t clusterId)
{
if (!fkESDInterface) {
AliError("the refitter must be connected to an ESDInterface containing the ESD event to reconstruct");
return 0x0;
}
AliMUONVCluster* cluster = fkESDInterface->FindCluster(clusterId);
if (!cluster) return 0x0;
if (cluster->GetNDigits() == 0) {
AliError(Form("no digit attached to cluster %d",clusterId));
return 0x0;
}
AliMUONVClusterStore* clusterStore = AliMUONESDInterface::NewClusterStore();
if (!clusterStore) return 0x0;
TIter next(fkESDInterface->CreateDigitIteratorInCluster(clusterId));
fClusterServer->UseDigits(next,fkESDInterface->GetDigits());
fClusterServer->Clusterize(cluster->GetChamberId(),*clusterStore,AliMpArea(),fkRecoParam);
TIter nextCl(clusterStore->CreateIterator());
AliMUONVCluster* newCluster = 0x0;
while ((newCluster = static_cast<AliMUONVCluster*>(nextCl())))
newCluster->SetUniqueID(AliMUONVCluster::BuildUniqueID(cluster->GetChamberId(), cluster->GetDetElemId(), nextClusterIndex++));
return clusterStore;
}
void AliMUONRefitter::CreateGeometryTransformer()
{
fGeometryTransformer = new AliMUONGeometryTransformer();
fGeometryTransformer->LoadGeometryData();
}
void AliMUONRefitter::CreateClusterServer(AliMUONGeometryTransformer& transformer)
{
AliMUONVClusterFinder* clusterFinder = AliMUONReconstructor::CreateClusterFinder(fkRecoParam->GetClusteringMode());
fClusterServer = clusterFinder ? new AliMUONSimpleClusterServer(clusterFinder,transformer) : 0x0;
}
AliMUONTrack* AliMUONRefitter::RetrackFromDigits(const AliMUONTrack& track)
{
UInt_t trackId = track.GetUniqueID();
if (!fkESDInterface->DigitsStored(trackId)) {
AliError(Form("no digit attached to track #%d",trackId));
return 0x0;
}
AliMUONVTrackStore* newTrackStore = AliMUONESDInterface::NewTrackStore();
if (!newTrackStore) return 0x0;
newTrackStore->Add(track)->Clear("C");
AliMUONVClusterStore* newClusterStore = AliMUONESDInterface::NewClusterStore();
if (!newClusterStore) {
delete newTrackStore;
return 0x0;
}
AliMUONVCluster* cluster;
TIter nextCluster(fkESDInterface->CreateClusterIterator(trackId));
while ((cluster = static_cast<AliMUONVCluster*>(nextCluster()))) {
newClusterStore->Clear();
TIter nextDigit(fkESDInterface->CreateDigitIterator(trackId, cluster->GetUniqueID()));
fClusterServer->UseDigits(nextDigit,fkESDInterface->GetDigits());
Int_t nNewClusters = fClusterServer->Clusterize(cluster->GetChamberId(),*newClusterStore,AliMpArea(),fkRecoParam);
if (nNewClusters == 0) {
AliWarning(Form("refit gave no cluster (chamber %d)",cluster->GetChamberId()));
AliInfo("initial ESD cluster:");
cluster->Print("FULL");
continue;
}
TIter nextCl(newClusterStore->CreateIterator());
AliMUONVCluster* newCluster = 0x0;
while ((newCluster = static_cast<AliMUONVCluster*>(nextCl())))
newCluster->SetUniqueID(AliMUONVCluster::BuildUniqueID(cluster->GetChamberId(), cluster->GetDetElemId(), nextClusterIndex++));
if (!AddClusterToTracks(*newClusterStore, *newTrackStore)) {
delete newClusterStore;
delete newTrackStore;
return 0x0;
}
}
if (newTrackStore->GetSize() > 1000) AliInfo(Form("%d tracks to refit... be patient!!",newTrackStore->GetSize()));
AliMUONTrack *currentTrack, *bestTrack = 0x0;
Double_t currentChi2, bestChi2 = AliMUONTrack::MaxChi2();
Int_t currentNCluster, bestNClusters = 0;
TIter next(newTrackStore->CreateIterator());
while ((currentTrack = static_cast<AliMUONTrack*>(next()))) {
AliMUONTrackParam* param = (AliMUONTrackParam*) currentTrack->GetTrackParamAtCluster()->First();
if (param) *param = *((AliMUONTrackParam*) track.GetTrackParamAtCluster()->First());
if (!fTracker->RefitTrack(*currentTrack)) continue;
if (fkRecoParam->ImproveTracks() && !currentTrack->IsImproved()) continue;
currentNCluster = currentTrack->GetNClusters();
currentChi2 = currentTrack->GetGlobalChi2();
if (currentNCluster > bestNClusters || (currentNCluster == bestNClusters && currentChi2 < bestChi2)) {
bestTrack = currentTrack;
bestNClusters = currentNCluster;
bestChi2 = currentChi2;
}
}
AliMUONTrack* newTrack = bestTrack ? new AliMUONTrack(*bestTrack) : 0x0;
delete newClusterStore;
delete newTrackStore;
return newTrack;
}
Bool_t AliMUONRefitter::AddClusterToTracks(const AliMUONVClusterStore &clusterStore, AliMUONVTrackStore &trackStore)
{
Int_t nClusters = clusterStore.GetSize();
if (nClusters < 1) return kTRUE;
if (nClusters * trackStore.GetSize() > fkRecoParam->GetMaxTrackCandidates()) {
AliError(Form("Too many track candidates (%d tracks). Stop refitting.", nClusters * trackStore.GetSize()));
return kFALSE;
}
AliMUONTrackParam dummyParam;
AliMUONTrack *currentTrack, *track;
AliMUONVCluster *newCluster;
Int_t nTracks = trackStore.GetSize();
Int_t iTrack = 0;
Int_t iCluster = 0;
TIter nextTrack(trackStore.CreateIterator());
while ((currentTrack = static_cast<AliMUONTrack*>(nextTrack())) && (iTrack < nTracks)) {
iTrack++;
iCluster = 0;
TIter nextCluster(clusterStore.CreateIterator());
while ((newCluster = static_cast<AliMUONVCluster*>(nextCluster())) && (iCluster < nClusters - 1)) {
iCluster++;
track = trackStore.Add(AliMUONTrack(*currentTrack));
dummyParam.SetZ(newCluster->GetZ());
track->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
}
dummyParam.SetZ(newCluster->GetZ());
currentTrack->AddTrackParamAtCluster(dummyParam, *newCluster, kTRUE);
}
return kTRUE;
}