G4ITStepProcessor.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// $Id: G4ITStepProcessor.cc 74551 2013-10-14 12:59:14Z gcosmo $
//
// Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr) 
//
// History:
// -----------
// 10 Oct 2011 M.Karamitros created
//
// -------------------------------------------------------------------

#include "G4ITStepProcessor.hh"
#include "G4UImanager.hh"
#include "G4ForceCondition.hh"
#include "G4GPILSelection.hh"
#include "G4ITTransportationManager.hh"
// #include "G4VSensitiveDetector.hh"    // Include from 'hits/digi'
#include "G4GeometryTolerance.hh"
#include "G4ParticleTable.hh"
#include "G4ITTrackingManager.hh"
#include "G4TrackingInformation.hh"
#include "G4IT.hh"
#include "G4ITNavigator.hh"             // Include from 'geometry'

#include "G4VITProcess.hh"
#include "G4VProcess.hh"
#include "G4ITTransportation.hh"

#include <iomanip>              // Include from 'system'
#include <vector>               // Include from 'system'

using namespace std;

static const size_t SizeOfSelectedDoItVector=100;
//static const size_t& gMaxNProcesses(G4VITProcess::GetMaxProcessIndex());

//____________________________________________________________________________________

G4ITStepProcessor::G4ITStepProcessor()
{
    verboseLevel = 0 ;
    //    fpUserSteppingAction = 0 ;
    fStoreTrajectory = 0;
    fpTrackingManager = 0;
    fpNavigator = 0;
    kCarTolerance = -1.;
    fInitialized = false;
    fPreviousTimeStep = DBL_MAX;
    CleanProcessor();
    ResetSecondaries();
}

G4ITStepProcessor::G4ITStepProcessorState::G4ITStepProcessorState() :
    G4ITStepProcessorState_Lock(),
//    fSelectedAtRestDoItVector (gMaxNProcesses,0),
//    fSelectedPostStepDoItVector (gMaxNProcesses,0)
    fSelectedAtRestDoItVector (G4VITProcess::GetMaxProcessIndex(),0),
    fSelectedPostStepDoItVector (G4VITProcess::GetMaxProcessIndex(),0)
{
    fPhysicalStep = -1.;
    fPreviousStepSize = -1.;

    fSafety = -1.;
    proposedSafety = -1.;
    endpointSafety = -1;

    fStepStatus = fUndefined;

    fTouchableHandle = 0;
}

// should not be used
G4ITStepProcessor::G4ITStepProcessorState::G4ITStepProcessorState(const G4ITStepProcessorState& ) :
    G4ITStepProcessorState_Lock(),
//    fSelectedAtRestDoItVector (gMaxNProcesses,0),
//    fSelectedPostStepDoItVector (gMaxNProcesses,0)
    fSelectedAtRestDoItVector (G4VITProcess::GetMaxProcessIndex(),0),
    fSelectedPostStepDoItVector (G4VITProcess::GetMaxProcessIndex(),0)
{
    fPhysicalStep = -1.;
    fPreviousStepSize = -1.;

    fSafety = -1.;
    proposedSafety = -1.;
    endpointSafety = -1;

    fStepStatus = fUndefined;

    fTouchableHandle = 0;
}

// should not be used
G4ITStepProcessor::G4ITStepProcessorState&  G4ITStepProcessor::G4ITStepProcessorState::operator=(const G4ITStepProcessorState& rhs)
{
    if(this == &rhs) return *this;

    fSelectedAtRestDoItVector.clear();
//    fSelectedAtRestDoItVector.resize(gMaxNProcesses,0);
    fSelectedAtRestDoItVector.resize(G4VITProcess::GetMaxProcessIndex(),0);
    fSelectedPostStepDoItVector.clear();
//    fSelectedPostStepDoItVector.resize(gMaxNProcesses,0);
    fSelectedPostStepDoItVector.resize(G4VITProcess::GetMaxProcessIndex(),0);

    fPhysicalStep = -1.;
    fPreviousStepSize = -1.;

    fSafety = -1.;
    proposedSafety = -1.;
    endpointSafety = -1;

    fStepStatus = fUndefined;

    fTouchableHandle = 0;
    return *this;
}
//____________________________________________________________________________________

G4ITStepProcessor::G4ITStepProcessorState::~G4ITStepProcessorState()
{;}
//____________________________________________________________________________________

