In This Package:

HistorianStepAction.cc

Go to the documentation of this file.

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

#include "Event/SimProcess.h"
#include "Event/SimVertex.h"
#include "Event/SimTrack.h"
#include "Event/SimParticleHistory.h"

#include "G4TouchableHistory.hh"
#include "G4TouchableHandle.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4SteppingManager.hh"
#include "G4MaterialPropertiesTable.hh"

#include "RuleParser/Rules.h"
#include "RuleParser/CreateRules.h"
#include "RuleParser/PrescaleRule.h"

#include "G4DataHelpers/G4DhPrimaryParticleInformation.h"
#include "G4DataHelpers/ITouchableToDetectorElement.h"
#include "G4DataHelpers/IHistoryKeeper.h"

#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/DeclareFactoryEntries.h"

#include <vector>


DECLARE_TOOL_FACTORY(HistorianStepAction);

using namespace DayaBay;

HistorianStepAction::HistorianStepAction ( 
  const std::string& type   , 
  const std::string& name   , 
  const IInterface*  parent ) 
  : QueriableStepAction( type, name, parent )
  , mHistoryKeeper(0)
  , mTrackRule(0)
  , mVertexRule(0)
{ 
  declareProperty("TrackSelection",m_TrackSelection = "none");
  declareProperty("VertexSelection",m_VertexSelection = "none");
  declareProperty("UseFastMuEnergyCut",m_UseFastMuEnergyCut = false);
}

StatusCode HistorianStepAction::initialize() 
{
  info() << "HistorianStepAction::initialize()" << endreq;

  StatusCode sc = QueriableStepAction::initialize();
  if (sc.isFailure()) return sc;
  
  // Set up default search path for detector elements. Note that this throws away PMT-specific information.
  // Also, change the default search path:
  m_DetectorElementSearchPath.clear();
  m_DetectorElementSearchPath.push_back("/dd/Structure/DayaBay");
  m_DetectorElementSearchPath.push_back("/dd/Structure/Sites");
  m_DetectorElementSearchPath.push_back("/dd/Structure/Pool");
  m_DetectorElementSearchPath.push_back("/dd/Structure/AD");
  m_DetectorElementSearchPath.push_back("/dd/Structure/RPC");
  
  sc = service("HistoryKeeper",mHistoryKeeper,true);
  if (sc.isFailure()) return sc;
  

  // Build the selection rules.
  
  using namespace RuleParser;
  using std::string;
    
  ParameterList tplist;
  QueriableStepAction::GetTrackParameterList(tplist);
  tplist.add<double>(kPar_ParentPdg,"ParentPdg","AncestorPdg","Ancestor");
  tplist.add<double>(kPar_ParentIndirection,"ParentIndirection","AncestorIndirection");
  tplist.add<double>(kPar_GrandParentPdg,"GrandParentPdg","GrandParent");
  tplist.add<double>(kPar_GrandParentIndirection,"GrandParentIndirection");


  ParameterList vplist;
  QueriableStepAction::GetVertexParameterList(vplist);
  vplist.add<double>(kPar_ParentPdg,"ParentPdg","AncestorPdg","Ancestor");
  vplist.add<double>(kPar_ParentIndirection,"ParentIndirection","AncestorIndirection");
  vplist.add<double>(kPar_GrandParentPdg,"GrandParentPdg","GrandParent");
  vplist.add<double>(kPar_GrandParentIndirection,"GrandParentIndirection");
  
  vplist.add<double>(kPar_distanceFromLastVertex,"distanceFromLastVertex");
  vplist.add<double>(kPar_timeSinceLastVertex,"TimeSinceLastVertex");
  vplist.add<double>(kPar_energyLossSinceLastVertex,"EnergyLostSinceLastVertex");
  vplist.add<double>(kPar_angleFromLastVertex,"AngleFromLastVertex");
                                                     
                                                     
  bool result;
  result = RuleParser::CreateRules(m_TrackSelection,tplist,mTrackRule);
  
  if(mTrackRule==0 || result==false) {
    err() << "Failed to interpret selection string: " << m_TrackSelection << endreq;
    return StatusCode::FAILURE;
  } else {
    info() << "Configured with track selection:  " << mTrackRule->name() << endreq;
  }
  

  result = RuleParser::CreateRules(m_VertexSelection,vplist,mVertexRule);
  
  if(mVertexRule==0 || result==false) {
    err() << "Failed to interpret selection string: " << m_VertexSelection << endreq;
    return StatusCode::FAILURE;
  } else {
    info() << "Configured with vertex selection: " << mVertexRule->name() << endreq;
  }
  

  return StatusCode::SUCCESS; 
}

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

