In This Package:

UnObserverStepAction.cc

Go to the documentation of this file.

#include "UnObserverStepAction.h"

#include "Event/SimUnobservableStatisticsHeader.h"

#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4Navigator.hh" // need this for debug (world->local position)

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

#include "G4DataHelpers/IHistoryKeeper.h"

#include "GaudiKernel/DeclareFactoryEntries.h"

DECLARE_TOOL_FACTORY(UnObserverStepAction);

using namespace DayaBay;


struct AddableVector : public std::vector<std::string>
{
  AddableVector& add(const std::string& s) { push_back(s); return *this; }
};

UnObserverStepAction::UnObserverStepAction ( 
  const std::string& type   , 
  const std::string& name   , 
  const IInterface*  parent ) 
  : QueriableStepAction( type , name , parent )
  , m_Stats()
{ 
   // Declare some good starting stats.
  m_Stats.push_back(AddableVector().add("photon_created_energy" ).add( "E" ).add( "StepNumber==1 and pdg==20022" ) );
  m_Stats.push_back(AddableVector().add("photon_backscatter_r" ).add( "local_r" ).add( "pdg==20022 and dAngle >= 90" ) );
  m_Stats.push_back(AddableVector().add("photon_forescatter_r" ).add( "local_r" ).add( "pdg==20022 and (dAngle <= 90 and dAngle > 0)" ) );
  m_Stats.push_back(AddableVector().add("photon_stop_r" ).add(        "local_r" ).add( "(pdg==20022) and ((IsStopping == 1) and (LogicalVolumeName == '/dd/Geometry/PMT/lvPmtHemi'))" ) );
  m_Stats.push_back(AddableVector().add("pmt_hit" ).add( "t" ).add( "pdg==20022 and (IsStopping == 1 and (LogicalVolumeName=='/dd/Geometry/PMT/lvPmtHemi'))" ) );
  m_Stats.push_back(AddableVector().add("edep-water").add(  "dE").add( "pdg!=20022 and (MaterialName == '/dd/Materials/Water')") );
  m_Stats.push_back(AddableVector().add("edep-oil").add(    "dE").add( "pdg!=20022 and (MaterialName == '/dd/Materials/Oil')") );
  m_Stats.push_back(AddableVector().add("edep-ls").add(     "dE").add( "pdg!=20022 and (MaterialName == '/dd/Materials/LiquidScintillator')") );
  m_Stats.push_back(AddableVector().add("edep-gdls").add(   "dE").add( "pdg!=20022 and (MaterialName == '/dd/Materials/GdDopedLS')") );
  m_Stats.push_back(AddableVector().add("edep-acryilc").add("dE").add( "pdg!=20022 and (MaterialName == '/dd/Materials/Acrylic')") );
  m_Stats.push_back(AddableVector().add("edep-ad1"  ).add("dE"  ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==1)"  ) );
  m_Stats.push_back(AddableVector().add("edep-ad2"  ).add("dE"  ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==2)"  ) );
  m_Stats.push_back(AddableVector().add("edep-ad3"  ).add("dE"  ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==3)"  ) );
  m_Stats.push_back(AddableVector().add("edep-ad4"  ).add("dE"  ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==4)"  ) );
  m_Stats.push_back(AddableVector().add("qedep-ad1" ).add("qdE" ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==1)"   ) );
  m_Stats.push_back(AddableVector().add("qedep-ad2" ).add("qdE" ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==2)"   ) );
  m_Stats.push_back(AddableVector().add("qedep-ad3" ).add("qdE" ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==3)"   ) );
  m_Stats.push_back(AddableVector().add("qedep-ad4" ).add("qdE" ).add("pdg!=20022 and ((MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS') and AD==4)"   ) );
  m_Stats.push_back(AddableVector().add("mu-path-water").add( "dx").add( "(pdg==13 || pdg==-13) and MaterialName == '/dd/Materials/Water'") );
  m_Stats.push_back(AddableVector().add("mu-path-oil").add(   "dx").add( "(pdg==13 || pdg==-13) and MaterialName == '/dd/Materials/Oil'") );
  m_Stats.push_back(AddableVector().add("mu-path-ls").add(    "dx").add( "(pdg==13 || pdg==-13) and MaterialName == '/dd/Materials/LiquidScintillator'") );
  m_Stats.push_back(AddableVector().add("mu-path-gdls").add(  "dx").add( "(pdg==13 || pdg==-13) and MaterialName == '/dd/Materials/GdDopedLS'") );
  m_Stats.push_back(AddableVector()
    .add("ad"    ).add( "AD" )
    .add("expos" ).add("Ex"  )
    .add("eypos" ).add("Ey"  )
    .add("ezpos" ).add("Ez"  )
    .add("et"    ).add("Et"  )
    .add("qexpos").add("qEx" )
    .add("qeypos").add("qEy" )
    .add("qezpos").add("qEz" )
    .add("qet"   ).add("qEt" )
    .add( "(pdg!=20022) and (MaterialName == '/dd/Materials/LiquidScintillator' or MaterialName == '/dd/Materials/GdDopedLS')"  ) );
   
   declareProperty("Stats",m_Stats);
   
   
   // Also, change the default search path:
   m_DetectorElementSearchPath =  AddableVector()
     .add("/dd/Structure/DayaBay")
     .add("/dd/Structure/Sites")
     .add("/dd/Structure/Pool")
     .add("/dd/Structure/AD")
     .add("/dd/Structure/RPC");
   
   
}

UnObserverStepAction::~UnObserverStepAction ()
{
  for(size_t i=0;i<mGroups.size();i++) {
    if(mGroups[i].mRule) delete mGroups[i].mRule;
  }   
}


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

  StatusCode sc = QueriableStepAction::initialize();
  if (sc.isFailure()) return sc;
  
  sc = service("HistoryKeeper",mKeeper,true);
  if (sc.isFailure()) return sc;

  // Build the selection rules.
  
  using namespace RuleParser;
  using std::string;

  ParameterList plist;
  QueriableStepAction::GetVertexParameterList(plist);

  // Now loop through the various statistics that have been offered.
  for(size_t istat = 0; istat < m_Stats.size(); istat++ ){
    std::vector<std::string>& stat = m_Stats[istat];

    // Check basic layout. Odd number of strings.
    if(stat.size()%2 != 1) {
      err() << "Bad configuration given to the UnObserverStepAction:" << endreq;
      err() << " Expected strings of the form <name> : <parameter> [: <name> : <parameter> ..]: <cut>" << endreq;
      err() << " but got " << stat.size() << " strings"
      //<< " \"" << line << << "\""
       <<  endreq;
      return StatusCode::FAILURE;
    }

    StatGroup statgroup;
    
    // Compile the cut
    statgroup.mRule = 0;
    bool compiled = RuleParser::CreateRules(stat[stat.size()-1],plist,statgroup.mRule);
    if(!compiled || !statgroup.mRule) {
      err() << "Bad configuration given to the UnObserverStepAction:" << endreq;
      err() << " Expected 3 strings of the form <name>,<parameter>,<cut>" << endreq;
      err() << " but I failed to interpret selection string: " << stat[stat.size()-1] << endreq;
      return StatusCode::FAILURE;      
    }
    info() << "Cut group : " << stat[stat.size()-1] << endreq;
    
    for(size_t isub=0; isub<(stat.size()-1)/2 ; isub+=2) {
      // Check that the paramter name matches something, and record the match.
      std::string parname = stat[isub+1];
      for(std::string::iterator c = parname.begin(); c!= parname.end(); ++c)  *c = tolower(*c);
    
      bool foundit = false;
      int parid = 0;
      for(size_t ipar = 0; ipar < plist.size(); ipar++) {
        std::string temp = plist[ipar].name();
        for(std::string::iterator c = temp.begin(); c!= temp.end(); ++c)  *c = tolower(*c);
        if(parname==temp) {
          foundit = true;
          parid = plist[ipar].id();
        }
      }
      if(!foundit) {   
        err() << "Bad configuration given to the UnObserverStepAction:" << endreq; 
        err() << "Didn't recognize parameter name " << parname << endreq;
        return StatusCode::FAILURE;
      }
      // Report success and record the info we'll need at runtime.
      info() << "   Statistic " << stat[isub] 
             << " param=" << parname << " (id=" << parid <<")"
             << endreq;
       
      NamedStat ns;
      ns.mName = stat[isub];
      ns.mParam = parid;

      statgroup.push_back(ns);
    }
    if(statgroup.size())
     mGroups.push_back(statgroup);
  }

  // keep track of bail-out messages
  mBail = 0; 

  return StatusCode::SUCCESS; 
}

StatusCode UnObserverStepAction::finalize() 
{
  info() << "UnObserverStepAction::finalize()" << endreq;
  info() << "There were "<<mBail<<" bail-out warnings encountered in UserSteppingAction."<<endreq;
  return QueriableStepAction::finalize(); 
}



void UnObserverStepAction::UserSteppingAction(const G4Step* g4step)
{

  // Use the PreStepPoint to correctly evaluate geometry-, material-, volume-related information
  // Use the PostStepPoint to evaluate 'track'-related infomation (ie, energy loss, step length)

  // Reference material:
  // The algorithm to handle a single step in G4:
  // http://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch05.html#sect.Track.Basic
  // G4 stepping info from http://geant4.cern.ch/support/faq.shtml#TRACK-1 :
  // "Hereafter we call current volume the volume where the step has just 
  //   gone through. Geometrical informations are available from preStepPoint."
  // 
  // Reset temporary things:
  QueriableStepAction::Reset(g4step,g4step->GetTrack(),g4step->GetPostStepPoint());
  if ( g4step->GetTrack()->GetCurrentStepNumber() == 1 )
    // This is the first step.  
    // The interesting point is the first one.
    QueriableStepAction::Reset(g4step,g4step->GetTrack(),g4step->GetPreStepPoint());

  // The following is useful for intensive debugging of each step on every track.
  // Make it so it does not get executed without re-compilation.
  if ( 0 && g4step->GetTrack()->GetDefinition()->GetParticleName() != "opticalphoton" )
    {
      G4StepPoint* pre = g4step->GetPreStepPoint();
      G4String mpre = "NONE", mpost = "NONE", status1 = "", status2 = "";
      G4String ppre = "NONE",ppost = "NONE",lpre = "NONE",lpost = "NONE";
      G4ThreeVector localPosition = 0, worldPosition = 0;
      if (pre) {
        if ( pre->GetStepStatus() == fGeomBoundary ) status1 = "Entering";
        G4Material* mm = pre->GetMaterial();
        if (mm) mpre = mm->GetName();
        G4TouchableHandle touch = pre->GetTouchableHandle();
        G4VPhysicalVolume* volume = touch->GetVolume();
        ppre = volume->GetName();
        G4LogicalVolume* lvolume = volume->GetLogicalVolume();
        lpre = lvolume->GetName();
      }
      else {
        verbose() << " THERE IS NO PRE-STEP-POINT FOR THIS STEP " << endreq;
      }
      pre = g4step->GetPostStepPoint();
      if (pre) {
        if (pre->GetStepStatus() == fGeomBoundary ) status2 = "Leaving";
        G4Material* mm = pre->GetMaterial();
        if (mm) mpost = mm->GetName();
        G4TouchableHandle touch = pre->GetTouchableHandle();
        if (touch) {
          worldPosition = pre->GetPosition();
          localPosition = touch->GetHistory()->GetTopTransform().TransformPoint( worldPosition );
          G4VPhysicalVolume* volume = touch->GetVolume();
          if (volume) {
            ppost = volume->GetName();
            G4LogicalVolume* lvolume = volume->GetLogicalVolume();
            if (lvolume)          lpost = lvolume->GetName();
          }
        }
      }
      else {
        verbose() << " THERE IS NO POST-STEP-POINT FOR THIS STEP " << endreq;
      }

      verbose() << "Track in " << g4step->GetTrack()->GetVolume()->GetName() 
              << ", particle is " << g4step->GetTrack()->GetDefinition()->GetParticleName()
              << " Etot(MeV) " << g4step->GetTrack()->GetTotalEnergy()/MeV
              << " Ekin(MeV) " << g4step->GetTrack()->GetKineticEnergy()/MeV
              << " worldX,Y,Z " << worldPosition.x() << ", " << worldPosition.y() << ", " << worldPosition.z()
              << " localX,Y,Z " << localPosition.x() << ", " << localPosition.y() << ", " << localPosition.z()
              << " DirX,Y,Z " << g4step->GetTrack()->GetMomentumDirection().x()
              << ", " << g4step->GetTrack()->GetMomentumDirection().y()
              << ", " << g4step->GetTrack()->GetMomentumDirection().z()
              << " G4TrackID "<< g4step->GetTrack()->GetTrackID()
              << " current step# " << g4step->GetTrack()->GetCurrentStepNumber()
              << " step length(mm) " << g4step->GetTrack()->GetStepLength()/mm
              << " status " << g4step->GetTrack()->GetTrackStatus() << endreq;
      verbose() << " ...." << status1 << " " << status2 
              << " prestep (current) material " << mpre << " pVol " << ppre << " lVol " << lpre 
              << " poststep material " << mpost << " pVol " << ppost << " lVol " << lpost
              << endreq;
      }
        
  // djaffe Bail with warning if we are in uppermost 'universe' volume
  if ( g4step->GetTrack()->GetVolume()->GetMotherLogical() == 0x0 ) {
    mBail++;
    if (mBail <= 10) { 
      warning() << "Bailing out because track is in " << g4step->GetTrack()->GetVolume()->GetName() 
                << ", particle is " << g4step->GetTrack()->GetDefinition()->GetParticleName()
                << " with total energy(MeV) " << g4step->GetTrack()->GetTotalEnergy()/MeV
                << " & kinetic energy(MeV) " << g4step->GetTrack()->GetKineticEnergy()/MeV
                << " G4TrackID "<< g4step->GetTrack()->GetTrackID()
                << " current step# " << g4step->GetTrack()->GetCurrentStepNumber()
                << " step length(mm) " << g4step->GetTrack()->GetStepLength()/mm
                << " status " << g4step->GetTrack()->GetTrackStatus()
                << endreq;
    }
    if (mBail > 10) {
      debug()   << "Bailing out because track is in " << g4step->GetTrack()->GetVolume()->GetName() 
                << ", particle is " << g4step->GetTrack()->GetDefinition()->GetParticleName()
                << " with total energy(MeV) " << g4step->GetTrack()->GetTotalEnergy()/MeV
                << " & kinetic energy(MeV) " << g4step->GetTrack()->GetKineticEnergy()/MeV
                << " G4TrackID "<< g4step->GetTrack()->GetTrackID()
                << " current step# " << g4step->GetTrack()->GetCurrentStepNumber()
                << " step length(mm) " << g4step->GetTrack()->GetStepLength()/mm
                << " status " << g4step->GetTrack()->GetTrackStatus()
                << endreq;
    }
    if (mBail == 10) { 
      warning() << "Additional bail-out messages will be suppressed" << endreq;
    }
    return;
  }
  
  // Get the current set of stats.
  SimUnobservableStatisticsHeader* ush = 0;
  StatusCode sc = mKeeper->GetCurrentUnobservable(ush);
  if(sc.isFailure() || (!ush)) {
    err() << "Problem with Keeper." << endreq;
    return;
  }
  assert(ush);
  
  SimUnobservableStatisticsHeader::stat_map& statmap = ush->stats();
  if(ush->stats().size() == 0) {

    // This is a first call.  Make sure that all the stats get intialized, by simply looking them up.
    for(size_t igrp = 0; igrp < mGroups.size(); igrp++) {
      StatGroup& statgroup = mGroups[igrp];
      for(size_t istat = 0; istat < statgroup.size(); istat++ ){
        statmap[statgroup[istat].mName];
        verbose() << "Initializing group " << igrp << " stat " << istat 
               << " name=" <<statgroup[istat].mName << " paramid=" << statgroup[istat].mParam << endreq;
      }
    }  
  }

  for(size_t igrp = 0; igrp < mGroups.size(); igrp++) {
    StatGroup& statgroup = mGroups[igrp];
    
    // The magic line.. does this step pass the cut?
    verbose() << " Group " << igrp << " statgroup.mRule->name() " << statgroup.mRule->name() << endreq;
    if(statgroup.mRule->select(this)) {
      // yes, so loop over the parameters.
      for(size_t istat = 0; istat < statgroup.size(); istat++ ){
        double x = 0;
        queryParam(statgroup[istat].mParam,x);
        SimStatistic& stat = statmap[statgroup[istat].mName];
        verbose() << "Querying group " << igrp << " stat " << istat 
                << " name=" <<statgroup[istat].mName << " paramid=" << statgroup[istat].mParam 
                << " and incrementing by " << x << endreq;
        stat.increment(x);
      }
    }
  }

}

---