void G4ITStepProcessor::ClearProcessInfo()
{
    std::map<const G4ParticleDefinition*, ProcessGeneralInfo*> ::iterator it;

    for(it = fProcessGeneralInfoMap.begin();it != fProcessGeneralInfoMap.end();it++)
    {
        if(it->second)
        {
            delete it->second;
            it->second = 0;
        }
    }

    fProcessGeneralInfoMap.clear();
}

//____________________________________________________________________________________

void G4ITStepProcessor::ForceReInitialization()
{
    fInitialized = false;
    ClearProcessInfo();
    Initialize();
}

//____________________________________________________________________________________

void G4ITStepProcessor::Initialize()
{
    CleanProcessor();
    if(fInitialized) return;
    //    ActiveOnlyITProcess();

    SetNavigator(G4ITTransportationManager::GetTransportationManager()
                 ->GetNavigatorForTracking());

    fPhysIntLength = DBL_MAX;
    kCarTolerance = 0.5*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();

    fInitialized = true;
}
//______________________________________________________________________________

G4ITStepProcessor::~G4ITStepProcessor()
{
    if(fpStep)
    {
        fpStep->DeleteSecondaryVector();
        delete fpStep;
    }

    if(fpSecondary)                      delete fpSecondary;
    ClearProcessInfo();
    G4ITTransportationManager::DeleteInstance();

    //    if(fpUserSteppingAction)             delete fpUserSteppingAction;
}
//______________________________________________________________________________
// should not be used
G4ITStepProcessor::G4ITStepProcessor(const G4ITStepProcessor& rhs)
{
    verboseLevel = rhs.verboseLevel ;
    fStoreTrajectory = rhs.fStoreTrajectory ;

    //    fpUserSteppingAction = 0 ;
    fpTrackingManager = 0;
    fpNavigator = 0;
    fInitialized = false;

    kCarTolerance = rhs.kCarTolerance;
    fInitialized = false;
    fPreviousTimeStep = DBL_MAX;

    CleanProcessor();
    ResetSecondaries();
}
//______________________________________________________________________________

G4ITStepProcessor& G4ITStepProcessor::operator=(const G4ITStepProcessor& rhs)
{
    if (this == &rhs) return *this; // handle self assignment
    //assignment operator
    return *this;
}
// ******************************************************************

void G4ITStepProcessor::ActiveOnlyITProcess()
{
    // Method not used for the time being
#ifdef debug
    G4cout<<"G4ITStepProcessor::CloneProcesses: is called"<<G4endl;
#endif

    G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
    G4ParticleTable::G4PTblDicIterator* theParticleIterator = theParticleTable->GetIterator();

    theParticleIterator->reset();
    // TODO : Ne faire la boucle que sur les IT **** !!!
    while( (*theParticleIterator)() )
    {
        G4ParticleDefinition* particle = theParticleIterator->value();
        G4ProcessManager* pm= particle->GetProcessManager();

        if(!pm)
        {
            G4cerr << "ERROR - G4ITStepProcessor::GetProcessNumber()" << G4endl
                   << "        ProcessManager is NULL for particle = "
                   << particle->GetParticleName() << ", PDG_code = "
                   << particle->GetPDGEncoding() << G4endl;
            G4Exception("G4ITStepProcessor::GetProcessNumber()", "ITStepProcessor0001",
                        FatalException, "Process Manager is not found.");
            return;
        }

        ActiveOnlyITProcess(pm);
    }
}
// ******************************************************************