// Callback functions for the RuleParser mechanism.


void  HistorianStepAction::queryParam(int id,double& output) const
{
  SimTrackReference ref;
  G4StepPoint* postPoint = mCurrentStep->GetPostStepPoint();
  const SimVertex* v = 0;
  Gaudi::XYZPoint pos;
  Gaudi::XYZVector dPos;
  
  switch(id) {
    case kPar_ParentPdg:
      ref = getAncestorTrack();
      if(ref.track()) output = ref.track()->particle();
      else output = 0;
      break;
      
    case kPar_ParentIndirection:
      output = getAncestorTrack().indirection();
      break;

    case kPar_GrandParentPdg:
      output = 0;
      ref = getAncestorTrack();
      if(ref.track()) ref = ref.track()->ancestorTrack();
      if(ref.track()) output = ref.track()->particle();
      break;
    
    case kPar_GrandParentIndirection:
      output = -999;
      ref = getAncestorTrack();
      if(ref.track()) output = ref.track()->ancestorTrack().indirection();
      break;

    // Special history stuff.
    case kPar_distanceFromLastVertex:      
    case kPar_timeSinceLastVertex:
    case kPar_energyLossSinceLastVertex:
    case kPar_angleFromLastVertex:
      assert(mCurrentTrackRef.track());
      assert(mCurrentTrackRef.track()->vertices().size());
      v = mCurrentTrackRef.track()->vertices().back();

      switch(id) {
        case kPar_distanceFromLastVertex:
          pos = Gaudi::XYZPoint(postPoint->GetPosition().x()
                                ,postPoint->GetPosition().y()
                                ,postPoint->GetPosition().z());
          dPos = pos - v->position();
          output = dPos.R();
          return;
        case kPar_timeSinceLastVertex:
          output = postPoint->GetGlobalTime() - v->time();
          return;
        case kPar_energyLossSinceLastVertex:
          output = postPoint->GetTotalEnergy() - v->totalEnergy();
          return;
        case kPar_angleFromLastVertex:
          output = ( postPoint->GetMomentumDirection().x() * v->momentum().x()
                   + postPoint->GetMomentumDirection().y() * v->momentum().y()
                   + postPoint->GetMomentumDirection().z() * v->momentum().z() )
                   / (v->momentum().R());
          output = acos(output)*CLHEP::radian;
          return;
      }
      break;

    default:
      QueriableStepAction::queryParam(id,output);
  }
  verbose() << " id " << id << " output " << output << endreq;

}

void  HistorianStepAction::queryParam(int id,std::string& output) const
{
  QueriableStepAction::queryParam(id,output);
  verbose() << " id " << id << " output " << output << endreq;

}


SimTrackReference HistorianStepAction::getAncestorTrack() const
{
  const G4HistoryUserTrackInfo* info = getUserTrackInfo();
  if(!info) {
    err() << "Found track id " << mCurrentTrack->GetTrackID() << " with no user info. Shouldn't happen." << endreq;
    err() << "(Are you running with both the HistorianTrackAction AND the HistorianStepAction?)" << endreq;
    assert(info);    
  }
  return info->track();
}

SimVertexReference HistorianStepAction::getAncestorVertex() const
{
  const G4HistoryUserTrackInfo* info = getUserTrackInfo();
  if(!info) {
    err() << "Found track id " << mCurrentTrack->GetTrackID() << " with no user info. Shouldn't happen." << endreq;
    err() << "(Are you running with both the HistorianTrackAction AND the HistorianStepAction?)" << endreq;
    assert(info);    
  }
  
  return info->vertex();
}


// Function to fill out the vertex.

