diff --git a/examples/ILDExSimu/CMakeLists.txt b/examples/ILDExSimu/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..faad5ea7613223b99750763488fdac95ba615b6e --- /dev/null +++ b/examples/ILDExSimu/CMakeLists.txt @@ -0,0 +1,26 @@ +cmake_minimum_required(VERSION 2.8.3 FATAL_ERROR) + +find_package(Geant4 REQUIRED ui_all vis_all) +##---Handle the case CLHEP is not included in Geant4------------------------------ +#if(NOT Geant4_clhep_FOUND) +# find_package(CLHEP REQUIRED) +# set(Geant4_INCLUDE_DIRS ${Geant4_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}) +# set(Geant4_LIBRARIES ${Geant4_LIBRARIES} ${CLHEP_LIBRARIES}) +#endif() +INCLUDE(${Geant4_USE_FILE}) # this also takes care of geant 4 definitions and include dirs + + +find_package(LCIO REQUIRED) + + +include_directories( ${CMAKE_SOURCE_DIR}/DDCore/include + ${CMAKE_SOURCE_DIR}/DDG4/include + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${ROOT_INCLUDE_DIR} + ${Geant4_INCLUDE_DIRS} + ${VGM_INCLUDE_DIR} + ${LCIO_INCLUDE_DIRS}) + +file(GLOB sources src/*.cpp) +add_executable(ILDExSimu ILDExSimu.cpp ${sources}) +target_link_libraries(ILDExSimu DD4hepCore DD4hepG4 ILDEx ${Geant4_LIBRARIES} ${LCIO_LIBRARIES}) diff --git a/examples/ILDExSimu/ILDExSimu.cpp b/examples/ILDExSimu/ILDExSimu.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b781d7334063d00bcf3483383b3c703907e0bc77 --- /dev/null +++ b/examples/ILDExSimu/ILDExSimu.cpp @@ -0,0 +1,170 @@ + +// This example is adapted from the Geant4 example N03 + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4RunManager.hh" +#include "G4UImanager.hh" +#include "G4UIsession.hh" +#include "Randomize.hh" +#include "G4VisExecutive.hh" +#include "G4UIterminal.hh" +#include "G4UIExecutive.hh" +#include "G4UItcsh.hh" +#ifdef G4VIS_USE_QT +#include "G4UIQt.hh" +#endif +#include "QGSP_BERT.hh" + +//#include "ILDExDetectorConstruction.hh" +#include "ILDExPhysicsList.h" +#include "ILDExPrimaryGeneratorAction.h" +#include "ILDExRunAction.h" +#include "ILDExEventAction.h" +#include "ILDExSteppingAction.h" +#include "ILDExSteppingVerbose.h" + +#include "DDG4/Geant4DetectorConstruction.h" +#include "DD4hep/LCDD.h" + +// -- lcio -- +// #include <UTIL/BitField64.h> +// #include <UTIL/ILDConf.h> +#include "lcio.h" +#include "IO/LCWriter.h" + + +using namespace DD4hep::Geometry; +using namespace DD4hep::Simulation; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +int main(int argc,char** argv) +{ + + std::string lcioOutFile("ILDExSimu.slcio") ; + + // -- open LCIO file ---- + lcio::LCWriter* lcWrt = lcio::LCFactory::getInstance()->createLCWriter() ; + + + // Choose the Random engine + // + CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine); + + LCDD& lcdd = LCDD::getInstance(); + + if( argc < 3 ){ + std::cout << " --- Usage: \n " + << " [1] ./bin/ILDExSimu file:../DDExamples/ILDExDet/compact/ILDEx.xml file:../DDExamples/ILDExDet/compact/geant4.xml run1_g4.mac \n" + << " [2] ./bin/ILDExSimu -i file:../DDExamples/ILDExDet/compact/ILDEx.xml file:../DDExamples/ILDExDet/compact/geant4.xml \n" + << " [1]: batch mode - [2]: interactive " << std::endl ; + exit( 0 ) ; + } + + + +//************************************************************** + + bool isBatchMode = ( std::string( argv[1] ) != "-i" ) ; + + int argStart = ( isBatchMode ? 1 : 2 ) ; + int argEnd = ( isBatchMode ? argc-1 : argc ) ; + + for(int i=argStart; i < argEnd ; ++i ) { + + // We need to construct the geometry at this level already + lcdd.fromCompact(argv[i]); + } + + //************************************************************** + + // Get the detector constructed + // + G4VUserDetectorConstruction* detector = new Geant4DetectorConstruction( lcdd ) ; + + // User Verbose output class + // + G4VSteppingVerbose::SetInstance(new ILDExSteppingVerbose); + + // Construct the default run manager + // + G4RunManager * runManager = new G4RunManager; + + + runManager->SetUserInitialization(detector); + + // + G4VUserPhysicsList* physics = new QGSP_BERT ; //new ILDExPhysicsList; + runManager->SetUserInitialization(physics); + + + // Set user action classes + // + G4VUserPrimaryGeneratorAction* gen_action = + new ILDExPrimaryGeneratorAction(lcdd); + runManager->SetUserAction(gen_action); + + //--- + ILDExRunAction* run_action = new ILDExRunAction; + run_action->runData.lcioWriter = lcWrt ; + run_action->runData.detectorName = lcdd.header().name() ; + runManager->SetUserAction(run_action); + + // + ILDExEventAction* event_action = new ILDExEventAction(run_action, lcdd ); + runManager->SetUserAction(event_action); + // + G4UserSteppingAction* stepping_action = + new ILDExSteppingAction(event_action); + runManager->SetUserAction(stepping_action); + + // Initialize G4 kernel + // + runManager->Initialize(); + + // Initialize visualization + // + G4VisManager* visManager = new G4VisExecutive; + visManager->Initialize(); + + // Get the pointer to the User Interface manager + + lcWrt->open( lcioOutFile , lcio::LCIO::WRITE_NEW ) ; + + G4UImanager* UImanager = G4UImanager::GetUIpointer(); + + if ( isBatchMode) { // batch mode + + G4String command = "/control/execute "; + G4String fileName = argv[argc-1]; + UImanager->ApplyCommand(command+fileName); + + } else { // interactive mode : define UI session + + // G4UIsession *ui = new G4UIterminal(new G4UItcsh()); + // G4UIsession* ui = new G4UIQt(argc, argv); + G4UIExecutive* ui = new G4UIExecutive(argc, argv); + + ui->SessionStart(); + // end ... + delete ui; + } + + // Job termination + // Free the store: user actions, physics_list and detector_description are + // owned and deleted by the run manager, so they should not + // be deleted in the main() program ! + + + delete visManager; + delete runManager; + + lcWrt->close() ; + delete lcWrt ; + + return 0; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/examples/ILDExSimu/include/ILDExEventAction.h b/examples/ILDExSimu/include/ILDExEventAction.h new file mode 100644 index 0000000000000000000000000000000000000000..dcc5be7d7ed24f157c97e1323e16c3f6a6f89bd3 --- /dev/null +++ b/examples/ILDExSimu/include/ILDExEventAction.h @@ -0,0 +1,51 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifndef ILDExEventAction_h +#define ILDExEventAction_h 1 + +#include "G4UserEventAction.hh" +#include "globals.hh" + +class ILDExRunAction; +class ILDExEventActionMessenger; + +namespace DD4hep{namespace Geometry{ class LCDD ; }} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExEventAction : public G4UserEventAction +{ +public: + ILDExEventAction(ILDExRunAction*, DD4hep::Geometry::LCDD& ); + + virtual ~ILDExEventAction(); + + void BeginOfEventAction(const G4Event*); + void EndOfEventAction(const G4Event*); + + void SumSupport(G4double de, G4double dl, G4double da) {EnergySupport += de; TrackLSupport += dl; AngleSupport += da; }; + void SumSensitive(G4double de, G4double dl, G4double da) {EnergySensitive += de; TrackLSensitive += dl; AngleSensitive += da; }; + + void SetPrintModulo(G4int val) {printModulo = val;}; + +private: + ILDExRunAction* runAct; + + G4double EnergySupport, EnergySensitive; + G4double TrackLSupport, TrackLSensitive; + G4double AngleSupport, AngleSensitive; + + G4int printModulo; + + ILDExEventActionMessenger* eventMessenger; + + DD4hep::Geometry::LCDD& _lcdd ; +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif + + diff --git a/examples/ILDExSimu/include/ILDExEventActionMessenger.h b/examples/ILDExSimu/include/ILDExEventActionMessenger.h new file mode 100644 index 0000000000000000000000000000000000000000..97191427ee980a2bb2f56599e053e277918e083a --- /dev/null +++ b/examples/ILDExSimu/include/ILDExEventActionMessenger.h @@ -0,0 +1,33 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifndef ILDExEventActionMessenger_h +#define ILDExEventActionMessenger_h 1 + +#include "globals.hh" +#include "G4UImessenger.hh" + +class ILDExEventAction; +class G4UIdirectory; +class G4UIcmdWithAnInteger; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExEventActionMessenger: public G4UImessenger +{ +public: + ILDExEventActionMessenger(ILDExEventAction*); + virtual ~ILDExEventActionMessenger(); + + void SetNewValue(G4UIcommand*, G4String); + +private: + ILDExEventAction* eventAction; + G4UIdirectory* eventDir; + G4UIcmdWithAnInteger* PrintCmd; +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif diff --git a/examples/ILDExSimu/include/ILDExPhysicsList.h b/examples/ILDExSimu/include/ILDExPhysicsList.h new file mode 100644 index 0000000000000000000000000000000000000000..0a658ba446d101a0ed28f53b62b37e2ac2bf2cd5 --- /dev/null +++ b/examples/ILDExSimu/include/ILDExPhysicsList.h @@ -0,0 +1,37 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifndef ILDExPhysicsList_h +#define ILDExPhysicsList_h 1 + +#include "G4VUserPhysicsList.hh" +#include "globals.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExPhysicsList: public G4VUserPhysicsList +{ +public: + ILDExPhysicsList(); + virtual ~ILDExPhysicsList(); + + // Construct particle and physics + void ConstructParticle(); + void ConstructProcess(); + + void SetCuts(); + +private: + + // these methods Construct physics processes and register them + void ConstructDecay(); + void ConstructEM(); +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif + + + diff --git a/examples/ILDExSimu/include/ILDExPrimaryGeneratorAction.h b/examples/ILDExSimu/include/ILDExPrimaryGeneratorAction.h new file mode 100644 index 0000000000000000000000000000000000000000..4d45243c659159485bfecc3564c61cbbcbf01b7f --- /dev/null +++ b/examples/ILDExSimu/include/ILDExPrimaryGeneratorAction.h @@ -0,0 +1,39 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifndef ILDExPrimaryGeneratorAction_h +#define ILDExPrimaryGeneratorAction_h 1 + +#include "G4VUserPrimaryGeneratorAction.hh" +#include "DD4hep/LCDD.h" +#include "globals.hh" + +class G4ParticleGun; +class G4Event; +class ILDExDetectorConstruction; +class ILDExPrimaryGeneratorMessenger; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction +{ +public: + ILDExPrimaryGeneratorAction(const DD4hep::Geometry::LCDD&); + virtual ~ILDExPrimaryGeneratorAction(); + + void GeneratePrimaries(G4Event*); + void SetRndmFlag(G4String val) { rndmFlag = val;} + +private: + G4ParticleGun* particleGun; //pointer a to G4 class + const DD4hep::Geometry::LCDD& ILDExDetector; //pointer to the geometry + ILDExPrimaryGeneratorMessenger* gunMessenger; //messenger of this class + G4String rndmFlag; //flag for a rndm impact point +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif + + diff --git a/examples/ILDExSimu/include/ILDExRunAction.h b/examples/ILDExSimu/include/ILDExRunAction.h new file mode 100644 index 0000000000000000000000000000000000000000..dd4e14d82933b82092fe0dbad197ad3dd4fe4fe6 --- /dev/null +++ b/examples/ILDExSimu/include/ILDExRunAction.h @@ -0,0 +1,59 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifndef ILDExRunAction_h +#define ILDExRunAction_h 1 + +#include "G4UserRunAction.hh" +#include "globals.hh" + +#include "IO/LCWriter.h" +#include <string> + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class G4Run; + +// helper structfor 'global' run data: +struct RunData{ + std::string detectorName ; + IO::LCWriter* lcioWriter ; + RunData() : detectorName("UNKNOWN") , lcioWriter(0) {} +} ; + + +class ILDExRunAction : public G4UserRunAction +{ +public: + ILDExRunAction(); + virtual ~ILDExRunAction(); + + void BeginOfRunAction(const G4Run*); + void EndOfRunAction(const G4Run*); + + void fillPerEvent(G4double ESupport, G4double ESensitive, + G4double LSupport, G4double LSensitive, + G4double AngleSupport, G4double AngleSensitive); + + const G4Run* g4run ; + RunData runData ; + +private: + + G4double sumESupport, sum2ESupport; + G4double sumESensitive, sum2ESensitive; + + G4double sumLSupport, sum2LSupport; + G4double sumLSensitive, sum2LSensitive; + + G4double sumAngleSupport, sum2AngleSupport; + G4double sumAngleSensitive, sum2AngleSensitive; + + +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif + diff --git a/examples/ILDExSimu/include/ILDExSteppingAction.h b/examples/ILDExSimu/include/ILDExSteppingAction.h new file mode 100644 index 0000000000000000000000000000000000000000..7a341e596cfd5bdcd4b5b22f6ff49d8b7c855d89 --- /dev/null +++ b/examples/ILDExSimu/include/ILDExSteppingAction.h @@ -0,0 +1,30 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifndef ILDExSteppingAction_h +#define ILDExSteppingAction_h 1 + +#include "G4UserSteppingAction.hh" + +class ILDExDetectorConstruction; +class ILDExEventAction; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExSteppingAction : public G4UserSteppingAction +{ +public: + ILDExSteppingAction(ILDExEventAction*); + virtual ~ILDExSteppingAction(); + + void UserSteppingAction(const G4Step*); + +private: + ILDExDetectorConstruction* detector; + ILDExEventAction* eventaction; +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif diff --git a/examples/ILDExSimu/include/ILDExSteppingVerbose.h b/examples/ILDExSimu/include/ILDExSteppingVerbose.h new file mode 100644 index 0000000000000000000000000000000000000000..cf7b1a4ee5602ac156c49a9b3fb260a53848fa67 --- /dev/null +++ b/examples/ILDExSimu/include/ILDExSteppingVerbose.h @@ -0,0 +1,28 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExSteppingVerbose; + +#ifndef ILDExSteppingVerbose_h +#define ILDExSteppingVerbose_h 1 + +#include "G4SteppingVerbose.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +class ILDExSteppingVerbose : public G4SteppingVerbose +{ + public: + + ILDExSteppingVerbose(); + ~ILDExSteppingVerbose(); + + void StepInfo(); + void TrackingStarted(); + +}; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif diff --git a/examples/ILDExSimu/run1.mac b/examples/ILDExSimu/run1.mac new file mode 100644 index 0000000000000000000000000000000000000000..1b84c5acbbd90e350885a6b9711eb7d028bdb50c --- /dev/null +++ b/examples/ILDExSimu/run1.mac @@ -0,0 +1,22 @@ +# $Id: run1.mac,v 1.2 2000-11-21 10:59:42 maire Exp $ +# +# Macro file for "exampleN03.cc" +# +# can be run in batch, without graphic +# or interactively: Idle> /control/execute run1.mac +# +/control/verbose 1 +/control/saveHistory +# +/run/verbose 2 +/event/verbose 0 +/tracking/verbose 1 +# +# muon 300 MeV to the direction (1.,0.,0.) +# 3 events +# +##/generator/select gun +/gun/direction 1. 1. 1. +/gun/particle pi+ +/gun/energy 100 MeV +/run/beamOn 3 diff --git a/examples/ILDExSimu/src/ILDExEventAction.cpp b/examples/ILDExSimu/src/ILDExEventAction.cpp new file mode 100644 index 0000000000000000000000000000000000000000..478b78d79f5db63fad4432509f5abdc189533e90 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExEventAction.cpp @@ -0,0 +1,289 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExEventAction.h" +#include "ILDExRunAction.h" +#include "ILDExEventActionMessenger.h" + +#include "G4Run.hh" +#include "G4Event.hh" +#include "G4TrajectoryContainer.hh" +#include "G4VTrajectory.hh" +#include "G4VVisManager.hh" +#include "G4UnitsTable.hh" + +#include "Randomize.hh" +#include <iomanip> + +#include "DDG4/Geant4Hits.h" +#include "DD4hep/VolumeManager.h" +#include "DD4hep/Volumes.h" + +//--- lcio +#include "lcio.h" +#include "IMPL/LCEventImpl.h" +#include "IMPL/LCCollectionVec.h" +#include "IMPL/SimTrackerHitImpl.h" +#include "IMPL/SimCalorimeterHitImpl.h" +#include "UTIL/Operators.h" +#include "UTIL/ILDConf.h" + + +#define DEBUG 0 + + +//------ helper functions ------------------ + +lcio::SimTrackerHitImpl* createSimTrackerHit( DD4hep::Simulation::Geant4TrackerHit* gh ){ + + lcio::SimTrackerHitImpl* lh = new lcio::SimTrackerHitImpl ; + + lh->setCellID0( ( gh->cellID >> 0 ) & 0xFFFFFFFF ) ; + lh->setCellID1( ( gh->cellID >> sizeof( int ) ) & 0xFFFFFFFF ) ; + + const double pos[3] = { gh->position.x() , gh->position.y() , gh->position.z() } ; + lh->setPosition( pos ) ; + + lh->setEDep( gh->energyDeposit ) ; + + lh->setTime( gh->truth.time ) ; + + lh->setMomentum( gh->momentum.x(), gh->momentum.y() , gh->momentum.z() ) ; + + lh->setPathLength( gh->length ) ; + + return lh ; +} + +//-------------- +lcio::SimCalorimeterHitImpl* createSimCalorimeterHit( DD4hep::Simulation::Geant4CalorimeterHit* gh ){ + + lcio::SimCalorimeterHitImpl* lh = new lcio::SimCalorimeterHitImpl ; + + lh->setCellID0( ( gh->cellID >> 0 ) & 0xFFFFFFFF ) ; + lh->setCellID1( ( gh->cellID >> sizeof( int ) ) & 0xFFFFFFFF ) ; + + const float pos[3] = { gh->position.x() , gh->position.y() , gh->position.z() } ; + lh->setPosition( pos ) ; + + lh->setEnergy( gh->energyDeposit ) ; + + return lh ; +} + +//------------------------------------------ + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExEventAction::ILDExEventAction(ILDExRunAction* run, DD4hep::Geometry::LCDD& lcdd ) + :runAct(run),printModulo(1),eventMessenger(0), _lcdd( lcdd ) +{ + eventMessenger = new ILDExEventActionMessenger(this); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExEventAction::~ILDExEventAction() +{ + delete eventMessenger; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExEventAction::BeginOfEventAction(const G4Event* evt) +{ + G4int evtNb = evt->GetEventID(); + if (evtNb%printModulo == 0) { + G4cout << "\n---> Begin of event: " << evtNb << G4endl; + CLHEP::HepRandom::showEngineStatus(); +} + + // initialisation per event + EnergySupport = EnergySensitive = 0.; + TrackLSupport = TrackLSensitive = 0.; + AngleSupport = AngleSensitive = 0.; + + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExEventAction::EndOfEventAction(const G4Event* evt) +{ + + + //accumulates statistic + // + runAct->fillPerEvent(EnergySupport, EnergySensitive, TrackLSupport, TrackLSensitive, AngleSupport, AngleSensitive); + + //print per event (modulo n) + // + G4int evtNb = evt->GetEventID(); + + + if (evtNb%printModulo == 0) { + +#if DEBUG + + G4cout << "---> End of event: " << evtNb << G4endl; + + G4cout + << " Support: total energy: " << std::setw(7) + << G4BestUnit(EnergySupport,"Energy") + << " total track length: " << std::setw(7) + << G4BestUnit(TrackLSupport,"Length") + << G4endl + << " Sensitive: total energy: " << std::setw(7) + << G4BestUnit(EnergySensitive,"Energy") + << " total track length: " << std::setw(7) + << G4BestUnit(TrackLSensitive,"Length") + << G4endl; + +#endif + } + + G4HCofThisEvent* hce = evt->GetHCofThisEvent() ; + + int nCol = hce->GetNumberOfCollections() ; + + lcio::LCEventImpl* lcEvt = new lcio::LCEventImpl ; + + lcEvt->setEventNumber( evt->GetEventID() ) ; + // lcEvt->setRunNumber( evt->GetRunID() ) ; + +#if DEBUG + G4cout << " ILDExEventAction::EndOfEventAction: HCE has " << nCol << " collections : " << G4endl ; +#endif + + for(int i=0 ; i<nCol ; ++i ){ + + G4VHitsCollection* hCol = hce->GetHC( i) ; + +#if DEBUG + G4cout << " --- " << i << ": " << hCol->GetName() << " from SD: " << hCol->GetSDname() << " size: " << hCol->GetSize() ; +#endif + + + bool isTracker = ( hCol->GetSize() ? dynamic_cast<DD4hep::Simulation::Geant4TrackerHit* >( hCol->GetHit(0) ) : 0 ) ; + bool isCalorimeter = ( hCol->GetSize() ? dynamic_cast<DD4hep::Simulation::Geant4CalorimeterHit*>( hCol->GetHit(0) ) : 0 ) ; + + if( hCol->GetSize() > 0 ){ + +#if DEBUG + if ( isTracker ) G4cout << " - type Geant4TrackerHit " << G4endl; + else if( isCalorimeter ) G4cout << " - type Geant4CalorimeterHit " << G4endl ; + else G4cout << " - type UNKNOWN " << G4endl; +#endif + + + + + + // ------ get the readout/cellId description string from first element: + + + std::string cellIDDesc("") ; + + DD4hep::Geometry::VolumeManager vm = _lcdd.volumeManager(); + + DD4hep::Geometry::VolumeManager::VolumeID volume_id = dynamic_cast<DD4hep::Simulation::Geant4Hit*>( hCol->GetHit(0) )->cellID ; + + std::cout << " looking up placed volume for id " << std::hex << volume_id << std::dec << std::endl ; + + if( volume_id ) { + const DD4hep::Geometry::PlacedVolume & pv = vm.lookupPlacement ( volume_id ) ; + + +#if DEBUG + const DD4hep::Geometry::DetElement & detElem = vm.lookupDetElement( volume_id) ; + if ( detElem.isValid() ) + std::cout << " ILDExEventAction::EndOfEventAction --- for detector element : " << detElem.name() << std::endl ; + else + std::cout << " ILDExEventAction::EndOfEventAction --- detector element not found " << std::endl ; +#endif + + + if( pv.isValid() && pv.volume().isSensitive() ) { + + DD4hep::Geometry::Volume vol = pv.volume(); + DD4hep::Geometry::SensitiveDetector sd = vol.sensitiveDetector(); + DD4hep::Geometry::Readout ro = sd.readout(); + DD4hep::Geometry::IDDescriptor iddesc = ro.idSpec(); + + + cellIDDesc = iddesc.fieldDescription() ; + + } else { + + std::cout << " **** WARNING: could not get sensitive placedVolume for cellID : " << std::hex << volume_id << std::dec << std::endl ; + } + + } + + if( isTracker ) { //----------------------------------------------------------------- + + lcio::LCCollectionVec* col = new lcio::LCCollectionVec( lcio::LCIO::SIMTRACKERHIT ) ; + + // the encoder sets the correct cellid encoding string + // ILDCellIDEncoder<SimTrackerHit> idDec( col ) ; + + UTIL::CellIDEncoder<SimTrackerHit> idDec( cellIDDesc, col ) ; + + for(int j=0,N= hCol->GetSize() ; j<N ; ++j) { + + lcio::SimTrackerHit* h = createSimTrackerHit( dynamic_cast<DD4hep::Simulation::Geant4TrackerHit*>( hCol->GetHit(j) ) ) ; + + col->addElement( h ) ; + +#if DEBUG + using namespace UTIL ; + std::cout << *h << G4endl ; +#endif + } + + lcEvt->addCollection( col , hCol->GetName() ) ; + } //----------------------------------------------------------------- + + + if( isCalorimeter ) { + + lcio::LCCollectionVec* col = new lcio::LCCollectionVec( lcio::LCIO::SIMCALORIMETERHIT ) ; + + UTIL::CellIDEncoder<SimCalorimeterHit> idDec( cellIDDesc, col ) ; + + col->setFlag( UTIL::make_bitset32( LCIO::CHBIT_LONG, LCIO::CHBIT_STEP ) ); + +#if DEBUG + std::cout << " setting collection flag: 0x" << std::hex << col->getFlag() << std::dec << std::endl ; +#endif + + for(int j=0,N= hCol->GetSize() ; j<N ; ++j) { + + lcio::SimCalorimeterHit* h = createSimCalorimeterHit( dynamic_cast<DD4hep::Simulation::Geant4CalorimeterHit*>( hCol->GetHit(j) ) ) ; + + col->addElement( h ) ; + +#if DEBUG + using namespace UTIL ; + std::cout << *h << G4endl ; +#endif + } + + lcEvt->addCollection( col , hCol->GetName() ) ; + } //----------------------------------------------------------------- + + } + } + + // --- write the event + lcEvt->setRunNumber( runAct->g4run->GetRunID() ) ; + lcEvt->setDetectorName( runAct->runData.detectorName ) ; + runAct->runData.lcioWriter->writeEvent( lcEvt ) ; + + +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/examples/ILDExSimu/src/ILDExEventActionMessenger.cpp b/examples/ILDExSimu/src/ILDExEventActionMessenger.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9089d329c10dbd7296fcd855629aa5a19c3f2280 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExEventActionMessenger.cpp @@ -0,0 +1,42 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExEventActionMessenger.h" +#include "ILDExEventAction.h" +#include "G4UIdirectory.hh" +#include "G4UIcmdWithAnInteger.hh" +#include "globals.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExEventActionMessenger::ILDExEventActionMessenger(ILDExEventAction* EvAct) +:eventAction(EvAct) +{ + eventDir = new G4UIdirectory("/ILDExDir/event/"); + eventDir->SetGuidance("event control"); + + PrintCmd = new G4UIcmdWithAnInteger("/ILDExDir/event/printModulo",this); + PrintCmd->SetGuidance("Print events modulo n"); + PrintCmd->SetParameterName("EventNb",false); + PrintCmd->SetRange("EventNb>0"); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExEventActionMessenger::~ILDExEventActionMessenger() +{ + delete PrintCmd; + delete eventDir; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExEventActionMessenger::SetNewValue( + G4UIcommand* command,G4String newValue) +{ + if(command == PrintCmd) + {eventAction->SetPrintModulo(PrintCmd->GetNewIntValue(newValue));} +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/examples/ILDExSimu/src/ILDExPhysicsList.cpp b/examples/ILDExSimu/src/ILDExPhysicsList.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3f7ec58e8637ce3a966341e8001bca196f6d96a9 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExPhysicsList.cpp @@ -0,0 +1,194 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExPhysicsList.h" + +#include "G4ProcessManager.hh" +#include "G4BosonConstructor.hh" +#include "G4LeptonConstructor.hh" +#include "G4MesonConstructor.hh" +#include "G4BosonConstructor.hh" +#include "G4BaryonConstructor.hh" +#include "G4IonConstructor.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExPhysicsList::ILDExPhysicsList(): G4VUserPhysicsList() +{ + defaultCutValue = 1.0*mm; + SetVerboseLevel(1); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExPhysicsList::~ILDExPhysicsList() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExPhysicsList::ConstructParticle() +{ + // In this method, static member functions should be called + // for all particles which you want to use. + // This ensures that objects of these particle types will be + // created in the program. + + G4BosonConstructor pBosonConstructor; + pBosonConstructor.ConstructParticle(); + + G4LeptonConstructor pLeptonConstructor; + pLeptonConstructor.ConstructParticle(); + +// G4MesonConstructor pMesonConstructor; +// pMesonConstructor.ConstructParticle(); +// +// G4BaryonConstructor pBaryonConstructor; +// pBaryonConstructor.ConstructParticle(); +// +// G4IonConstructor pIonConstructor; +// pIonConstructor.ConstructParticle(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExPhysicsList::ConstructProcess() +{ + AddTransportation(); + ConstructEM(); + ConstructDecay(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4ComptonScattering.hh" +#include "G4GammaConversion.hh" +#include "G4PhotoElectricEffect.hh" + +#include "G4eMultipleScattering.hh" +#include "G4eIonisation.hh" +#include "G4eBremsstrahlung.hh" +#include "G4eplusAnnihilation.hh" + +#include "G4MuMultipleScattering.hh" +#include "G4MuIonisation.hh" +#include "G4MuBremsstrahlung.hh" +#include "G4MuPairProduction.hh" + +#include "G4hMultipleScattering.hh" +#include "G4hIonisation.hh" +#include "G4hBremsstrahlung.hh" +#include "G4hPairProduction.hh" + +#include "G4ionIonisation.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExPhysicsList::ConstructEM() +{ + theParticleIterator->reset(); + while( (*theParticleIterator)() ){ + G4ParticleDefinition* particle = theParticleIterator->value(); + G4ProcessManager* pmanager = particle->GetProcessManager(); + G4String particleName = particle->GetParticleName(); + + if (particleName == "gamma") { + // gamma +// pmanager->AddDiscreteProcess(new G4PhotoElectricEffect); +// pmanager->AddDiscreteProcess(new G4ComptonScattering); +// pmanager->AddDiscreteProcess(new G4GammaConversion); + + } else if (particleName == "e-") { + //electron + pmanager->AddProcess(new G4eMultipleScattering,-1, 1, 1); + pmanager->AddProcess(new G4eIonisation, -1, 2, 2); +// pmanager->AddProcess(new G4eBremsstrahlung, -1, 3, 3); + + } else if (particleName == "e+") { + //positron + pmanager->AddProcess(new G4eMultipleScattering,-1, 1, 1); + pmanager->AddProcess(new G4eIonisation, -1, 2, 2); +// pmanager->AddProcess(new G4eBremsstrahlung, -1, 3, 3); +// pmanager->AddProcess(new G4eplusAnnihilation, 0,-1, 4); + + } else if( particleName == "mu+" || + particleName == "mu-" ) { + //muon + pmanager->AddProcess(new G4MuMultipleScattering,-1, 1, 1); + pmanager->AddProcess(new G4MuIonisation, -1, 2, 2); + // pmanager->AddProcess(new G4MuBremsstrahlung, -1, 3, 3); + // pmanager->AddProcess(new G4MuPairProduction, -1, 4, 4); + + } else if( particleName == "proton" || + particleName == "pi-" || + particleName == "pi+" ) { + //proton +// pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1); +// pmanager->AddProcess(new G4hIonisation, -1, 2, 2); +// pmanager->AddProcess(new G4hBremsstrahlung, -1, 3, 3); +// pmanager->AddProcess(new G4hPairProduction, -1, 4, 4); + + } else if( particleName == "alpha" || + particleName == "He3" ) { + //alpha +// pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1); +// pmanager->AddProcess(new G4ionIonisation, -1, 2, 2); + + } else if( particleName == "GenericIon" ) { + //Ions +// pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1); +// pmanager->AddProcess(new G4ionIonisation, -1, 2, 2); + + } else if ((!particle->IsShortLived()) && + (particle->GetPDGCharge() != 0.0) && + (particle->GetParticleName() != "chargedgeantino")) { + //all others charged particles except geantino +// pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1); +// pmanager->AddProcess(new G4hIonisation, -1, 2, 2); + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4Decay.hh" + +void ILDExPhysicsList::ConstructDecay() +{ + // Add Decay Process + G4Decay* theDecayProcess = new G4Decay(); + theParticleIterator->reset(); + while( (*theParticleIterator)() ){ + G4ParticleDefinition* particle = theParticleIterator->value(); + G4ProcessManager* pmanager = particle->GetProcessManager(); + if (theDecayProcess->IsApplicable(*particle)) { + pmanager ->AddProcess(theDecayProcess); + // set ordering for PostStepDoIt and AtRestDoIt + pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep); + pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest); + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExPhysicsList::SetCuts() +{ + if (verboseLevel >0){ + G4cout << "ILDExPhysicsList::SetCuts:"; + G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl; + } + + // set cut values for gamma at first and for e- second and next for e+, + // because some processes for e+/e- need cut values for gamma + // + SetCutValue(defaultCutValue, "gamma"); + SetCutValue(defaultCutValue, "e-"); + SetCutValue(defaultCutValue, "e+"); + // SetCutValue(defaultCutValue, "proton"); + + if (verboseLevel>0) DumpCutValuesTable(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + diff --git a/examples/ILDExSimu/src/ILDExPrimaryGeneratorAction.cpp b/examples/ILDExSimu/src/ILDExPrimaryGeneratorAction.cpp new file mode 100644 index 0000000000000000000000000000000000000000..890cbcd010f6890d66ed9dfa52e5860c7d34fb71 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExPrimaryGeneratorAction.cpp @@ -0,0 +1,62 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExPrimaryGeneratorAction.h" +#include "DD4hep/Volumes.h" +#include "TGeoBBox.h" +//#include "ILDExPrimaryGeneratorMessenger.h" + +#include "G4Event.hh" +#include "G4ParticleGun.hh" +#include "G4ParticleTable.hh" +#include "G4ParticleDefinition.hh" +#include "Randomize.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExPrimaryGeneratorAction::ILDExPrimaryGeneratorAction(const DD4hep::Geometry::LCDD& ILDExDC) +:ILDExDetector(ILDExDC),rndmFlag("off") +{ + G4int n_particle = 1; + particleGun = new G4ParticleGun(n_particle); + + // default particle kinematic + + G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); + G4String particleName; + G4ParticleDefinition* particle + = particleTable->FindParticle(particleName="e-"); + particleGun->SetParticleDefinition(particle); + particleGun->SetParticleMomentumDirection(G4ThreeVector(1.,0.,0.)); + particleGun->SetParticleEnergy(50.*MeV); + DD4hep::Geometry::Box worldbox = ILDExDetector.worldVolume().solid(); + G4double position = -0.5*(worldbox->GetDX()); + particleGun->SetParticlePosition(G4ThreeVector(position,0.*cm,0.*cm)); + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExPrimaryGeneratorAction::~ILDExPrimaryGeneratorAction() +{ + delete particleGun; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) +{ + //this function is called at the begining of event + // + G4double x0 = 0.*cm; + G4double y0 = 0.*cm; + G4double z0 = 0.*cm; + + particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0)); + + particleGun->GeneratePrimaryVertex(anEvent); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + diff --git a/examples/ILDExSimu/src/ILDExRunAction.cpp b/examples/ILDExSimu/src/ILDExRunAction.cpp new file mode 100644 index 0000000000000000000000000000000000000000..645ee75dd52394d837bd948b5e3533a68dbe5101 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExRunAction.cpp @@ -0,0 +1,134 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExRunAction.h" + +#include "G4Run.hh" +#include "G4RunManager.hh" +#include "G4UnitsTable.hh" + +//--- lcio -- +#include "lcio.h" +#include "IMPL/LCRunHeaderImpl.h" + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExRunAction::ILDExRunAction() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExRunAction::~ILDExRunAction() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExRunAction::BeginOfRunAction(const G4Run* aRun) +{ + G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; + + // keep the run around: + g4run = aRun ; + + //inform the runManager to save random number seed + G4RunManager::GetRunManager()->SetRandomNumberStore(true); + + + //initialize cumulative quantities + // + + sumESupport = sum2ESupport = sumESensitive = sum2ESensitive = 0.; + sumLSupport = sum2LSupport = sumLSensitive = sum2LSensitive = 0.; + sumAngleSupport = sum2AngleSupport = sumAngleSensitive = sum2AngleSensitive = 0.; + + + // --- write an lcio::RunHeader --------- + lcio::LCRunHeaderImpl* rh = new lcio::LCRunHeaderImpl ; + rh->setRunNumber( aRun->GetRunID() ) ; + rh->setDetectorName( runData.detectorName ) ; + runData.lcioWriter->writeRunHeader( rh ) ; + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExRunAction::fillPerEvent(G4double ESupport, G4double ESensitive, + G4double LSupport, G4double LSensitive, + G4double AngleSupport, G4double AngleSensitive) +{ + //accumulate statistic + // + sumESupport += ESupport; sum2ESupport += ESupport*ESupport; + sumESensitive += ESensitive; sum2ESensitive += ESensitive*ESensitive; + + sumLSupport += LSupport; sum2LSupport += LSupport*LSupport; + sumLSensitive += LSensitive; sum2LSensitive += LSensitive*LSensitive; + + sumAngleSupport += AngleSupport; sum2AngleSupport += AngleSupport*AngleSupport; + sumAngleSensitive += AngleSensitive; sum2AngleSensitive += AngleSensitive*AngleSensitive; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExRunAction::EndOfRunAction(const G4Run* aRun) +{ + G4int NbOfEvents = aRun->GetNumberOfEvent(); + if (NbOfEvents == 0) return; + + //compute statistics: mean and rms + // + sumESupport /= NbOfEvents; sum2ESupport /= NbOfEvents; + G4double rmsESupport = sum2ESupport - sumESupport*sumESupport; + if (rmsESupport >0.) rmsESupport = std::sqrt(rmsESupport); else rmsESupport = 0.; + + sumESensitive /= NbOfEvents; sum2ESensitive /= NbOfEvents; + G4double rmsESensitive = sum2ESensitive - sumESensitive*sumESensitive; + if (rmsESensitive >0.) rmsESensitive = std::sqrt(rmsESensitive); else rmsESensitive = 0.; + + sumLSupport /= NbOfEvents; sum2LSupport /= NbOfEvents; + G4double rmsLSupport = sum2LSupport - sumLSupport*sumLSupport; + if (rmsLSupport >0.) rmsLSupport = std::sqrt(rmsLSupport); else rmsLSupport = 0.; + + sumLSensitive /= NbOfEvents; sum2LSensitive /= NbOfEvents; + G4double rmsLSensitive = sum2LSensitive - sumLSensitive*sumLSensitive; + if (rmsLSensitive >0.) rmsLSensitive = std::sqrt(rmsLSensitive); else rmsLSensitive = 0.; + + sumAngleSupport /= NbOfEvents; sum2AngleSupport /= NbOfEvents; + G4double rmsAngleSupport = sum2AngleSupport - sumAngleSupport*sumAngleSupport; + if (rmsAngleSupport >0.) rmsAngleSupport = std::sqrt(rmsAngleSupport); else rmsAngleSupport = 0.; + + sumAngleSensitive /= NbOfEvents; sum2AngleSensitive /= NbOfEvents; + G4double rmsAngleSensitive = sum2AngleSensitive - sumAngleSensitive*sumAngleSensitive; + if (rmsAngleSensitive >0.) rmsAngleSensitive = std::sqrt(rmsAngleSensitive); else rmsAngleSensitive = 0.; + + + //print + // + G4cout + << "\n--------------------End of Run------------------------------\n" + << "\n mean Energy in Support : " << G4BestUnit(sumESupport,"Energy") + << " +- " << G4BestUnit(rmsESupport,"Energy") + << "\n mean Energy in Sensitive: " << G4BestUnit(sumESensitive,"Energy") + << " +- " << G4BestUnit(rmsESensitive,"Energy") + << G4endl; + + G4cout + << "\n mean trackLength in Support : " << G4BestUnit(sumLSupport,"Length") + << " +- " << G4BestUnit(rmsLSupport,"Length") + << "\n mean trackLength in Sensitive : " << G4BestUnit(sumLSensitive,"Length") + << " +- " << G4BestUnit(rmsLSensitive,"Length") + << G4endl; + + G4cout + << "\n mean Angle in Support : " << G4BestUnit(sumAngleSupport,"Angle") + << " +- " << G4BestUnit(rmsAngleSupport,"Angle") + << "\n mean Angle in Sensitive: " << G4BestUnit(sumAngleSensitive,"Angle") + << " +- " << G4BestUnit(rmsAngleSensitive,"Angle") + << "\n------------------------------------------------------------\n" + << G4endl; + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/examples/ILDExSimu/src/ILDExSteppingAction.cpp b/examples/ILDExSimu/src/ILDExSteppingAction.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c9d097e0e5eb5db6e21523daea742a5432196273 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExSteppingAction.cpp @@ -0,0 +1,52 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExSteppingAction.h" +#include "ILDExEventAction.h" + +#include "G4Step.hh" + +////#include "G4RunManager.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExSteppingAction::ILDExSteppingAction(ILDExEventAction* evt) +:eventaction(evt) +{ } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExSteppingAction::~ILDExSteppingAction() +{ } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExSteppingAction::UserSteppingAction(const G4Step* aStep) +{ + + static G4Material* SiMaterial = G4Material::GetMaterial("Silicon"); + static G4Material* TPCGasMaterial = G4Material::GetMaterial("Argon"); + + // get volume of the current step + G4VPhysicalVolume* volume = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume(); + G4Material* material = volume->GetLogicalVolume()->GetMaterial(); + + // collect energy and track length step by step + G4double edep = aStep->GetTotalEnergyDeposit(); + + G4double stepl = 0.; + if (aStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) stepl = aStep->GetStepLength(); + + if (material == SiMaterial || material == TPCGasMaterial) { + eventaction->SumSensitive(edep, stepl, 0.0); + } + else { + eventaction->SumSupport(edep, stepl, 0.0); + } + + //example of saving random number seed of this event, under condition + //// if (condition) G4RunManager::GetRunManager()->rndmSaveThisEvent(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/examples/ILDExSimu/src/ILDExSteppingVerbose.cpp b/examples/ILDExSimu/src/ILDExSteppingVerbose.cpp new file mode 100644 index 0000000000000000000000000000000000000000..33fc810296d0636177d944fe6d4e4daf90000410 --- /dev/null +++ b/examples/ILDExSimu/src/ILDExSteppingVerbose.cpp @@ -0,0 +1,154 @@ + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "ILDExSteppingVerbose.h" + +#include "G4SteppingManager.hh" +#include "G4UnitsTable.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExSteppingVerbose::ILDExSteppingVerbose() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +ILDExSteppingVerbose::~ILDExSteppingVerbose() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExSteppingVerbose::StepInfo() +{ + CopyState(); + + G4int prec = G4cout.precision(3); + + if( verboseLevel >= 1 ){ + if( verboseLevel >= 4 ) VerboseTrack(); + if( verboseLevel >= 3 ){ + G4cout << G4endl; + G4cout << std::setw( 5) << "#Step#" << " " + << std::setw( 6) << "X" << " " + << std::setw( 6) << "Y" << " " + << std::setw( 6) << "Z" << " " + << std::setw( 9) << "KineE" << " " + << std::setw( 9) << "dEStep" << " " + << std::setw(10) << "StepLeng" + << std::setw(10) << "TrakLeng" + << std::setw(10) << "Volume" << " " + << std::setw(10) << "Process" << G4endl; + } + + G4cout << std::setw(5) << fTrack->GetCurrentStepNumber() << " " + << std::setw(6) << G4BestUnit(fTrack->GetPosition().x(),"Length") + << std::setw(6) << G4BestUnit(fTrack->GetPosition().y(),"Length") + << std::setw(6) << G4BestUnit(fTrack->GetPosition().z(),"Length") + << std::setw(6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy") + << std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy") + << std::setw(6) << G4BestUnit(fStep->GetStepLength(),"Length") + << std::setw(6) << G4BestUnit(fTrack->GetTrackLength(),"Length") + << " "; + + // if( fStepStatus != fWorldBoundary){ + if( fTrack->GetNextVolume() != 0 ) { + G4cout << std::setw(10) << fTrack->GetVolume()->GetName(); + } else { + G4cout << std::setw(10) << "OutOfWorld"; + } + + if(fStep->GetPostStepPoint()->GetProcessDefinedStep() != 0){ + G4cout << " " + << std::setw(10) + << fStep->GetPostStepPoint()->GetProcessDefinedStep() + ->GetProcessName(); + } else { + G4cout << " UserLimit"; + } + + G4cout << G4endl; + + if( verboseLevel == 2 ){ + G4int tN2ndariesTot = fN2ndariesAtRestDoIt + + fN2ndariesAlongStepDoIt + + fN2ndariesPostStepDoIt; + if(tN2ndariesTot>0){ + G4cout << " :----- List of 2ndaries - " + << "#SpawnInStep=" << std::setw(3) << tN2ndariesTot + << "(Rest=" << std::setw(2) << fN2ndariesAtRestDoIt + << ",Along=" << std::setw(2) << fN2ndariesAlongStepDoIt + << ",Post=" << std::setw(2) << fN2ndariesPostStepDoIt + << "), " + << "#SpawnTotal=" << std::setw(3) << (*fSecondary).size() + << " ---------------" + << G4endl; + + for(size_t lp1=(*fSecondary).size()-tN2ndariesTot; + lp1<(*fSecondary).size(); lp1++){ + G4cout << " : " + << std::setw(6) + << G4BestUnit((*fSecondary)[lp1]->GetPosition().x(),"Length") + << std::setw(6) + << G4BestUnit((*fSecondary)[lp1]->GetPosition().y(),"Length") + << std::setw(6) + << G4BestUnit((*fSecondary)[lp1]->GetPosition().z(),"Length") + << std::setw(6) + << G4BestUnit((*fSecondary)[lp1]->GetKineticEnergy(),"Energy") + << std::setw(10) + << (*fSecondary)[lp1]->GetDefinition()->GetParticleName(); + G4cout << G4endl; + } + + G4cout << " :-----------------------------" + << "----------------------------------" + << "-- EndOf2ndaries Info ---------------" + << G4endl; + } + } + + } + G4cout.precision(prec); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void ILDExSteppingVerbose::TrackingStarted() +{ + + CopyState(); +G4int prec = G4cout.precision(3); + if( verboseLevel > 0 ){ + + G4cout << std::setw( 5) << "Step#" << " " + << std::setw( 6) << "X" << " " + << std::setw( 6) << "Y" << " " + << std::setw( 6) << "Z" << " " + << std::setw( 9) << "KineE" << " " + << std::setw( 9) << "dEStep" << " " + << std::setw(10) << "StepLeng" + << std::setw(10) << "TrakLeng" + << std::setw(10) << "Volume" << " " + << std::setw(10) << "Process" << G4endl; + + G4cout << std::setw( 5) << fTrack->GetCurrentStepNumber() << " " + << std::setw( 6) << G4BestUnit(fTrack->GetPosition().x(),"Length") + << std::setw( 6) << G4BestUnit(fTrack->GetPosition().y(),"Length") + << std::setw( 6) << G4BestUnit(fTrack->GetPosition().z(),"Length") + << std::setw( 6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy") + << std::setw( 6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy") + << std::setw( 6) << G4BestUnit(fStep->GetStepLength(),"Length") + << std::setw( 6) << G4BestUnit(fTrack->GetTrackLength(),"Length") + << " "; + + if(fTrack->GetNextVolume()){ + G4cout << std::setw(10) << fTrack->GetVolume()->GetName(); + } else { + G4cout << "OutOfWorld"; + } + G4cout << " initStep" << G4endl; + } + G4cout.precision(prec); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......