void G4ITStepProcessor::ActiveOnlyITProcess(G4ProcessManager* processManager)
{
    // Method not used for the time being
    G4ProcessVector* processVector = processManager->GetProcessList();

    G4VITProcess* itProcess = 0 ;
    for(int i = 0 ; i < processVector->size() ; i++)
    {
        G4VProcess* base_process = (*processVector)[i];
        itProcess = dynamic_cast<G4VITProcess*>(base_process);

        if(!itProcess)
        {
            processManager->SetProcessActivation(base_process, false);
        }
    }
}
// ******************************************************************
void G4ITStepProcessor::SetupGeneralProcessInfo(G4ParticleDefinition* particle,
                                                G4ProcessManager* pm)
{

#ifdef debug
    G4cout<<"G4ITStepProcessor::GetProcessNumber: is called track"<<G4endl;
#endif
    if(!pm)
    {
        G4cerr << "ERROR - G4SteppingManager::GetProcessNumber()" << G4endl
               << "        ProcessManager is NULL for particle = "
               << particle->GetParticleName() << ", PDG_code = "
               << particle->GetPDGEncoding() << G4endl;
        G4Exception("G4SteppingManager::GetProcessNumber()", "ITStepProcessor0002",
                    FatalException, "Process Manager is not found.");
        return;
    }

    std::map<const G4ParticleDefinition*, ProcessGeneralInfo*>::iterator it = fProcessGeneralInfoMap.find(particle);
    if(it != fProcessGeneralInfoMap.end())
    {
        G4Exception("G4SteppingManager::SetupGeneralProcessInfo()", "ITStepProcessor0003",
                    FatalException, "Process info already registered.");
        return;
    }

    // here used as temporary
    fpProcessInfo = new ProcessGeneralInfo();

    // AtRestDoits
    fpProcessInfo->MAXofAtRestLoops =        pm->GetAtRestProcessVector()->entries();
    fpProcessInfo->fpAtRestDoItVector =       pm->GetAtRestProcessVector(typeDoIt);
    fpProcessInfo->fpAtRestGetPhysIntVector = pm->GetAtRestProcessVector(typeGPIL);
#ifdef debug
    G4cout << "G4ITStepProcessor::GetProcessNumber: #ofAtRest="
           << fpProcessInfo->MAXofAtRestLoops << G4endl;
#endif

    // AlongStepDoits
    fpProcessInfo->MAXofAlongStepLoops = pm->GetAlongStepProcessVector()->entries();
    fpProcessInfo->fpAlongStepDoItVector = pm->GetAlongStepProcessVector(typeDoIt);
    fpProcessInfo->fpAlongStepGetPhysIntVector = pm->GetAlongStepProcessVector(typeGPIL);
#ifdef debug
    G4cout << "G4ITStepProcessor::GetProcessNumber:#ofAlongStp="
           << fpProcessInfo->MAXofAlongStepLoops << G4endl;
#endif

    // PostStepDoits
    fpProcessInfo->MAXofPostStepLoops = pm->GetPostStepProcessVector()->entries();
    fpProcessInfo->fpPostStepDoItVector = pm->GetPostStepProcessVector(typeDoIt);
    fpProcessInfo->fpPostStepGetPhysIntVector = pm->GetPostStepProcessVector(typeGPIL);
#ifdef debug
    G4cout << "G4ITStepProcessor::GetProcessNumber: #ofPostStep="
           << fpProcessInfo->MAXofPostStepLoops << G4endl;
#endif

    if (SizeOfSelectedDoItVector<fpProcessInfo->MAXofAtRestLoops    ||
            SizeOfSelectedDoItVector<fpProcessInfo->MAXofAlongStepLoops ||
            SizeOfSelectedDoItVector<fpProcessInfo->MAXofPostStepLoops  )
    {
        G4cerr << "ERROR - G4ITStepProcessor::GetProcessNumber()" << G4endl
               << "        SizeOfSelectedDoItVector= " << SizeOfSelectedDoItVector
               << " ; is smaller then one of MAXofAtRestLoops= "
               << fpProcessInfo->MAXofAtRestLoops << G4endl
               << "        or MAXofAlongStepLoops= " << fpProcessInfo->MAXofAlongStepLoops
               << " or MAXofPostStepLoops= " << fpProcessInfo->MAXofPostStepLoops << G4endl;
        G4Exception("G4ITStepProcessor::GetProcessNumber()",
                    "ITStepProcessor0004", FatalException,
                    "The array size is smaller than the actual No of processes.");
    }

    if(!fpProcessInfo->fpAtRestDoItVector       &&
            !fpProcessInfo->fpAlongStepDoItVector    &&
            !fpProcessInfo->fpPostStepDoItVector)
    {
        G4ExceptionDescription exceptionDescription ;
        exceptionDescription << "No DoIt process found " ;
        G4Exception("G4ITStepProcessor::DoStepping","ITStepProcessor0005",
                    FatalErrorInArgument,exceptionDescription);
        return ;
    }

    if(fpProcessInfo->fpAlongStepGetPhysIntVector && fpProcessInfo->MAXofAlongStepLoops>0)
    {
        fpProcessInfo->fpTransportation = dynamic_cast<G4ITTransportation*>
                ((*fpProcessInfo->fpAlongStepGetPhysIntVector)[fpProcessInfo->MAXofAlongStepLoops-1]);

        if(fpProcessInfo->fpTransportation == 0)
        {
            G4ExceptionDescription exceptionDescription ;
            exceptionDescription << "No transportation process found " ;
            G4Exception("G4ITStepProcessor::SetupGeneralProcessInfo","ITStepProcessor0006",
                        FatalErrorInArgument,exceptionDescription);
        }
    }
    fProcessGeneralInfoMap[particle] = fpProcessInfo;
    //    fpProcessInfo = 0;
}