StatusCode HistorianStepAction::CreateVertex(SimVertex* &vertex,
                                           const SimTrackReference& parent)
{ 
     
  const G4ThreeVector& pos = mCurrentStepPoint->GetPosition();
  const G4ThreeVector mom = mCurrentStepPoint->GetMomentum();

  vertex = new SimVertex(
                      parent,
                      getProcess(),
                      mCurrentStepPoint->GetGlobalTime(),
                      Gaudi::XYZPoint(pos.x(),pos.y(),pos.z()),
                      mCurrentStepPoint->GetTotalEnergy(),
                      Gaudi::XYZVector(mom.x(),mom.y(),mom.z())
                     );

  return StatusCode::SUCCESS;
}

// Function to fill out the track.


StatusCode HistorianStepAction::CreateTrack(const G4Track* g4track,SimTrack* &track)
{
  SimTrackReference parentTrack;
  SimVertexReference parentVertex;
  int parentPdg = 0;
 
  const HepMC::GenParticle*  genPart = 0;

  if(g4track->GetParentID()==0){
    //info() << "Building primary particle" << endreq;
    // Primary track is a special case.
    parentTrack = SimTrackReference(0,0);
    parentVertex= SimVertexReference(0,0);

    const G4DhPrimaryParticleInformation* primaryinfo = 
        dynamic_cast<const G4DhPrimaryParticleInformation*>(
             g4track->GetDynamicParticle()->GetPrimaryParticle()->GetUserInformation() );
     if(primaryinfo) {
       genPart  = primaryinfo->GetHepParticle();      
     }

  } else {
    const G4HistoryUserTrackInfo* userinfo = getUserTrackInfo();
    if(userinfo) {
      parentTrack = userinfo->track();
      parentVertex = userinfo->vertex();
      parentPdg = userinfo->parentPdg();
    } else {
      err() << "Found track id " << mCurrentTrack->GetTrackID() << " with no user info. Shouldn't happen." << endreq;
      err() << "(Are you running with both the HistorianTrackAction AND the HistorianStepAction?)" << endreq;
      assert(userinfo);    
    }
    genPart = parentTrack.track()->primaryParticle();
  }

  int pdg = pdgFromG4Track(g4track);
  if(pdg == 0 ) {
    warning() << "Found particle with zero pdg code. Geant claims this particle is: " 
              << g4track->GetDefinition()->GetParticleName() << endreq;
  }
 
  track = new SimTrack();
  track->setTrackId(g4track->GetTrackID());
  track->setParentParticle(parentPdg);
  track->setAncestorTrack(parentTrack);
  track->setAncestorVertex(parentVertex);
  track->setPrimaryParticle(genPart);
  track->setParticle(pdg);
    
  // Create the start vertex.
  const G4ThreeVector& pos = mCurrentStepPoint->GetPosition();
  const G4ThreeVector mom = mCurrentStepPoint->GetMomentum();

  // Figure out the process.
  SimProcess proc;
  if(g4track->GetParentID()==0){
    // This is a primary.
    proc = SimProcess(SimProcess::kPrimaryVertex,"");
  } else {
    if(mCurrentTrack->GetCreatorProcess())
      proc = SimProcess(SimProcess::kParticleStart,mCurrentTrack->GetCreatorProcess()->GetProcessName());
    else
      proc = SimProcess(SimProcess::kParticleStart,"");
  }

  SimVertex* vertex = new SimVertex(
                       SimTrackReference(track,0),
                       proc,
                       mCurrentStepPoint->GetGlobalTime(),
                       Gaudi::XYZPoint(pos.x(),pos.y(),pos.z()),
                       mCurrentStepPoint->GetTotalEnergy(),
                       Gaudi::XYZVector(mom.x(),mom.y(),mom.z())
                      );
    
  track->addVertex(vertex);
  mCurrentHistory->addTrack(track);
  mCurrentHistory->addTrackRef(g4track->GetTrackID(),SimTrackReference(track,0));
  mCurrentHistory->addVertex(vertex);
  if(g4track->GetParentID()==0){ 
    // This is a primary.
    mCurrentHistory->addPrimaryTrack(track);
  } else {
    // This track is not a primary, so try to link from our parent vertex to us.
    // Invert the parent-vertex indirection to find out the parentVertex->this track reference.
    SimTrackReference downref(track,parentVertex.indirection());
    if(!parentVertex.vertex())
      err() << "SimTrack being created, but no Parent vertex is stored. This shouldn't happen." << endreq;
    else
      parentVertex.vertex()->addSecondary(downref);
  }  
  return StatusCode::SUCCESS;
}


