In This Package:

QueriableStepAction.cc

Go to the documentation of this file.

#include "QueriableStepAction.h"
#include "G4ToHistoryUtensils.h"
#include "G4HistoryUserTrackInfo.h"
#include "HistorianCustomRules.h"

#include "G4TouchableHistory.hh"
#include "G4TouchableHandle.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4StepPoint.hh"
#include "G4SteppingManager.hh"
#include "G4MaterialPropertiesTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4LossTableManager.hh"
#include "G4MaterialCutsCouple.hh"
#include "G4Gamma.hh"
#include "G4Electron.hh"

#include "RuleParser/ParameterDescription.h"
#include "RuleParser/PrescaleRule.h"

#include "G4DataHelpers/ITouchableToDetectorElement.h"
#include "DetDesc/IDetectorElement.h"
#include "DetDesc/GeometryInfoPlus.h"
#include "DetDesc/IPVolume.h"

#include <vector>


using namespace DayaBay;



QueriableStepAction::QueriableStepAction ( 
  const std::string& type   , 
  const std::string& name   , 
  const IInterface*  parent ) 
  : GiGaStepActionBase( type , name , parent )
  , m_DetectorElementSearchPath(1,"/dd/Structure")
  , mCurrentStep(0)
  , mCurrentTrack(0)
  , mCurrentStepPoint(0)
  , mTouchToDetElem(0)
{ 
  declareProperty("TouchableToDetelem",m_TouchableToDetelem_name = "TH2DE");
  declareProperty("DetectorElementSearchPath",m_DetectorElementSearchPath);

  m_IdParameterNames.push_back("PmtID");
  m_IdParameterNames.push_back("DetectorID");
  m_IdParameterNames.push_back("SiteID");
  declareProperty("IdParameterNames",m_IdParameterNames,
                  "Names of the user parameters attached to detector elements used to give ID numbers.");
  
  // Liang Zhan, Jan 18, 2009.
  declareProperty("BirksConstant1", m_BirksConstant1 = 6.5e-3*g/cm2/MeV, "Birks constant C1");
  declareProperty("BirksConstant2", m_BirksConstant2 = 2.1e-6*(g/cm2/MeV)*(g/cm2/MeV), "Birks constant C2");
}

StatusCode QueriableStepAction::initialize() 
{
  StatusCode sc =  GiGaStepActionBase::initialize();
  if(sc.isFailure()) return sc;

  std::string t2de_name = std::string("t2de_") + name();
  mTouchToDetElem = tool<ITouchableToDetectorElement>(
     m_TouchableToDetelem_name,
     name(),
     this);

  // Get the Neutron Capture tool which is updated in DsPhysConHadron.cc
  // to gain access to the neutron capture information.
  // --- Wei, Aug 15, 2008
  m_capinfo = tool<INeutronCaptureInfo>("G4DhNeutronCaptureInfoTool");

  verbose() << " QueriableStepAction::initialize() done" << endreq;

  return StatusCode::SUCCESS;
}

StatusCode QueriableStepAction::finalize() 
{
  return  GiGaStepActionBase::finalize();
}

void QueriableStepAction::Reset(const G4Step* step, 
                                const G4Track* track, 
                                const G4StepPoint* point)
                                
{
  mCurrentStep = step;
  mCurrentTrack = track;
  mCurrentStepPoint = point;
  
  mDetElement=0;
  mDetElementMatch=-1;
  mLogicVolume=0;
  mPhysVolume=0;
  mProcess=SimProcess();
  mQuenchedEnergy = -1e9;
  mUserTrackInfo = 0;
  mDetectorId = 0xFFFFFFF;

  verbose() << " The current track has been reset " << endreq;
}