// ******************************************************************

void G4ITStepProcessor::SetTrack(G4Track* track)
{
    fpTrack = track ;
    if(fpTrack)
    {
        fpITrack = GetIT(fpTrack) ;
        fpStep = const_cast<G4Step*>(fpTrack -> GetStep());

        if(fpITrack)
        {
            fpTrackingInfo = fpITrack->GetTrackingInfo() ;
        }
        else
        {
            fpTrackingInfo = 0;
            G4cerr << "Track ID : " << fpTrack->GetTrackID() << G4endl;

            G4ExceptionDescription exceptionDescription ("No IT pointer was attached to the track you try to process.");
            G4Exception("G4ITStepProcessor::SetTrack","ITStepProcessor0007",
                        FatalErrorInArgument,exceptionDescription);
        }
    }
    else
    {
        fpITrack = 0;
        fpStep = 0 ;
    }
}
//______________________________________________________________________________

void G4ITStepProcessor::GetProcessInfo()
{
    G4ParticleDefinition* particle = fpTrack->GetDefinition();
    std::map<const G4ParticleDefinition*, ProcessGeneralInfo*>::iterator it = fProcessGeneralInfoMap.find(particle);

    if(it == fProcessGeneralInfoMap.end())
    {
        SetupGeneralProcessInfo(particle,fpTrack->GetDefinition()->GetProcessManager());
        if(fpProcessInfo == 0)
        {
            G4ExceptionDescription exceptionDescription ("...");
            G4Exception("G4ITStepProcessor::GetProcessNumber","ITStepProcessor0008",
                        FatalErrorInArgument,exceptionDescription);
            return;
        }
    }
    else
    {
        fpProcessInfo = it->second;
    }
}
//______________________________________________________________________________

void G4ITStepProcessor::SetupMembers()
{
    fpSecondary      = fpStep->GetfSecondary();
    fpPreStepPoint   = fpStep->GetPreStepPoint();
    fpPostStepPoint  = fpStep->GetPostStepPoint();

    fpState = (G4ITStepProcessorState*) fpITrack->GetTrackingInfo()->GetStepProcessorState();

    GetProcessInfo();
    ResetSecondaries();
}
//______________________________________________________________________________

void G4ITStepProcessor::ResetSecondaries()
{
    // Reset the secondary particles
    fN2ndariesAtRestDoIt    = 0;
    fN2ndariesAlongStepDoIt = 0;
    fN2ndariesPostStepDoIt  = 0;
}
//______________________________________________________________________________

void G4ITStepProcessor::GetAtRestIL()
{
    // Select the rest process which has the shortest time before
    // it is invoked. In rest processes, GPIL()
    // returns the time before a process occurs.
    G4double lifeTime (DBL_MAX), shortestLifeTime (DBL_MAX);

    fAtRestDoItProcTriggered = 0;
    shortestLifeTime = DBL_MAX;

    unsigned int NofInactiveProc=0;

    for( size_t ri=0 ; ri < fpProcessInfo->MAXofAtRestLoops ; ri++ )
    {
        fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpAtRestGetPhysIntVector)[ri];
        if (fpCurrentProcess== 0)
        {
            (fpState->fSelectedAtRestDoItVector)[ri] = InActivated;
            NofInactiveProc++;
            continue;
        }   // NULL means the process is inactivated by a user on fly.

        fCondition=NotForced;
        fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
        lifeTime = fpCurrentProcess->AtRestGPIL( *fpTrack, &fCondition );
        fpCurrentProcess->SetProcessState(0);

        if(fCondition==Forced)
        {
            (fpState->fSelectedAtRestDoItVector)[ri] = Forced;
        }
        else
        {
            (fpState->fSelectedAtRestDoItVector)[ri] = InActivated;
            if(lifeTime < shortestLifeTime )
            {
                shortestLifeTime = lifeTime;
                fAtRestDoItProcTriggered =  G4int(ri);
                (fpState->fSelectedAtRestDoItVector)[fAtRestDoItProcTriggered] = NotForced;
            }
        }
    }

    fTimeStep = shortestLifeTime ;

    // at least one process is necessary to destroy the particle
    // exit with warning
    if(NofInactiveProc==fpProcessInfo->MAXofAtRestLoops)
    {
        G4cerr << "ERROR - G4ITStepProcessor::InvokeAtRestDoItProcs()" << G4endl
               << "        No AtRestDoIt process is active!" << G4endl;
    }
}
//___________________________________________________________________________