bool HistorianStepAction::IsInterestingVertex(const G4Step*)
{
  // All the magic is done by the RuleParser system, and the Rules it has procuded:
  return mVertexRule->select(this);
}

bool HistorianStepAction::IsInterestingTrack(const G4Track*)
{
  // All the magic is done by the RuleParser system, and the Rules it has procuded:
  return mTrackRule->select(this);
}


void HistorianStepAction::UserSteppingAction(const G4Step* g4step)
{
  
  // Reset temporary things:
  mHistoryKeeper->GetCurrentHistory(mCurrentHistory);
  assert(mCurrentHistory);
 
  //Add the fast muon simulation energy cut for different particle, intend to increase the simulation speed. By Haoqi Jan,27, 2011.
  
  G4Track* g4track = g4step->GetTrack();
  if( g4track->GetGlobalTime()> 3e16 )  {  //3e16 /* a year */
    G4HistoryUserTrackInfo* userinfo = dynamic_cast<G4HistoryUserTrackInfo*>
      ( g4track->GetUserInformation() );
    int parentPdg = 0;
    if(userinfo) {
      parentPdg = userinfo->parentPdg();
    }
    warning()<<"A 1 year old track orignated by particle "<<parentPdg<<endreq;
  }

  G4int trackId = g4track ->GetTrackID();
  G4StepPoint* thePrePoint;
  G4VPhysicalVolume* physV;
  G4Material* material;
  G4ParticleDefinition* particle = g4track->GetDefinition();
  G4String pname = particle->GetParticleName();
  G4ProcessVector* processVector; 
  thePrePoint = g4step->GetPreStepPoint();
  physV = thePrePoint->GetPhysicalVolume();
  material = physV->GetLogicalVolume()->GetMaterial(); 
  G4double EnergyCut;
 
if(g4track->GetCurrentStepNumber()==1&&m_UseFastMuEnergyCut==true)
  { idflag=true;}
 //G4cout<<"FastMuCut ="<<m_UseFastMuEnergyCut <<G4endl;
  if(idflag==true)
    {
       
      if(material->GetName()=="/dd/Materials/Rock"||material->GetName()=="/dd/Materials/Air")
      {  
       //  G4cout<<"Rock or air aterial Name ="<<material->GetName()<<G4endl;
         if(particle->GetPDGCharge()!=0)
         {
            EnergyCut=500.0;
         }
         else
         {
            EnergyCut=30.0; 
         }
      }//rock
      else if(material->GetName()=="/dd/Materials/OwsWater"||material->GetName()=="/dd/Materials/IwsWater")
      {
         // G4cout<<"Ows or IWs aterial Name ="<<material->GetName()<<G4endl;
         if(particle->GetPDGCharge()!=0)
         {
            EnergyCut=50.0;
         }
         else
         {
            EnergyCut=3.0;   
          }
       }
       else 
       {
        EnergyCut=0.;
     //    if(material->GetName()=="/dd/Materials/LiquidScintillator"||material->GetName()=="/dd/Materials/GdDopedLS") G4cout<<"Material ="<<material->GetName()<<G4endl;

        }//material
      if(g4track->GetKineticEnergy()/MeV<EnergyCut)
      {
          
          //aTrack->SetKineticEnergy(0.00001/MeV);
          g4track->SetTrackStatus(fStopAndKill);
          idflag=false;
      }
   }//idflag
  
  // First, do things related to SimTracks.

  if(g4track->GetCurrentStepNumber() == 1) {
    // This is the first time we've seen this track. 
    // The interesting point is the first one, the track initial vertex.
    QueriableStepAction::Reset(g4step,g4track,g4step->GetPreStepPoint());

    // Does a SimTrack already exist for this G4track?
    int id = g4track->GetTrackID();
    SimTrackReference ref = mCurrentHistory->track(id);
    SimTrack* track = ref.track();

    if(track==0) { // There is no track or reference.
      // Maybe we need to make one. Is it a primary track?
      if(g4track->GetParentID()==0) {
        // This is a primary track, so this action is responsible for starting it.
        // Create the SimTrack
        CreateTrack(g4track,track);
        ref = SimTrackReference(track,0);
      } else {
        if(IsInterestingTrack(g4track)) {
          CreateTrack(g4track,track);        
          ref = SimTrackReference(track,0);
        }
      }
    }

    if(ref.track()==0) {
      // OK, there is still no valid reference. Let's build an indirect one.
      SimTrackReference  parentTrack = getAncestorTrack();
      SimVertexReference parentVertex = getAncestorVertex();

       // Increment indirection.
       ref = SimTrackReference(parentTrack.track(),parentTrack.indirection()+1);
       mCurrentHistory->addTrackRef(g4track->GetTrackID(),ref);
     }


    // currentTrack should now be set correctly, one way or another
    if(! ref.isDirect()) {
      // The reference is indirect, meaning we aren't recording the current one.
      // So, instead, we want to store some 'missing' information.
      ref.track()->incrementUnrecordedDescendants(pdgFromG4Track(g4track));
    }

    // The current track reference, for convenience.
    //mCurrentHistory->setCurrentTrack(ref);
  }
  
  // Second, do things related to SimVertex.
  
  // The point of interest is the post-step if we're looking at the step/vertex
  QueriableStepAction::Reset(g4step,g4track,g4step->GetPostStepPoint());
  
  mCurrentTrackRef = mCurrentHistory->track(g4track->GetTrackID());
  if(!mCurrentTrackRef.track()) {
    // This is an error and shouldn't happen.
    err() << "Somehow we got a bad track reference on track id " 
          << g4track->GetTrackID() << endreq;
    err() << " This shouldn't happen - everything should have a parent." << endreq;
    assert(0);
  }
  

  // Is the current track getting recorded?
  if(mCurrentTrackRef.isDirect()){
    // Yes. Is this interaction interesting on it's own?
    if(IsInterestingVertex(g4step)) {
      verbose() << " InterestingVertex. Track in " << g4step->GetTrack()->GetVolume()->GetName() 
              << " particle " << g4step->GetTrack()->GetDefinition()->GetParticleName()
              << " G4TrackID "<< g4step->GetTrack()->GetTrackID()
              << " current step# " << g4step->GetTrack()->GetCurrentStepNumber()
              << " step length(mm) " << g4step->GetTrack()->GetStepLength()/mm
              << endreq;
      SimVertex* vertex = 0;
      CreateVertex(vertex,mCurrentTrackRef);
      // Give ownership to History.
      mCurrentHistory->addVertex(vertex); 
      // We have a recorded track and recorded vertex. So, add the vertex.
      mCurrentTrackRef.track()->addVertex(vertex);
    }
  }

  // Get current PDG code.
  int pdg = pdgFromG4Track(g4track);
  
  // Get the current best reference to a vertex
  SimVertexReference curVertexRef(mCurrentTrackRef.track()->vertices().back(),mCurrentTrackRef.indirection());

  assert(fpSteppingManager);
  assert(fpSteppingManager->GetStep() == g4step);
  
  // Loop over secondaries, push information forward.
  int nSecAtRest = fpSteppingManager->GetfN2ndariesAtRestDoIt();
  int nSecAlong  = fpSteppingManager->GetfN2ndariesAlongStepDoIt();
  int nSecPost   = fpSteppingManager->GetfN2ndariesPostStepDoIt();
  int nSecTotal  = nSecAtRest+nSecAlong+nSecPost;
  G4TrackVector* secVec = fpSteppingManager->GetfSecondary();

  size_t end = (*secVec).size();
  size_t start = end-nSecTotal;
  for(size_t itrack=start; itrack < end; itrack++)
  {
    G4Track* g4track = (*secVec)[itrack];
    G4HistoryUserTrackInfo* user
      = new G4HistoryUserTrackInfo(curVertexRef,mCurrentTrackRef,pdg);
    if(g4track->GetUserInformation()) {
      warning() << "Already found user info!  itrack=" << itrack
        << " out of " << (*secVec).size() << endreq;
    } else {
      g4track->SetUserInformation(user);
    }
  }

}

---