StatusCode  QueriableStepAction::GetTrackParameterList(RuleParser::ParameterList& tplist)
{
  using namespace RuleParser;
  using std::string;
  
  tplist.add<double>(kPar_t,"time","t");
  tplist.add<double>(kPar_x,"x","global_x");
  tplist.add<double>(kPar_y,"y","global_y");
  tplist.add<double>(kPar_z,"z","global_z");
  tplist.add<double>(kPar_r,"r","radius","pos_r");
  tplist.add<double>(kPar_local_x,"lx","local_x","det_x");             
  tplist.add<double>(kPar_local_y,"ly","local_y","det_y");
  tplist.add<double>(kPar_local_z,"lz","local_z","det_z");
  tplist.add<double>(kPar_local_r,"lr","local_r","det_r");
  tplist.add<string>(kPar_LogicalVolumeName,"Volume","LogicalVolumeName","LogicalVolume","VolumeName");
  tplist.add<string>(kPar_MaterialName,"Material","MaterialName");
  tplist.add<string>(kPar_DetectorElementName,"DetectorElementName");
  tplist.add<double>(kPar_DetectorElementMatch,"Match","DetectorElementMatch");
  tplist.add<double>(kPar_NicheId,"NicheId","Niche");
  tplist.add<double>(kPar_DetectorId,"DetectorId");
  tplist.add<double>(kPar_SiteId,"SiteId");
  tplist.add<double>(kPar_Site,"Site");
  tplist.add<double>(kPar_AD,"AD","AdNumber");

  tplist.add<double>(kPar_momentum,"momentum","p");
  tplist.add<double>(kPar_TotEnergy,"E","totEnergy","TotalEnergy");
  tplist.add<double>(kPar_KineticEnergy,"KE","kineticEnergy");
  tplist.add<double>(kPar_vx,"vx","global_dir_x","global_u");
  tplist.add<double>(kPar_vy,"vy","global_dir_y","global_v");
  tplist.add<double>(kPar_vz,"vz","global_dir_z","global_w");
  tplist.add<double>(kPar_local_vx,"lvx","local_dir_x","local_u");
  tplist.add<double>(kPar_local_vy,"lvy","local_dir_y","local_v");
  tplist.add<double>(kPar_local_vz,"lvz","local_dir_z","local_w");
  tplist.add<double>(kPar_ProcessType,"ProcessType");
  tplist.add<string>(kPar_ProcessName,"Process","ProcessName");
  tplist.add<double>(kPar_Pdg,"pdg","pdgcode","particle");
  tplist.add<double>(kPar_Charge,"charge","ParticleCharge","q");
  tplist.add<double>(kPar_TrackId,"id","trackid");
  tplist.add<double>(kPar_CreatorPdg,"creatorPdg","creator");
  tplist.add<double>(kPar_AncestorPdg,"ParentPdg","AncestorPdg","Ancestor");
  tplist.add<double>(kPar_mass,"mass","m");
  tplist.add<string>(kPar_ParticleName,"ParticleName");
  tplist.add<string>(kPar_CreatorProcessName,"CreatorProcessName","CreatorProcess");



  boost::shared_ptr<RuleFactory> prescaleFactory(new PrescaleRuleFactory);
  boost::shared_ptr<RuleFactory> customFactory(new DetElemContainsRuleFactory<QueriableStepAction>(this));
  tplist.add<int>(kPar_Prescale,"Prescale","by",prescaleFactory);
  tplist.add<const IDetectorElement*>(kPar_DetectorElement,"DetElem","in",customFactory);
  tplist.add<const IDetectorElement*>(kPar_DetectorElement,"DetectorElement","in",customFactory);


  return StatusCode::SUCCESS;
}