void G4ITStepProcessor::DefinePhysicalStepLength(G4Track* track)
{
    SetTrack(track);
    DoDefinePhysicalStepLength();
}
//______________________________________________________________________________


void G4ITStepProcessor::SetInitialStep()
{
    // DEBUG
    //    G4cout << "SetInitialStep for : " << fpITrack-> GetName() << G4endl;
    //________________________________________________________
    // Initialize geometry


    if ( ! fpTrack->GetTouchableHandle())
    {
        G4ThreeVector direction= fpTrack->GetMomentumDirection();
        fpNavigator->LocateGlobalPointAndSetup( fpTrack->GetPosition(),
                                                &direction, false, false );
        fpState->fTouchableHandle = fpNavigator->CreateTouchableHistory();

        fpTrack->SetTouchableHandle( fpState->fTouchableHandle );
        fpTrack->SetNextTouchableHandle( fpState->fTouchableHandle );
    }
    else
    {
        fpState->fTouchableHandle = fpTrack->GetTouchableHandle();
        fpTrack->SetNextTouchableHandle( fpState->fTouchableHandle );
        G4VPhysicalVolume* oldTopVolume= fpTrack->GetTouchableHandle()->GetVolume();
        G4VPhysicalVolume* newTopVolume=
                fpNavigator->ResetHierarchyAndLocate( fpTrack->GetPosition(),
                                                      fpTrack->GetMomentumDirection(),
                                                      *((G4TouchableHistory*)fpTrack->GetTouchableHandle()()) );
        if(newTopVolume != oldTopVolume || oldTopVolume->GetRegularStructureId() == 1 )
        {
            fpState->fTouchableHandle = fpNavigator->CreateTouchableHistory();
            fpTrack->SetTouchableHandle( fpState->fTouchableHandle );
            fpTrack->SetNextTouchableHandle( fpState->fTouchableHandle );
        }
    }

    fpCurrentVolume = fpState->fTouchableHandle->GetVolume();

    //________________________________________________________
    // If the primary track has 'Suspend' or 'PostponeToNextEvent' state,
    // set the track state to 'Alive'.
    if( (fpTrack->GetTrackStatus()==fSuspend) ||
            (fpTrack->GetTrackStatus()==fPostponeToNextEvent) )
    {
        fpTrack->SetTrackStatus(fAlive);
    }

    // If the primary track has 'zero' kinetic energy, set the track
    // state to 'StopButAlive'.
    if(fpTrack->GetKineticEnergy() <= 0.0)
    {
        fpTrack->SetTrackStatus( fStopButAlive );
    }
    //________________________________________________________
    // Set vertex information of G4Track at here
    if ( fpTrack->GetCurrentStepNumber() == 0 )
    {
        fpTrack->SetVertexPosition( fpTrack->GetPosition() );
        fpTrack->SetVertexMomentumDirection( fpTrack->GetMomentumDirection() );
        fpTrack->SetVertexKineticEnergy( fpTrack->GetKineticEnergy() );
        fpTrack->SetLogicalVolumeAtVertex( fpTrack->GetVolume()->GetLogicalVolume() );
    }
    //________________________________________________________
    // If track is already outside the world boundary, kill it
    if( fpCurrentVolume==0 )
    {
        // If the track is a primary, stop processing
        if(fpTrack->GetParentID()==0)
        {
            G4cerr << "ERROR - G4ITStepProcessor::SetInitialStep()" << G4endl
                   << "        Primary particle starting at - "
                   << fpTrack->GetPosition()
                   << " - is outside of the world volume." << G4endl;
            G4Exception("G4ITStepProcessor::SetInitialStep()", "ITStepProcessor0011",
                        FatalException, "Primary vertex outside of the world!");
        }

        fpTrack->SetTrackStatus( fStopAndKill );
        G4cout << "WARNING - G4ITStepProcessor::SetInitialStep()" << G4endl
               << "          Initial track position is outside world! - "
               << fpTrack->GetPosition() << G4endl;
    }
    else{
        // Initial set up for attribues of 'Step'
        fpStep->InitializeStep( fpTrack );
    }


    if( fpTrack->GetTrackStatus() == fStopAndKill ) return ;

    fpTrackingManager->StartTracking(fpTrack);

    fpState->fStepStatus = fUndefined;
}
//______________________________________________________________________________

