MdEdepSimulatorAG/SteppingAction.cc

Go to the documentation of this file.

#include "SteppingAction.h"

#include <fwk/CentralConfig.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/MathConstants.h>

//G4 classes
#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 <string>
using std::string;

using namespace fwk;
using namespace utl;
using namespace EdepSimulatorAG;

SteppingAction::SteppingAction () : fScintID(-1){}

SteppingAction::~SteppingAction (  )
{;}

void SteppingAction::UserSteppingAction ( const G4Step * theStep )
{
  G4Track* theTrack = theStep->GetTrack (  );
  G4double edep     = theStep->GetTotalEnergyDeposit();
 
  if( !edep )
    return;

  G4StepPoint       * thePostPoint = theStep->GetPostStepPoint (  );
  //G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume (  ); // Frank Schröder: commented since unused
  G4ParticleDefinition * particleType = theTrack->GetDefinition (  );

  G4String process = "Undefined";
  if ( thePostPoint->GetProcessDefinedStep() != NULL)
   process = thePostPoint->GetProcessDefinedStep()->GetProcessName();

  G4double x =  theStep->GetPostStepPoint()->GetPosition().x();
  G4double y =  theStep->GetPostStepPoint()->GetPosition().y();
  G4double z =  theStep->GetPostStepPoint()->GetPosition().z();
  G4VPhysicalVolume *thePrePV = theStep->GetPreStepPoint (  )->GetPhysicalVolume (  );
  string  prePVName  = string (thePrePV->GetName());
  string  postPVName = "OutOfWorld";

  G4int    trackid  = theTrack->GetTrackID();
  G4int    parentid = theTrack->GetParentID();
 
  int sId = fScintID; // fScintID is [1-64] by construction 

  if ( prePVName     == "RightScint" ) 
    sId += 1;
  else if ( prePVName == "LeftScint" ) 
    sId -= 1;

  //do not add virtual strips
  if( sId == 65 || sId == 0 )// strips 0 and 65 are always virtual
    return;
  else if( fScintID == 32 && sId == 33 ) // 33 is virtual for strips in the range [1-32]
    return;
  else if( fScintID == 33 && sId == 32 ) // 32 is virtual for strips in the range [33-64]
    return;

  //If scintillator Id (key) not found, then creates the key
  if ( fScintMap.find( sId ) == fScintMap.end()){
    DataCollection  nuovo;
    //        G4cerr << "Scinti ID not found... I will creat it" << G4endl;
    fScintMap[sId] = nuovo;
    //TODO: keep the identity of the "primary" impinging particle.
    //should have partenID == 0 for central scintillator but not so for right and left scintillators
  }   

  Data stepd;
  //Fill obj with relevant step data
  stepd.Sname   = prePVName;
  stepd.Pname   = particleType->GetParticleName();
  stepd.PDGCode = particleType->GetPDGEncoding();
  stepd.ParentID= parentid;
  stepd.TrackID = trackid;
  stepd.x       = x;
  stepd.y       = y;
  stepd.z       = z;
  stepd.px      = theStep->GetPreStepPoint()->GetMomentumDirection().x();
  stepd.py      = theStep->GetPreStepPoint()->GetMomentumDirection().y();
  stepd.pz      = theStep->GetPreStepPoint()->GetMomentumDirection().z();
  stepd.KinE    = theStep->GetPreStepPoint()->GetKineticEnergy();
  stepd.Edep    = edep;
  stepd.Time    = theStep->GetPreStepPoint()->GetGlobalTime();

  //Fill the vector with the data of a given scintillator
  fScintMap[sId].push_back( stepd );
/*
  G4cerr << "\n===========================================\n";
  G4cerr << "        PARTICLE CROSSING A SCINTILLATOR STRIP!!!!            \n";
  G4cerr <<  particleType->GetParticleName() << " ";
  G4cerr <<  "PDG code: " << particleType->GetPDGEncoding()  << " ";
  G4cerr << "prePV: " << prePVName << " ";
  G4cerr << "trackID " << trackid << G4endl;
  G4cerr << "Edep " << edep/CLHEP::MeV << G4endl;
  G4cerr << "X:" << x/CLHEP::m << " m ";
  G4cerr << "Y:" << y/CLHEP::m << " m ";
  G4cerr << "Z:" << z/CLHEP::m << " m " << G4endl;
  G4cerr << "===========================================\n";
*/                                 
  return;

}