StatusCode  QueriableStepAction::GetVertexParameterList(RuleParser::ParameterList& vplist)
{
  GetTrackParameterList(vplist); // Everything valid for a track is valid for a vertex.

  using namespace RuleParser;
  using std::string;


  // And add vertex-specific stuff:
  vplist.add<double>(kPar_Step_dE,"Step_dE","dE");
  vplist.add<double>(kPar_Step_dE_Ion,"Step_dE_Ion","de_ion","ionization");
  vplist.add<double>(kPar_Step_qdE,"Step_qDE","quenched_dE","qdE");
  vplist.add<double>(kPar_Step_dx,"Step_dx","StepLength","dx");
  vplist.add<double>(kPar_Step_dt,"Step_dt","StepDuration","dt");
  vplist.add<double>(kPar_Step_dAngle,"Step_dAngle","dAngle");
  vplist.add<double>(kPar_capTargetZ,"capTargetZ");
  vplist.add<double>(kPar_capTargetA,"capTargetA");
  vplist.add<double>(kPar_Step_E_weighted_x,"Ex","E_weighted_x");
  vplist.add<double>(kPar_Step_E_weighted_y,"Ey","E_weighted_y"); 
  vplist.add<double>(kPar_Step_E_weighted_z,"Ez","E_weighted_z");    
  vplist.add<double>(kPar_Step_E_weighted_t,"Et","E_weighted_t");    
  vplist.add<double>(kPar_Step_qE_weighted_x,"qEx","qE_weighted_x","quenched_weighted_x");
  vplist.add<double>(kPar_Step_qE_weighted_y,"qEy","qE_weighted_y","quenched_weighted_y");    
  vplist.add<double>(kPar_Step_qE_weighted_z,"qEz","qE_weighted_z","quenched_weighted_z");    
  vplist.add<double>(kPar_Step_qE_weighted_t,"qEt","qE_weighted_t","quenched_weighted_t");                                     


  vplist.add<double>(kPar_IsStopping,"IsStopping","stop","End");
  vplist.add<double>(kPar_IsStarting,"IsStarting","start","begin");
  vplist.add<double>(kPar_StepNumber,"StepNumber");                    
  vplist.add<double>(kPar_VolumeChanged,"VolumeChanged","NewVolume");
  vplist.add<double>(kPar_MaterialChanged,"MaterialChanged","NewMaterial");


  verbose() << " kPar_MaterialName " << kPar_MaterialName << " kPar_MaterialChanged " << kPar_MaterialChanged
          << " kPar_VolumeChanged " << kPar_VolumeChanged << " kPar_Pdg " << kPar_Pdg 
          << " kPar_LogicalVolumeName " << kPar_LogicalVolumeName 
          << endreq;
  
  return StatusCode::SUCCESS;
}                                                   
      