void G4ITStepProcessor::InitDefineStep()
{

    if(!fpStep)
    {

        // Create new Step and give it to the track
        fpStep = new G4Step();
        fpTrack->SetStep(fpStep);
        fpSecondary = fpStep->NewSecondaryVector();

        // Create new state and set it in the trackingInfo
        fpState = new G4ITStepProcessorState();
        fpITrack->GetTrackingInfo()->SetStepProcessorState((G4ITStepProcessorState_Lock*)fpState);

        SetupMembers();
        fpNavigator->NewNavigatorState();

        SetInitialStep();
    }
    else
    {
        SetupMembers();

        fpState->fPreviousStepSize = fpTrack->GetStepLength();
/***
        // Send G4Step information to Hit/Dig if the volume is sensitive
        fpCurrentVolume = fpStep->GetPreStepPoint()->GetPhysicalVolume();
        StepControlFlag =  fpStep->GetControlFlag();
        if( fpCurrentVolume != 0 && StepControlFlag != AvoidHitInvocation)
        {
            fpSensitive = fpStep->GetPreStepPoint()->
                    GetSensitiveDetector();

            //            if( fSensitive != 0 ) {
            //              fSensitive->Hit(fStep);
            //            }
        }
***/
        // Store last PostStepPoint to PreStepPoint, and swap current and next
        // volume information of G4Track. Reset total energy deposit in one Step.
        fpStep->CopyPostToPreStepPoint();
        fpStep->ResetTotalEnergyDeposit();

        //JA Set the volume before it is used (in DefineStepLength() for User Limit)
        fpCurrentVolume = fpStep->GetPreStepPoint()->GetPhysicalVolume();
/*
        G4cout << G4endl;
        G4cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!" << G4endl;
        G4cout << "PreStepPoint Volume : " << fpCurrentVolume->GetName() << G4endl;
        G4cout << "Track Touchable : " << fpTrack->GetTouchableHandle()->GetVolume()->GetName() << G4endl;
        G4cout << "Track NextTouchable : " << fpTrack->GetNextTouchableHandle()->GetVolume()->GetName() << G4endl;
*/
        // Reset the step's auxiliary points vector pointer
        fpStep->SetPointerToVectorOfAuxiliaryPoints(0);

        // Switch next touchable in track to current one
        fpTrack->SetTouchableHandle(fpTrack->GetNextTouchableHandle());
        fpState->fTouchableHandle = fpTrack->GetTouchableHandle();
        fpTrack->SetNextTouchableHandle( fpState->fTouchableHandle );
        G4VPhysicalVolume* oldTopVolume= fpTrack->GetTouchableHandle()->GetVolume();
        fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()->GetNavigatorState());

        G4VPhysicalVolume* newTopVolume=
                fpNavigator->ResetHierarchyAndLocate( fpTrack->GetPosition(),
                                                      fpTrack->GetMomentumDirection(),
                                                      *((G4TouchableHistory*)fpTrack->GetTouchableHandle()()) );

        //        G4cout << "New Top Volume : " << newTopVolume->GetName() << G4endl;

        if(newTopVolume != oldTopVolume || oldTopVolume->GetRegularStructureId() == 1 )
        {
            fpState->fTouchableHandle = fpNavigator->CreateTouchableHistory();
            fpTrack->SetTouchableHandle( fpState->fTouchableHandle );
            fpTrack->SetNextTouchableHandle( fpState->fTouchableHandle );
        }

        fpNavigator->SetNavigatorState(fpITrack->GetTrackingInfo()->GetNavigatorState());
    }
}

//______________________________________________________________________________


