source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/G4Detector/optics/src/OptSteppingAction.cc @ 114

#include "OptSteppingAction.hh"
#include "OptUserStackingAction.hh"
#include "Photon.hh"

#include <G4SteppingManager.hh>
#include <G4Track.hh>
#include <G4Step.hh>
#include <G4StepPoint.hh>
#include <G4TrackStatus.hh>
#include <G4VPhysicalVolume.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleTypes.hh>
#include <G4String.hh>
#include <G4ProcessManager.hh>
#include <G4ProcessVector.hh>
#include <G4OpBoundaryProcess.hh>
#include <G4TrackStatus.hh>
#include <G4EventManager.hh>
#include <G4Event.hh>
#include <G4MaterialPropertiesTable.hh>
#include <G4Box.hh>
#include <G4IntersectionSolid.hh>
#include <G4Tubs.hh>
#include <TTree.h>
#include <TVector3.h>
#include <TObject.h>
#include <TObjString.h>

#include "SimuApplication.hh"
#include "ConfigFileParser.hh"
#include "EsafConfigurable.hh"
#include "Config.hh"
#include <iostream>
#include "VirtualTelParm.hh"
#include <G4GeometryTolerance.hh>

using namespace std;
using namespace NTraceLens;
const char* oprocess[] ={ "Undefined", "FresnelRefraction", 
                          "FresnelReflection", "TotalInternalReflection",
                          "LambertianReflection", "LobeReflection",
                          "SpikeReflection", "BackScattering",
                          "Absorption", "Detection", "NotAtBoundary",
                          "SameMaterial", "StepTooSmall", "NoRINDEX" };

// Detector subsystem identifier
enum EDetectorSystem {
    keUndefined               = 0,
    kNotInside                = 1,
    kBaffle                   = 2,
    kOptics                   = 3,
    kFocalPlane               = 4,
    kWalls                    = 5,
    kEarth                    = 6,
    kG4FocalSurface           = 7
};

OptSteppingAction::OptSteppingAction(OptUserStackingAction* stack) : fStepCounter(0), fMaxIterations(1000)
{
    //theStatus = Undefined;
    fBoundary = NULL;
    fProcessManager = NULL;
    fProcessVector = NULL;
    fStatus = 0;
    fStacking = stack;
    ConfigFileParser* conf = Config::Get()->GetCF("G4","G4Detector");
    string det=conf->GetStr("G4Detector.DetectorType");
    //double fGeomTolerance=G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
    if (det.find("Tus")!=string::npos){
        fTus=true;
        cover=0; 
    }
    else {
        fTus=false;
        tel_param* telParm = VirtualTelParm::GetOpticalSystemParam();
        r_cut = telParm->surface[1].r_cut;
        r_wall = telParm->r_wall;
        
// virtual cover without baffle

        G4Box* cbox = new G4Box("cbox",(r_wall+1.)*mm,r_wall*mm,3.5*m);//(r_cut+10.)
        G4Tubs* ctub =new G4Tubs("ctub",0.,r_wall*mm,3.*m,0.*deg,360.*deg);
        G4IntersectionSolid* intersec = new G4IntersectionSolid("cbox*ctub", cbox, ctub);
        
        cover = new G4DisplacedSolid ("cover",intersec,0,G4ThreeVector(0,0,2.5*m));
        
    }
    SetSender("OptSteppingAction");

}

OptSteppingAction::~OptSteppingAction()
{

}