// Callback functions for the RuleParser mechanism.
void  QueriableStepAction::queryParam(int id,double& output) const
{

  switch(id) {
    case kPar_IsStopping: 
      output = (  ( mCurrentTrack->GetTrackStatus() == fStopAndKill)
               || ( mCurrentTrack->GetTrackStatus() == fKillTrackAndSecondaries ) );
      break;

    case kPar_IsStarting: 
      output = (mCurrentTrack->GetCurrentStepNumber() == 1);
      break;

    case kPar_StepNumber:
      output = mCurrentTrack->GetCurrentStepNumber();
      break;
 
    case kPar_MaterialChanged:
      output = (mCurrentStep->GetPostStepPoint()->GetMaterial() != mCurrentStep->GetPreStepPoint()->GetMaterial());
      break;

    case kPar_VolumeChanged:
      output = (mCurrentStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
      break;
    
    case kPar_DetectorElementMatch:
      output = getDetectorElementMatch();
      break;

    case kPar_ProcessType:
      output = getProcess().type();
      break;

    case kPar_Pdg:
      output = mCurrentTrack->GetDefinition()->GetPDGEncoding();
      break;

    case kPar_Charge:
      output = mCurrentTrack->GetDefinition()->GetPDGCharge();
      break;
  
    case kPar_TrackId:
      output = mCurrentStep->GetTrack()->GetTrackID();
      break;
      
    // wangzhe
    // Help me track things.
    case kPar_CreatorPdg: {
      G4HistoryUserTrackInfo* userinfo = dynamic_cast<G4HistoryUserTrackInfo*>
        ( mCurrentStep->GetTrack()->GetUserInformation() );
      if(userinfo) 
        output = userinfo->parentPdg();
      break;
      }
      
    case kPar_AncestorPdg: {
      G4HistoryUserTrackInfo* userinfo = dynamic_cast<G4HistoryUserTrackInfo*>
        ( mCurrentStep->GetTrack()->GetUserInformation() );
      if(userinfo)
        output = userinfo->track().track()->particle(); // This will bring me to ancestor track
      break;
    }
    // wz
      
    case kPar_Step_dE:
      output = mCurrentStep->GetTotalEnergyDeposit();
      break;
    case kPar_Step_dE_Ion:
      output = mCurrentStep->GetTotalEnergyDeposit() - mCurrentStep->GetNonIonizingEnergyDeposit();
      break;
    case kPar_Step_qdE:
      output = getQuenchedEnergy();
      break;
    case kPar_Step_dx:
      output = mCurrentStep->GetStepLength();
      break;
    case kPar_Step_dt:
      output = mCurrentStep->GetDeltaTime();
      break;
    case kPar_Step_dAngle:
      output = acos(mCurrentStep->GetPostStepPoint()->GetMomentumDirection().dot(
                     mCurrentStep->GetPreStepPoint()->GetMomentumDirection())
                     ) * CLHEP::radian;
      break;
    case kPar_capTargetZ:
      // Return the neutron capture target Z
      output = m_capinfo->getCapture().GetCapTargetZ();
      //info()<<"stepAction: Z "<<m_capinfo->getCapture().GetCapTargetZ()
       //        <<" A: "<<m_capinfo->getCapture().GetCapTargetA()
        //       <<" time: "<<m_capinfo->getCapture().GetCapTime()
         //      <<" gammaNum: "<<m_capinfo->getCapture().GetCapGammaN()
          //     <<endreq;
      break;
    case kPar_capTargetA:
      // Return the neutron capture target A
      output = m_capinfo->getCapture().GetCapTargetA();
      //info()<<"stepAction: A"<<m_capinfo->getCapture().GetCapTargetA()
       //        <<endreq;               
      break;
      // For energy-weighted quantities, use the average position or time of the step rather than
      // the final position or time of the step.
    case kPar_Step_E_weighted_x:
      output = mCurrentStep->GetPreStepPoint()->GetPosition().x();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetPosition().x());
      output = mCurrentStep->GetTotalEnergyDeposit() * output;
      break;
    case kPar_Step_E_weighted_y:
      output = mCurrentStep->GetPreStepPoint()->GetPosition().y();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetPosition().y());
      output = mCurrentStep->GetTotalEnergyDeposit() * output;
      break;
    case kPar_Step_E_weighted_z:
      output = mCurrentStep->GetPreStepPoint()->GetPosition().z();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetPosition().z());
      output = mCurrentStep->GetTotalEnergyDeposit() * output;
      break;
    case kPar_Step_E_weighted_t:
      output = mCurrentStep->GetTotalEnergyDeposit() * mCurrentStepPoint->GetGlobalTime();
      output = mCurrentStep->GetPreStepPoint()->GetGlobalTime();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetGlobalTime());
      output = mCurrentStep->GetTotalEnergyDeposit() * output;
      break;

    case kPar_Step_qE_weighted_x:
      output = mCurrentStep->GetPreStepPoint()->GetPosition().x();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetPosition().x());
      output = getQuenchedEnergy() * output ; 
      break;
    case kPar_Step_qE_weighted_y:
      output = mCurrentStep->GetPreStepPoint()->GetPosition().y();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetPosition().y());
      output = getQuenchedEnergy() * output ; 
      break;    
    case kPar_Step_qE_weighted_z:
      output = mCurrentStep->GetPreStepPoint()->GetPosition().z();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetPosition().z());
      output = getQuenchedEnergy() * output ; 
      break;
    case kPar_Step_qE_weighted_t:
      output = mCurrentStep->GetPreStepPoint()->GetGlobalTime();
      if ( mCurrentStep->GetPostStepPoint() ) output = 0.5*(output + mCurrentStep->GetPostStepPoint()->GetGlobalTime());
      output = getQuenchedEnergy() * output ; 
      break;

    case kPar_t:
      output = mCurrentStepPoint->GetGlobalTime();
      break;
    case kPar_x:
      output = mCurrentStepPoint->GetPosition().x();
      break;    
    case kPar_y:
      output = mCurrentStepPoint->GetPosition().y();
      break;
    case kPar_z:
      output = mCurrentStepPoint->GetPosition().z();
      break;
    case kPar_r:
      output = sqrt(mCurrentStepPoint->GetPosition().x()*mCurrentStepPoint->GetPosition().x() +
                    mCurrentStepPoint->GetPosition().y()*mCurrentStepPoint->GetPosition().y());
      break;

    case kPar_local_x:
    case kPar_local_y:
    case kPar_local_z:
    case kPar_local_r:
    {
      const G4AffineTransform transformation =
        mCurrentStepPoint->GetTouchable()->GetHistory()->GetTopTransform();
      G4ThreeVector local = transformation.TransformPoint(mCurrentStepPoint->GetPosition());
      
      switch(id){
        case kPar_local_x: output = local.x(); break;
        case kPar_local_y: output = local.y(); break;
        case kPar_local_z: output = local.z(); break;
        case kPar_local_r: output = sqrt(local.x()*local.x() + local.y()*local.y()); break;        
      }
    }
    break;
    
    case kPar_NicheId:
      output = getDetectorId();
      break;

    case kPar_DetectorId:
      output = getDetectorId() & 0xFFFF0000;
      break;
    case kPar_SiteId:
      output = getDetectorId() & 0xFF000000;
      break;
    case kPar_AD:
      output = (getDetectorId() & 0x00FF0000) >> 16;
      break;
    case kPar_Site:
      output = (getDetectorId() & 0xFF000000) >> 24;
      break;

    case kPar_momentum:
       output = mCurrentStepPoint->GetMomentum().mag();
       break;
    case kPar_TotEnergy:
      output = mCurrentStepPoint->GetTotalEnergy();
      break;
    case kPar_KineticEnergy:
      output = mCurrentStepPoint->GetKineticEnergy();
      break;
    case kPar_mass:
      output = mCurrentStep->GetTrack()->GetDynamicParticle()->GetMass();
      break;
    case kPar_vx:
      output = mCurrentStepPoint->GetMomentumDirection().x();
      break;
    case kPar_vy:
      output = mCurrentStepPoint->GetMomentumDirection().y();
      break;
    case kPar_vz:
      output = mCurrentStepPoint->GetMomentumDirection().z();
      break;

    case kPar_local_vx:
    case kPar_local_vy:
    case kPar_local_vz:
    {
      const G4AffineTransform transformationV =
        mCurrentStepPoint->GetTouchable()->GetHistory()->GetTopTransform();
      G4ThreeVector localV = transformationV.TransformAxis(mCurrentStepPoint->GetMomentumDirection());
      
      switch(id){
        case kPar_local_vx: output = localV.x(); break;
        case kPar_local_vy: output = localV.y(); break;
        case kPar_local_vz: output = localV.z(); break;
      }
    }
    break;
        
    default:
      err() << "Unknown real parameter " << id << endreq;
      output = 0;
  }
  verbose() << " id " << id << " output " << output << endreq;

}