// ************************************************************************
//      Compute Interaction Length
// ************************************************************************
void G4ITStepProcessor::DoDefinePhysicalStepLength()
{

    InitDefineStep();

    G4TrackStatus trackStatus = fpTrack -> GetTrackStatus()  ;

    if(trackStatus == fStopAndKill)
    {
        return ;
    }

    if(trackStatus == fStopButAlive)
    {
        fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator->GetNavigatorState());
        fpNavigator->SetNavigatorState(0);
        return GetAtRestIL() ;
    }


    // Find minimum Step length and corresponding time
    // demanded by active disc./cont. processes

    // ReSet the counter etc.
    fpState->fPhysicalStep  = DBL_MAX;          // Initialize by a huge number
    fPhysIntLength = DBL_MAX;          // Initialize by a huge number

    double proposedTimeStep = DBL_MAX;
    G4VProcess* processWithPostStepGivenByTimeStep(0);

    // GPIL for PostStep
    fPostStepDoItProcTriggered = fpProcessInfo->MAXofPostStepLoops;
    fPostStepAtTimeDoItProcTriggered = fpProcessInfo->MAXofPostStepLoops;

    //    G4cout << "fpProcessInfo->MAXofPostStepLoops : " << fpProcessInfo->MAXofPostStepLoops
    //           << " mol : " << fpITrack -> GetName() << " id : " << fpTrack->GetTrackID()
    //           << G4endl;

    for(size_t np=0; np < fpProcessInfo->MAXofPostStepLoops; np++)
    {
        fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpPostStepGetPhysIntVector)[np];
        if (fpCurrentProcess== 0)
        {
            (fpState->fSelectedPostStepDoItVector)[np] = InActivated;
            continue;
        }   // NULL means the process is inactivated by a user on fly.

        fCondition=NotForced;
        fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));

        //        G4cout << "Is going to call : " << fpCurrentProcess -> GetProcessName() << G4endl;
        fPhysIntLength = fpCurrentProcess->
                PostStepGPIL( *fpTrack,
                              fpState->fPreviousStepSize,
                              &fCondition );
        fpCurrentProcess->SetProcessState(0);

        switch (fCondition)
        {
        case ExclusivelyForced: // Will need special treatment
            (fpState->fSelectedPostStepDoItVector)[np] = ExclusivelyForced;
            fpState->fStepStatus = fExclusivelyForcedProc;
            fpStep->GetPostStepPoint()
                    ->SetProcessDefinedStep(fpCurrentProcess);
            break;

        case Conditionally:
            //       (fpState->fSelectedPostStepDoItVector)[np] = Conditionally;
            G4Exception("G4ITStepProcessor::DefinePhysicalStepLength()", "ITStepProcessor0008",
                        FatalException, "This feature is no more supported");
            break;

        case Forced:
            (fpState->fSelectedPostStepDoItVector)[np] = Forced;
            break;

        case StronglyForced:
            (fpState->fSelectedPostStepDoItVector)[np] = StronglyForced;
            break;

        default:
            (fpState->fSelectedPostStepDoItVector)[np] = InActivated;
            break;
        }

        if (fCondition==ExclusivelyForced)
        {
            for(size_t nrest=np+1; nrest < fpProcessInfo->MAXofPostStepLoops; nrest++)
            {
                (fpState->fSelectedPostStepDoItVector)[nrest] = InActivated;
            }
            return;  // Please note the 'return' at here !!!
        }
        else
        {
            if(fPhysIntLength < fpState->fPhysicalStep )
            {
                // To avoid checking whether the process is actually
                // proposing a time step, the returned time steps are
                // negative (just for tagging)
                if(fpCurrentProcess->ProposesTimeStep())
                {
                    fPhysIntLength *= -1;
                    if(fPhysIntLength < proposedTimeStep)
                    {
                        proposedTimeStep = fPhysIntLength;
                        fPostStepAtTimeDoItProcTriggered = np;
                        processWithPostStepGivenByTimeStep = fpCurrentProcess;
                    }
                }
                else
                {
                    fpState->fPhysicalStep = fPhysIntLength;
                    fpState->fStepStatus = fPostStepDoItProc;
                    fPostStepDoItProcTriggered = G4int(np);
                    fpStep->GetPostStepPoint()
                            ->SetProcessDefinedStep(fpCurrentProcess);
                }
            }
        }
    }

    // GPIL for AlongStep
    fpState->proposedSafety = DBL_MAX;
    G4double safetyProposedToAndByProcess = fpState->proposedSafety;

    for(size_t kp=0; kp < fpProcessInfo->MAXofAlongStepLoops; kp++)
    {
        fpCurrentProcess = (G4VITProcess*) (*fpProcessInfo->fpAlongStepGetPhysIntVector)[kp];
        if (fpCurrentProcess== 0) continue;
        // NULL means the process is inactivated by a user on fly.

        fpCurrentProcess->SetProcessState(fpTrackingInfo->GetProcessState(fpCurrentProcess->GetProcessID()));
        fPhysIntLength = fpCurrentProcess-> AlongStepGPIL( *fpTrack,
                                                           fpState->fPreviousStepSize,
                                                           fpState->fPhysicalStep,
                                                           safetyProposedToAndByProcess,
                                                           &fGPILSelection );

        if(fPhysIntLength < fpState->fPhysicalStep)
        {
            fpState->fPhysicalStep      = fPhysIntLength;
            // Should save PS and TS in IT

            // Check if the process wants to be the GPIL winner. For example,
            // multi-scattering proposes Step limit, but won't be the winner.
            if(fGPILSelection==CandidateForSelection)
            {
                fpState->fStepStatus = fAlongStepDoItProc;
                fpStep->GetPostStepPoint()
                        ->SetProcessDefinedStep(fpCurrentProcess);
            }

            // Transportation is assumed to be the last process in the vector
            if(kp == fpProcessInfo->MAXofAlongStepLoops-1)
            {
                fpTransportation = dynamic_cast<G4ITTransportation*>(fpCurrentProcess);

                if(! fpTransportation)
                {
                    G4ExceptionDescription exceptionDescription ;
                    exceptionDescription << "No transportation process found " ;
                    G4Exception("G4ITStepProcessor::DoDefinePhysicalStepLength","ITStepProcessor0009",
                                FatalErrorInArgument,exceptionDescription);
                }

                fTimeStep               = fpTransportation->GetInteractionTimeLeft();


                if (fpTrack->GetNextVolume() != 0)
                    fpState->fStepStatus = fGeomBoundary;
                else
                    fpState->fStepStatus = fWorldBoundary;
            }
        }
        else
        {
            if(kp == fpProcessInfo->MAXofAlongStepLoops-1)
            {
                fpTransportation = dynamic_cast<G4ITTransportation*>(fpCurrentProcess);

                if(! fpTransportation)
                {
                    G4ExceptionDescription exceptionDescription ;
                    exceptionDescription << "No transportation process found " ;
                    G4Exception("G4ITStepProcessor::DoDefinePhysicalStepLength","ITStepProcessor0010",
                                FatalErrorInArgument,exceptionDescription);
                }

                fTimeStep               = fpTransportation->GetInteractionTimeLeft();
            }
        }

        if(proposedTimeStep < fTimeStep)
        {
            if(fPostStepAtTimeDoItProcTriggered<fpProcessInfo->MAXofPostStepLoops)
            {
                if ((fpState->fSelectedPostStepDoItVector)[fPostStepAtTimeDoItProcTriggered] ==
                        InActivated)
                {
                    (fpState->fSelectedPostStepDoItVector)[fPostStepAtTimeDoItProcTriggered] = NotForced;
                    //                    (fpState->fSelectedPostStepDoItVector)[fPostStepDoItProcTriggered] = InActivated;

                    fpState->fStepStatus = fPostStepDoItProc;
                    fpStep->GetPostStepPoint()->SetProcessDefinedStep(processWithPostStepGivenByTimeStep);

                    fTimeStep = proposedTimeStep;

                    fpTransportation->ComputeStep(*fpTrack,*fpStep,fTimeStep,fpState->fPhysicalStep);
                }
            }
        }
        else
        {
            if (fPostStepDoItProcTriggered<fpProcessInfo->MAXofPostStepLoops)
            {
                if ((fpState->fSelectedPostStepDoItVector)[fPostStepDoItProcTriggered] ==
                        InActivated)
                {
                    (fpState->fSelectedPostStepDoItVector)[fPostStepDoItProcTriggered] =
                            NotForced;
                }
            }
        }

        fpCurrentProcess->SetProcessState(0);

        // Make sure to check the safety, even if Step is not limited
        //  by this process.                      J. Apostolakis, June 20, 1998
        //
        if (safetyProposedToAndByProcess < fpState->proposedSafety)
            // proposedSafety keeps the smallest value:
            fpState->proposedSafety               = safetyProposedToAndByProcess;
        else
            // safetyProposedToAndByProcess always proposes a valid safety:
            safetyProposedToAndByProcess = fpState->proposedSafety;

    }

    fpITrack->GetTrackingInfo()->SetNavigatorState(fpNavigator->GetNavigatorState());
    fpNavigator->SetNavigatorState(0);
}

//______________________________________________________________________________