void OptSteppingAction::UserSteppingAction(const G4Step * theStep)
{
    G4VPhysicalVolume* prevol=theStep->GetPreStepPoint()->GetPhysicalVolume();
    G4VPhysicalVolume* postvol=theStep->GetPostStepPoint()->GetPhysicalVolume();
    G4String PreStepVolume,PostStepVolume;
    if(prevol){
        PreStepVolume = prevol->GetName();
    }
    if(postvol){
        PostStepVolume = postvol->GetName();
    }
    
    G4StepPoint *PostStep = 0;
    G4StepPoint *PreStep = 0;
    G4ThreeVector pre_pos, pre_dir;
    G4ThreeVector post_pos, post_dir;
    
    PreStep = theStep->GetPreStepPoint();
    pre_pos = PreStep->GetPosition();
    pre_dir = PreStep->GetMomentumDirection();
    PostStep = theStep->GetPostStepPoint();
    post_pos =  PostStep->GetPosition();
    post_dir =  PostStep->GetMomentumDirection();

    
    fStepCounter++;
    if ( fStepCounter>fMaxIterations ){
        theStep->GetTrack()->SetTrackStatus(fStopAndKill);
        MsgForm(EsafMsg::Warning,"This photon exceeded the maximum number of iterations (%i), killed.",fMaxIterations);
    }

    //find the boundary process only once
    if( !fBoundary ){
        fProcessManager = theStep->GetTrack()->GetDefinition()->GetProcessManager();
        G4int nprocesses = fProcessManager->GetProcessListLength();
        fProcessVector = fProcessManager->GetProcessList();

        for(int i=0;i<nprocesses;i++){
            if( (*fProcessVector)[i]->GetProcessName()=="OpBoundary" ){
                fBoundary = (G4OpBoundaryProcess*)(*fProcessVector)[i];
                break;
            }
        }
    }

    theStatus = Undefined;
    if(fBoundary)theStatus = fBoundary->GetStatus();

    double rho = post_pos.perp();
    if(!fTus && (fabs(post_pos.y())>=r_cut || rho>=r_wall)){
        if(cover->Inside (pre_pos)==kInside ){
            theStep->GetTrack()->SetTrackStatus(fStopAndKill);
            fStatus=-4;
            double dist = cover->DistanceToOut (pre_pos, pre_dir);
            // intersection with virtual wall
            
            if(dist!=kInfinity ){
              //printf("dist %g\n",dist);
                pre_pos+=pre_dir.unit()*dist;
                PostStep->SetPosition(pre_pos);
                PostStep->SetMomentumDirection(pre_dir);
                post_pos=pre_pos;
                post_dir=pre_dir;
                double t = PostStep->GetLocalTime();
                t+=dist/c_light;
                PostStep->SetLocalTime(t);
            }
        }
        else if(cover->Inside (pre_pos)==kSurface){
            
            //double newdist= cover->DistanceToIn (pre_pos, pre_dir);
            //if(newdist==kInfinity ){
              PostStep->SetPosition(pre_pos);
              //printf("dist %g\n",newdist);
            //}
//          if(newdist!=kInfinity ){
//              pre_pos-=pre_dir.unit()*newdist;
                //pre_pos.z()=0.;
//              PostStep->SetPosition(pre_pos);
//              PostStep->SetMomentumDirection(pre_dir);
//              post_pos=pre_pos;
//              post_dir=pre_dir;
//              double t = PostStep->GetLocalTime();
//              t-=newdist/c_light;
//              PostStep->SetLocalTime(t);
//              
//          
//          
//          
//          }
//          
//      
        }
        else PostStep->SetPosition(pre_pos);
        
    }

  
    if ( PostStep && fStatus!=-4 && PostStep->GetStepStatus()==fPostStepDoItProc &&
         theStep->GetTrack()->GetTrackStatus()==fStopAndKill ) fStatus=-3;

    if(postvol){
      
        if ( PostStepVolume.contains("FresnelLens") || PreStepVolume.contains("FresnelLens")||
            PreStepVolume.contains("mirror_phys") ||PostStepVolume.contains("mirror_phys") ){
            if(fStatus==0) fStatus++;
        }

        if(PostStepVolume=="absorber"){
            
            theStep->GetTrack()->SetTrackStatus(fStopAndKill);
            fStatus=(fStatus>0?-2:-1);
            if(fStatus==-1||fStatus==-2)PostStep->SetPosition(pre_pos);
        }
    }

   
    if ( theStatus==Absorption || theStatus==Detection ){
        fStatus=-3;
        
    }
    G4TrackStatus TrackStatus = theStep->GetTrack()->GetTrackStatus();

    if(TrackStatus==fStopAndKill){
        fStepCounter=0;
        int history;
        if ( post_pos.z()<0 && post_dir.z()<0  ){
            history=kEarth;
            double dist = cover->DistanceToOut (pre_pos, pre_dir);
        // intersection with virtual wall
            
            if(dist!=kInfinity && dist!=0){

                pre_pos+=pre_dir.unit()*dist;
                PostStep->SetPosition(pre_pos);
                PostStep->SetMomentumDirection(pre_dir);
                post_pos=pre_pos;
                post_dir=pre_dir;
                double t = PostStep->GetLocalTime();
                t+=dist/c_light;
                PostStep->SetLocalTime(t);
            }
         
         }
         else if ( fStatus==-2 ) history=kG4FocalSurface;
         else if ( fStatus==-3) history=kOptics;
         else if ( fStatus==-4) history=kWalls;
         else history=keUndefined;

         fStacking->SetTrackHistory(history);
         fStatus=0;
    }
}