void  QueriableStepAction::queryParam(int id,std::string& output) const
{
  const ILVolume* volume = 0;
  const IDetectorElement* ide = 0;

  switch(id) {
  case kPar_CreatorProcessName:
    if(mCurrentTrack->GetParentID()==0) // This is a primary.
        mProcess = SimProcess(SimProcess::kPrimaryVertex,"");
    else if(mCurrentTrack->GetCreatorProcess())
        mProcess = SimProcess(SimProcess::kParticleStart,mCurrentTrack->GetCreatorProcess()->GetProcessName());
    else
        mProcess = SimProcess(SimProcess::kParticleStart,""); 
    output = mProcess.name();
    verbose()<<" id " << id << " creator process name "<<output<<endreq;
    return;
    
  case kPar_ProcessName:
    switch(mCurrentStep->GetPostStepPoint()->GetStepStatus()) {
        case fWorldBoundary:   mProcess = SimProcess(SimProcess::kWorldBoundary,""); break;
        case fGeomBoundary:    mProcess = SimProcess(SimProcess::kGeomBoundary,""); break;
        default:
        SimProcessFromG4Process(mCurrentStep->GetPostStepPoint()->GetProcessDefinedStep(),mProcess);
    }
    output = mProcess.name();
    verbose()<<" id " << id << " process name "<<output<<endreq;
    return;

  case kPar_DetectorElementName:
    ide = getDetectorElement();
    if(ide)  output = ide->name();
    else output = "";
    verbose() << " id " << id << " detector element output " << output << endreq;
    return;

  case kPar_LogicalVolumeName:
    volume = getLogicalVolume();
    if(volume)  output = volume->name();
    else output = "";
    verbose() << " id " << id << " logical volume name output " << output << endreq;
    return;
    
  case kPar_MaterialName:
    volume = getLogicalVolume();
    if(volume)output = volume->materialName();
    else output ="";
    verbose() << " id " << id << " material name output " << output << endreq;
    return;
    
  case kPar_ParticleName:
    output = mCurrentStep->GetTrack()->GetDefinition()->GetParticleName();
    verbose() << " id " << id << " particle name output " << output << endreq;
    return;
    
  default:
    err() << "Unknown string parameter " << id << endreq;
    output = "";
  }
}


// Functions that get useful things and cache them
// Note that current step point should be the pre step point for geometry, material, volume and detector element info

const IDetectorElement* QueriableStepAction::getDetectorElement() const
{
  if(!mDetElement) {
    // wangzhe:
    // (GetPhysicalVolume() == 0) at the boundary of the world.
    // (MotherLogical pointer == 0) means that current step point is
    // in the outmost volume, "Universe". 
    // Currently no touchable history info is available here, so neglect
    // any query in this volume.
    if(mCurrentStep->GetPreStepPoint()->GetPhysicalVolume() == 0) {
      return 0;
    }
    if(mCurrentStep->GetPreStepPoint()->GetPhysicalVolume()->GetMotherLogical() == 0) {
      return 0;
    }
    // wz
    G4VTouchable* touch = mCurrentStep->GetPreStepPoint()->GetTouchableHandle()();
    G4TouchableHistory* hist = dynamic_cast<G4TouchableHistory*>(touch);
    assert(hist);

    if (!hist->GetHistoryDepth()) {
        err() << "getDetectorElement(): given empty touchable history" << endreq;
        return 0;
    }

    StatusCode sc = mTouchToDetElem->GetBestDetectorElement(hist,
                                                            m_DetectorElementSearchPath,
                                                            mDetElement,
                                                            mDetElementMatch);
    if (sc.isFailure()) {
        err() << "Failed to get best detector element" << endreq;
        return 0;
    }
  }
  return mDetElement;  
}

int QueriableStepAction::getDetectorElementMatch() const
{
  getDetectorElement();
  return mDetElementMatch;
}

const ILVolume* QueriableStepAction::getLogicalVolume() const
{
  if(!mLogicVolume) {
    const IPVolume* pvol = getPhysicalVolume();
    if(pvol) mLogicVolume = pvol->lvolume();
  }
  return mLogicVolume; 
}

const IPVolume* QueriableStepAction::getPhysicalVolume() const
{
  if(!mPhysVolume) {
     mTouchToDetElem->G4VolumeToDetDesc(mCurrentStep->GetPreStepPoint()->GetPhysicalVolume(),
                                        mPhysVolume
                                        );
   }
   return mPhysVolume;  
}

const SimProcess& QueriableStepAction::getProcess() const
{
  if(!mProcess.isValid()){
    if(mCurrentStepPoint == mCurrentStep->GetPreStepPoint()) {
      if(mCurrentTrack->GetParentID()==0) // This is a primary.
         mProcess = SimProcess(SimProcess::kPrimaryVertex,"");
      else if(mCurrentTrack->GetCreatorProcess())
        mProcess = SimProcess(SimProcess::kParticleStart,mCurrentTrack->GetCreatorProcess()->GetProcessName());
      else
        mProcess = SimProcess(SimProcess::kParticleStart,"");
    } else {
      switch(mCurrentStep->GetPostStepPoint()->GetStepStatus()) {
        case fWorldBoundary:   mProcess = SimProcess(SimProcess::kWorldBoundary,""); break;
        case fGeomBoundary:    mProcess = SimProcess(SimProcess::kGeomBoundary,""); break;
        default:
        SimProcessFromG4Process(mCurrentStepPoint->GetProcessDefinedStep(),mProcess);
      }
    }
  }
  return mProcess;
}

double QueriableStepAction::getQuenchedEnergy() const
{
  if(mQuenchedEnergy>=0) return mQuenchedEnergy;
  G4Track *aTrack = mCurrentStep->GetTrack();
  G4ParticleDefinition* aParticle = aTrack->GetDefinition();
  G4String particleName = aParticle->GetParticleName();
  //info()<<"particle name "<<particleName<<endreq;

  double dE = mCurrentStep->GetTotalEnergyDeposit();
  double dx = mCurrentStep->GetStepLength();

  // Find quenched energy deposit.
  mQuenchedEnergy = 0;
  if(dE > 0) {
    if(aParticle == G4Gamma::Gamma()) // It is a gamma
    {
      G4LossTableManager* manager = G4LossTableManager::Instance();
      dx = manager->GetRange(G4Electron::Electron(), dE, aTrack->GetMaterialCutsCouple());
      //info()<<"dE_dx = "<<dE/dx/(MeV/mm)<<"MeV/mm"<<endreq;
    } 
    G4Material* aMaterial = mCurrentStep->GetPreStepPoint()->GetMaterial();
    G4MaterialPropertiesTable* aMaterialPropertiesTable =
        aMaterial->GetMaterialPropertiesTable();
    if (aMaterialPropertiesTable) {

      // There are some properties. Is there a scintillator property?
      const G4MaterialPropertyVector* Fast_Intensity = 
          aMaterialPropertiesTable->GetProperty("FASTCOMPONENT"); 
      const G4MaterialPropertyVector* Slow_Intensity =
          aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
  
      if (Fast_Intensity || Slow_Intensity ) {
        // It's a scintillator.
        double delta = dE/dx/aMaterial->GetDensity(); 
        //double birk1 = 0.0125*g/cm2/MeV;
        double birk1 = m_BirksConstant1; 
        if(mCurrentTrack->GetDefinition()->GetPDGCharge()>1.1)//for particle charge greater than 1.
           birk1 = 0.57*birk1;
        //double birk2 = (0.0031*g/MeV/cm2)*(0.0031*g/MeV/cm2);
        double birk2 = m_BirksConstant2;
        mQuenchedEnergy= dE /(1+birk1*delta+birk2*delta*delta);
      }
    }
  }
  return mQuenchedEnergy;  
}

const G4HistoryUserTrackInfo* QueriableStepAction::getUserTrackInfo() const
{
  if(!mUserTrackInfo) {
    mUserTrackInfo = dynamic_cast<G4HistoryUserTrackInfo*>(mCurrentTrack->GetUserInformation());
  }
  return mUserTrackInfo;  
}

unsigned int QueriableStepAction::getDetectorId() const
{
  if(mDetectorId==0xFFFFFFF) {
    const IDetectorElement* de = getDetectorElement();
    mDetectorId = getDetectorId(de);
  }
  return mDetectorId;
}

unsigned int QueriableStepAction::getDetectorId(const IDetectorElement* de) const
{
  DetectorIdCache_t::iterator it = mDetectorIdCache.find(de);
  if(it!=mDetectorIdCache.end()) {
    return it->second;
  }

  const ParamValidDataObject* params = de->params();
  
  // Check each type of user param that the user could have provided:
  for(unsigned int ip = 0; ip<m_IdParameterNames.size(); ip++) {
    if (params->exists(m_IdParameterNames[ip])) {
      if (unsigned int id = (unsigned int)(params->param<int>(m_IdParameterNames[ip]))) {
        // Got it!
        mDetectorIdCache[de] = id;  // cache it for next time
        //info() << "Found id for " << getDetectorElement()->name() << " --> id " << std::hex << mDetectorId << std::dec << endreq;
        return id;
      }
    }
  }
  
  // Ok, we didnt' find it.
  // So, move up the support tree recursively.

  // --this does't work because geo->detElem() is private
  // const IGeometryInfo* igeo = de->geometry();
  //   ILVolume::ReplicaPath dummy;
  //   IGeometryInfo* nextGeo = 0;
  //   igeo->location(nextGeo, dummy); 
  //   if(nextGeo==0) return 0x0; // We failed to find any kind of ID.
  // const GeometryInfoPlus* geoplus = dynamic_cast<GeometryInfoPlus*>(nextGeo);
  //   assert(geoplus); // What's up with this?
  //   return getDetectorId(geoplus->detElem());

  // -- this works if the child elements are all set correctly in the XML.
  const IDetectorElement* nextElem = de->parentIDetectorElement();
  if(nextElem==0) {
    const IDetectorElement* ide = getDetectorElement();
    std::string DeName;
    if(ide) DeName = ide->name();
    else DeName = "";
    warning() << "Found no next-element for " << DeName << " --> id " << std::hex << mDetectorId << std::dec << endreq;
    mDetectorIdCache[de] = 0;
    return 0;
  }
  unsigned int id = getDetectorId(nextElem);
  mDetectorIdCache[de] = id;
  return id;
}

---