From b9c1117b08badbb5f68cecaa90394e288f8773df Mon Sep 17 00:00:00 2001 From: Lorenzo Pezzotti Date: Thu, 15 Jul 2021 15:52:54 +0200 Subject: [PATCH 1/6] first impelementation of FullOptic parser option --- + | 141 +++++++++++++++++++ DREMTubes.cc | 14 +- include/DREMTubesActionInitialization.hh | 7 +- include/DREMTubesSteppingAction.hh | 9 +- src/DREMTubesActionInitialization.cc | 8 +- src/DREMTubesSteppingAction.cc | 166 ++++++++++++----------- 6 files changed, 256 insertions(+), 89 deletions(-) create mode 100644 + diff --git a/+ b/+ new file mode 100644 index 0000000..83a7efa --- /dev/null +++ b/+ @@ -0,0 +1,141 @@ +//************************************************** +// \file DREMTubes.cc +// \brief: main() of DREMTubes project +// \author: Lorenzo Pezzotti (CERN EP-SFT-sim) @lopezzot +// \start date: 7 July 2021 +//************************************************** + +// Includers from project files +// +#include "DREMTubesDetectorConstruction.hh" +#include "DREMTubesActionInitialization.hh" +#include "DREMTubesPhysicsList.hh" + +// Includers from Geant4 +// +#ifdef G4MULTITHREADED +#include "G4MTRunManager.hh" +#else +#include "G4RunManager.hh" +#endif +#include "G4UImanager.hh" +#include "G4UIcommand.hh" +#include "FTFP_BERT.hh" +#include "Randomize.hh" +#include "G4VisExecutive.hh" +#include "G4UIExecutive.hh" + +// G4err output for usage error +// +namespace PrintUsageError { + void UsageError() { + G4cerr << "->DREMTubes usage: " << G4endl; + G4cerr << "DREMTubes [-m macro ] [-u UIsession] [-t nThreads] [-pl PhysicsList]" + << G4endl; + G4cerr << " [-opt FullOptic]" << G4endl; + } +} + +// main() function +// +int main(int argc, char** argv) { + + // Error in argument numbers + // + if ( argc > 11 ) { + PrintUsageError::UsageError(); + return 1; + } + + // Convert arguments in G4string and G4int + // + G4String macro; + G4String session; + G4String custom_pl = "FTFP_BERT"; //default physics list + G4bool FullOptic = false; + #ifdef G4MULTITHREADED + G4int nThreads = 0; + #endif + + for ( G4int i=1; i Run with full optical description"< 0 ) { + runManager->SetNumberOfThreads(nThreads); + } + #else + G4RunManager * runManager = new G4RunManager; + #endif + + // Set mandatory initialization classes + // + auto DetConstruction = new DREMTubesDetectorConstruction(); + runManager->SetUserInitialization(DetConstruction); + + runManager->SetUserInitialization(new DREMTubesPhysicsList(custom_pl)); + + auto actionInitialization = new DREMTubesActionInitialization(); + runManager->SetUserInitialization(actionInitialization); + + // Initialize visualization + // + auto visManager = new G4VisExecutive("Quiet"); + visManager->Initialize(); + + // Get the pointer to the User Interface manager + auto UImanager = G4UImanager::GetUIpointer(); + + // Process macro or start UI session + // + if ( macro.size() ) { //macro card mode + G4String command = "/control/execute "; + UImanager->ApplyCommand(command+macro); + } + else { //start UI session + UImanager->ApplyCommand("/control/execute DREMTubes_init_vis.mac"); + if (ui->IsGUI()) { + UImanager->ApplyCommand("/control/execute DREMTubes_gui.mac"); + } + ui->SessionStart(); + delete ui; + } + + // Program termination (user actions deleted by run manager) + // + delete visManager; + delete runManager; + +} + +//************************************************** diff --git a/DREMTubes.cc b/DREMTubes.cc index 5196e12..7eda123 100644 --- a/DREMTubes.cc +++ b/DREMTubes.cc @@ -32,6 +32,7 @@ namespace PrintUsageError { G4cerr << "->DREMTubes usage: " << G4endl; G4cerr << "DREMTubes [-m macro ] [-u UIsession] [-t nThreads] [-pl PhysicsList]" << G4endl; + G4cerr << " [-opt FullOptic]" << G4endl; } } @@ -41,7 +42,7 @@ int main(int argc, char** argv) { // Error in argument numbers // - if ( argc > 9 ) { + if ( argc > 11 ) { PrintUsageError::UsageError(); return 1; } @@ -51,6 +52,7 @@ int main(int argc, char** argv) { G4String macro; G4String session; G4String custom_pl = "FTFP_BERT"; //default physics list + G4bool FullOptic = false; #ifdef G4MULTITHREADED G4int nThreads = 0; #endif @@ -59,6 +61,8 @@ int main(int argc, char** argv) { if ( G4String(argv[i]) == "-m" ) macro = argv[i+1]; else if ( G4String(argv[i]) == "-u" ) session = argv[i+1]; else if ( G4String(argv[i]) == "-pl") custom_pl = argv[i+1]; + else if ( G4String(argv[i]) == "-opt") FullOptic = + G4UIcommand::ConvertToBool(argv[i+1]); #ifdef G4MULTITHREADED else if ( G4String(argv[i]) == "-t" ) { nThreads = G4UIcommand::ConvertToInt(argv[i+1]); @@ -68,7 +72,11 @@ int main(int argc, char** argv) { PrintUsageError::UsageError(); return 1; } - } + } + + //Print if FullOptic option is on + // + if (FullOptic){ G4cout<<"DREMTubes-> Run with full optical description"<SetUserInitialization(new DREMTubesPhysicsList(custom_pl)); - auto actionInitialization = new DREMTubesActionInitialization(); + auto actionInitialization = new DREMTubesActionInitialization( FullOptic ); runManager->SetUserInitialization(actionInitialization); // Initialize visualization diff --git a/include/DREMTubesActionInitialization.hh b/include/DREMTubesActionInitialization.hh index 761c5f2..42cd7bc 100644 --- a/include/DREMTubesActionInitialization.hh +++ b/include/DREMTubesActionInitialization.hh @@ -13,18 +13,23 @@ //Includers from Geant4 // #include "G4VUserActionInitialization.hh" +#include "G4Types.hh" class DREMTubesActionInitialization : public G4VUserActionInitialization { public: //Constructor // - DREMTubesActionInitialization(); + DREMTubesActionInitialization( const G4bool FullOptic ); virtual ~DREMTubesActionInitialization(); virtual void BuildForMaster() const; virtual void Build() const; + private: + + G4bool fFullOptic; + }; #endif diff --git a/include/DREMTubesSteppingAction.hh b/include/DREMTubesSteppingAction.hh index d2d32fe..4a2ca76 100644 --- a/include/DREMTubesSteppingAction.hh +++ b/include/DREMTubesSteppingAction.hh @@ -13,6 +13,7 @@ //Includers from Geant4 // #include "G4UserSteppingAction.hh" +#include "G4Types.hh" class G4OpBoundaryProcess; class DREMTubesDetectorConstruction; @@ -23,14 +24,15 @@ class DREMTubesSteppingAction : public G4UserSteppingAction { public: //Constructor // - DREMTubesSteppingAction(DREMTubesEventAction* eventAction); + DREMTubesSteppingAction(DREMTubesEventAction* eventAction, + const G4bool FullOptic ); //De-constructor // virtual ~DREMTubesSteppingAction(); //User impementation of SteppingAction // - virtual void UserSteppingAction( const G4Step* step ); + virtual void UserSteppingAction( const G4Step* step); private: //Data members @@ -38,7 +40,8 @@ class DREMTubesSteppingAction : public G4UserSteppingAction { DREMTubesEventAction* fEventAction; G4OpBoundaryProcess* fOpProcess; - + + G4bool fFullOptic; }; #endif diff --git a/src/DREMTubesActionInitialization.cc b/src/DREMTubesActionInitialization.cc index de4d7e5..d268a6e 100644 --- a/src/DREMTubesActionInitialization.cc +++ b/src/DREMTubesActionInitialization.cc @@ -15,8 +15,10 @@ //Constructor // -DREMTubesActionInitialization::DREMTubesActionInitialization() - : G4VUserActionInitialization() {} +DREMTubesActionInitialization::DREMTubesActionInitialization( const G4bool FullOptic ) + : G4VUserActionInitialization(), + fFullOptic( FullOptic ) +{} //De-constructor // @@ -39,7 +41,7 @@ void DREMTubesActionInitialization::Build() const { auto eventAction = new DREMTubesEventAction; SetUserAction(new DREMTubesRunAction( eventAction )); SetUserAction(eventAction); - SetUserAction(new DREMTubesSteppingAction(eventAction)); + SetUserAction(new DREMTubesSteppingAction(eventAction, fFullOptic)); } diff --git a/src/DREMTubesSteppingAction.cc b/src/DREMTubesSteppingAction.cc index 94a0816..29fc1ff 100644 --- a/src/DREMTubesSteppingAction.cc +++ b/src/DREMTubesSteppingAction.cc @@ -26,9 +26,11 @@ //Define constructor // DREMTubesSteppingAction::DREMTubesSteppingAction( - DREMTubesEventAction* eventAction) + DREMTubesEventAction* eventAction, + const G4bool FullOptic) : G4UserSteppingAction(), - fEventAction(eventAction) + fEventAction(eventAction), + fFullOptic(FullOptic) {} //Define de-constructor @@ -37,8 +39,8 @@ DREMTubesSteppingAction::~DREMTubesSteppingAction() {} //Define UserSteppingAction() method // -void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step ) { - +void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step) { + //Random seed and random number generator // unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); @@ -114,88 +116,94 @@ void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step ) { //Store information from Scintillation and Cherenkov //signals //-------------------------------------------------- - - std::string Fiber; - std::string S_fiber = "S_fiber"; - std::string C_fiber = "C_fiber"; - Fiber = PreStepVolume->GetName(); //name of current step fiber - - G4int copynumber; - - if ( strstr( Fiber.c_str(), S_fiber.c_str() ) ) { //scintillating fiber - G4double saturatedenergydeposited = 0.; - copynumber = step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); - if(step->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) { - if (steplength != 0) { - saturatedenergydeposited = - (energydeposited/steplength) / - ( 1+k_B*(energydeposited/steplength) ) * steplength; + + if (!fFullOptic){ + std::string Fiber; + std::string S_fiber = "S_fiber"; + std::string C_fiber = "C_fiber"; + Fiber = PreStepVolume->GetName(); //name of current step fiber + + G4int copynumber; + + if ( strstr( Fiber.c_str(), S_fiber.c_str() ) ) { //scintillating fiber + G4double saturatedenergydeposited = 0.; + copynumber = step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); + if(step->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) { + if (steplength != 0) { + saturatedenergydeposited = + (energydeposited/steplength) / + ( 1+k_B*(energydeposited/steplength) ) * steplength; + } } + fEventAction->AddScin(energydeposited); + std::poisson_distribution scin_distribution( + saturatedenergydeposited*3.78); + int s_signal = scin_distribution(generator); + fEventAction->AddVectorScinEnergy(s_signal,copynumber); + //energy deposited in any scintillating fiber (saturated) } - fEventAction->AddScin(energydeposited); - std::poisson_distribution scin_distribution(saturatedenergydeposited*3.78); - int s_signal = scin_distribution(generator); - fEventAction->AddVectorScinEnergy(s_signal,copynumber); - //energy deposited in any scintillating fiber (saturated) - } - - if ( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fiber - fEventAction->AddCher(energydeposited); - } - - G4OpBoundaryProcessStatus theStatus = Undefined; - - G4ProcessManager* OpManager = - G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); - - if (OpManager) { - G4int MAXofPostStepLoops = - OpManager->GetPostStepProcessVector()->entries(); - G4ProcessVector* fPostStepDoItVector = - OpManager->GetPostStepProcessVector(typeDoIt); - - for ( G4int i=0; i(fCurrentProcess); - if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;} + if ( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fiber + fEventAction->AddCher(energydeposited); } - } - - if( particlename == "opticalphoton" ) { //optical photons - - switch ( theStatus ){ - - case TotalInternalReflection: //photon reflected inside fiber - Fiber = PreStepVolume->GetName(); - - if(strstr(Fiber.c_str(),S_fiber.c_str())){ //scintillating fibre - step->GetTrack()->SetTrackStatus(fStopAndKill); - } - - if( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fibre - copynumber = - step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); - int c_signal = cher_distribution(generator); - fEventAction->AddCherenkov(c_signal); - fEventAction->AddVectorCherPE(copynumber, c_signal); - step->GetTrack()->SetTrackStatus(fStopAndKill); + + G4OpBoundaryProcessStatus theStatus = Undefined; + + G4ProcessManager* OpManager = + G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); + + if (OpManager) { + G4int MAXofPostStepLoops = + OpManager->GetPostStepProcessVector()->entries(); + G4ProcessVector* fPostStepDoItVector = + OpManager->GetPostStepProcessVector(typeDoIt); + + for ( G4int i=0; i(fCurrentProcess); + if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;} } - break; - - case Detection: - //To be used for SiPM simulation (optional) - // - break; - - default: - step->GetTrack()->SetTrackStatus(fStopAndKill); - break; - - } //end of switch cases + } + + if( particlename == "opticalphoton" ) { //optical photons + + switch ( theStatus ){ + + case TotalInternalReflection: //photon reflected inside fiber + Fiber = PreStepVolume->GetName(); + + if(strstr(Fiber.c_str(),S_fiber.c_str())){ //scintillating fibre + step->GetTrack()->SetTrackStatus(fStopAndKill); + } + + if( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fibre + copynumber = + step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); + int c_signal = cher_distribution(generator); + fEventAction->AddCherenkov(c_signal); + fEventAction->AddVectorCherPE(copynumber, c_signal); + step->GetTrack()->SetTrackStatus(fStopAndKill); + } + break; + + case Detection: + //To be used for SiPM simulation (optional) + // + break; + + default: + step->GetTrack()->SetTrackStatus(fStopAndKill); + break; + + } //end of switch cases + + }//end of optical photon loop - }//end of optical photon loop + }//end of FullOptic == false option + if (fFullOptic){ + if ( particlename == "opticalphoton" ){G4cout<<"fotone"< Date: Thu, 15 Jul 2021 16:35:04 +0200 Subject: [PATCH 2/6] remove + file --- + | 141 -------------------------------------------------------------- 1 file changed, 141 deletions(-) delete mode 100644 + diff --git a/+ b/+ deleted file mode 100644 index 83a7efa..0000000 --- a/+ +++ /dev/null @@ -1,141 +0,0 @@ -//************************************************** -// \file DREMTubes.cc -// \brief: main() of DREMTubes project -// \author: Lorenzo Pezzotti (CERN EP-SFT-sim) @lopezzot -// \start date: 7 July 2021 -//************************************************** - -// Includers from project files -// -#include "DREMTubesDetectorConstruction.hh" -#include "DREMTubesActionInitialization.hh" -#include "DREMTubesPhysicsList.hh" - -// Includers from Geant4 -// -#ifdef G4MULTITHREADED -#include "G4MTRunManager.hh" -#else -#include "G4RunManager.hh" -#endif -#include "G4UImanager.hh" -#include "G4UIcommand.hh" -#include "FTFP_BERT.hh" -#include "Randomize.hh" -#include "G4VisExecutive.hh" -#include "G4UIExecutive.hh" - -// G4err output for usage error -// -namespace PrintUsageError { - void UsageError() { - G4cerr << "->DREMTubes usage: " << G4endl; - G4cerr << "DREMTubes [-m macro ] [-u UIsession] [-t nThreads] [-pl PhysicsList]" - << G4endl; - G4cerr << " [-opt FullOptic]" << G4endl; - } -} - -// main() function -// -int main(int argc, char** argv) { - - // Error in argument numbers - // - if ( argc > 11 ) { - PrintUsageError::UsageError(); - return 1; - } - - // Convert arguments in G4string and G4int - // - G4String macro; - G4String session; - G4String custom_pl = "FTFP_BERT"; //default physics list - G4bool FullOptic = false; - #ifdef G4MULTITHREADED - G4int nThreads = 0; - #endif - - for ( G4int i=1; i Run with full optical description"< 0 ) { - runManager->SetNumberOfThreads(nThreads); - } - #else - G4RunManager * runManager = new G4RunManager; - #endif - - // Set mandatory initialization classes - // - auto DetConstruction = new DREMTubesDetectorConstruction(); - runManager->SetUserInitialization(DetConstruction); - - runManager->SetUserInitialization(new DREMTubesPhysicsList(custom_pl)); - - auto actionInitialization = new DREMTubesActionInitialization(); - runManager->SetUserInitialization(actionInitialization); - - // Initialize visualization - // - auto visManager = new G4VisExecutive("Quiet"); - visManager->Initialize(); - - // Get the pointer to the User Interface manager - auto UImanager = G4UImanager::GetUIpointer(); - - // Process macro or start UI session - // - if ( macro.size() ) { //macro card mode - G4String command = "/control/execute "; - UImanager->ApplyCommand(command+macro); - } - else { //start UI session - UImanager->ApplyCommand("/control/execute DREMTubes_init_vis.mac"); - if (ui->IsGUI()) { - UImanager->ApplyCommand("/control/execute DREMTubes_gui.mac"); - } - ui->SessionStart(); - delete ui; - } - - // Program termination (user actions deleted by run manager) - // - delete visManager; - delete runManager; - -} - -//************************************************** From e457cd45d6901b6a8418e9ca243becb43e7ac23b Mon Sep 17 00:00:00 2001 From: Lorenzo Pezzotti Date: Mon, 19 Jul 2021 21:19:10 +0200 Subject: [PATCH 3/6] towards full optics as parser option --- include/DREMTubesSteppingAction.hh | 16 +- src/DREMTubesSteppingAction.cc | 240 +++++++++++++++-------------- 2 files changed, 142 insertions(+), 114 deletions(-) diff --git a/include/DREMTubesSteppingAction.hh b/include/DREMTubesSteppingAction.hh index 4a2ca76..3713f83 100644 --- a/include/DREMTubesSteppingAction.hh +++ b/include/DREMTubesSteppingAction.hh @@ -32,7 +32,21 @@ class DREMTubesSteppingAction : public G4UserSteppingAction { //User impementation of SteppingAction // - virtual void UserSteppingAction( const G4Step* step); + virtual void UserSteppingAction( const G4Step* step ); + + //Retrieve auxialiry info from Step + // + void AuxSteppingAction( const G4Step* step ); + + //Fast signal simulation (no optical photon propagation) + //fFullOptic == false + // + void FastSteppingAction( const G4Step* step ); + + //Slow signal simulation (optical photon propagation) + //fFullOptic == true + // + void SlowSteppingAction( const G4Step* step ); private: //Data members diff --git a/src/DREMTubesSteppingAction.cc b/src/DREMTubesSteppingAction.cc index 29fc1ff..e6ccfdb 100644 --- a/src/DREMTubesSteppingAction.cc +++ b/src/DREMTubesSteppingAction.cc @@ -39,17 +39,35 @@ DREMTubesSteppingAction::~DREMTubesSteppingAction() {} //Define UserSteppingAction() method // -void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step) { +void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step ) { - //Random seed and random number generator + //Save auxiliary information // - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - std::default_random_engine generator(seed); - - //Shift from 200 Cp.e./GeV to 50 Cp.e./GeV + AuxSteppingAction( step ); + + //Save fast signal information // - std::poisson_distribution cher_distribution(0.155); + if( fFullOptic == false ){ + FastSteppingAction( step ); + } + //Save slow signal information + // + else { + SlowSteppingAction( step ); + } +} + +//Define SlowSteppingAction() method +// +void DREMTubesSteppingAction::SlowSteppingAction( const G4Step* step ){ + +} + +//Define AuxSteppingAction() method +// +void DREMTubesSteppingAction::AuxSteppingAction( const G4Step* step ) { + // Get step info // G4VPhysicalVolume* PreStepVolume @@ -58,20 +76,11 @@ void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step) { //= step->GetPostStepPoint()->GetTouchableHandle()->GetVolume(); G4double energydeposited = step->GetTotalEnergyDeposit(); G4String particlename = step->GetTrack()->GetDefinition()->GetParticleName(); - G4double steplength = step->GetStepLength(); - - double k_B = 0.126; //Birks constant - + //-------------------------------------------------- //Store auxiliary information from event steps //-------------------------------------------------- - /*if (PreStepVolume->GetName() == "module"){ - //Function to save absorber material name - fEventAction->SaveAbsorberMaterial - (PreStepVolume->GetLogicalVolume()->GetMaterial()->GetName()); - }*/ - if (PreStepVolume->GetName() == "leakageabsorber" ){ fEventAction->AddEscapedEnergy(step->GetTrack()->GetKineticEnergy()); step->GetTrack()->SetTrackStatus(fStopAndKill); @@ -89,121 +98,126 @@ void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step) { fEventAction->Addem(energydeposited); //energy deposited by em-component } } - - /*//em2 fraction estimation (deprecated) - // - if (particlename == "pi0" && step->GetTrack()->GetCurrentStepNumber() == 1){ - //if it's a neutral pion at first step - fEventAction->Addem2(step->GetTrack()->GetTotalEnergy()); - } - - if (particlename == "gamma" && - step->GetTrack()->GetCreatorProcess()->GetProcessName() != "eBrem" && - step->GetTrack()->GetCreatorProcess()->GetProcessName() != "annihil"){ - if(step->GetTrack()->GetCurrentStepNumber() == 1){ - fEventAction->Addem2(step->GetTrack()->GetTotalEnergy()); - } - }*/ - + if ( step->GetTrack()->GetTrackID() == 1 && step->GetTrack()->GetCurrentStepNumber() == 1){ //Save primary particle energy and name fEventAction->SavePrimaryParticle(particlename); fEventAction->SavePrimaryEnergy(step->GetTrack()->GetKineticEnergy()); } + +} + +//Define FastSteppingAction() method +// +void DREMTubesSteppingAction::FastSteppingAction( const G4Step* step ) { + + //Random seed and random number generator + // + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::default_random_engine generator(seed); + + //Shift from 200 Cp.e./GeV to 50 Cp.e./GeV + // + std::poisson_distribution cher_distribution(0.155); + + // Get step info + // + G4VPhysicalVolume* PreStepVolume + = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume(); + //G4VPhysicalVolume* PostStepVolume + //= step->GetPostStepPoint()->GetTouchableHandle()->GetVolume(); + G4double energydeposited = step->GetTotalEnergyDeposit(); + G4String particlename = step->GetTrack()->GetDefinition()->GetParticleName(); + G4double steplength = step->GetStepLength(); + double k_B = 0.126; //Birks constant //-------------------------------------------------- //Store information from Scintillation and Cherenkov //signals //-------------------------------------------------- - if (!fFullOptic){ - std::string Fiber; - std::string S_fiber = "S_fiber"; - std::string C_fiber = "C_fiber"; - Fiber = PreStepVolume->GetName(); //name of current step fiber - - G4int copynumber; - - if ( strstr( Fiber.c_str(), S_fiber.c_str() ) ) { //scintillating fiber - G4double saturatedenergydeposited = 0.; - copynumber = step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); - if(step->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) { - if (steplength != 0) { - saturatedenergydeposited = - (energydeposited/steplength) / - ( 1+k_B*(energydeposited/steplength) ) * steplength; - } + std::string Fiber; + std::string S_fiber = "S_fiber"; + std::string C_fiber = "C_fiber"; + Fiber = PreStepVolume->GetName(); //name of current step fiber + + G4int copynumber; + + if ( strstr( Fiber.c_str(), S_fiber.c_str() ) ) { //scintillating fiber + G4double saturatedenergydeposited = 0.; + copynumber = step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); + if(step->GetTrack()->GetDefinition()->GetPDGCharge() != 0.) { + if (steplength != 0) { + saturatedenergydeposited = + (energydeposited/steplength) / + ( 1+k_B*(energydeposited/steplength) ) * steplength; } - fEventAction->AddScin(energydeposited); - std::poisson_distribution scin_distribution( - saturatedenergydeposited*3.78); - int s_signal = scin_distribution(generator); - fEventAction->AddVectorScinEnergy(s_signal,copynumber); - //energy deposited in any scintillating fiber (saturated) - } - - if ( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fiber - fEventAction->AddCher(energydeposited); } - - G4OpBoundaryProcessStatus theStatus = Undefined; - - G4ProcessManager* OpManager = - G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); - - if (OpManager) { - G4int MAXofPostStepLoops = - OpManager->GetPostStepProcessVector()->entries(); - G4ProcessVector* fPostStepDoItVector = - OpManager->GetPostStepProcessVector(typeDoIt); + fEventAction->AddScin(energydeposited); + std::poisson_distribution scin_distribution( + saturatedenergydeposited*3.78); + int s_signal = scin_distribution(generator); + fEventAction->AddVectorScinEnergy(s_signal,copynumber); + //energy deposited in any scintillating fiber (saturated) + } + + if ( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fiber + fEventAction->AddCher(energydeposited); + } - for ( G4int i=0; i(fCurrentProcess); - if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;} - } + G4OpBoundaryProcessStatus theStatus = Undefined; + + G4ProcessManager* OpManager = + G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); + + if (OpManager) { + G4int MAXofPostStepLoops = + OpManager->GetPostStepProcessVector()->entries(); + G4ProcessVector* fPostStepDoItVector = + OpManager->GetPostStepProcessVector(typeDoIt); + + for ( G4int i=0; i(fCurrentProcess); + if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;} } - - if( particlename == "opticalphoton" ) { //optical photons - - switch ( theStatus ){ - - case TotalInternalReflection: //photon reflected inside fiber - Fiber = PreStepVolume->GetName(); - - if(strstr(Fiber.c_str(),S_fiber.c_str())){ //scintillating fibre - step->GetTrack()->SetTrackStatus(fStopAndKill); - } - - if( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fibre - copynumber = - step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); - int c_signal = cher_distribution(generator); - fEventAction->AddCherenkov(c_signal); - fEventAction->AddVectorCherPE(copynumber, c_signal); - step->GetTrack()->SetTrackStatus(fStopAndKill); - } - break; - - case Detection: - //To be used for SiPM simulation (optional) - // - break; - - default: + } + + if( particlename == "opticalphoton" ) { //optical photons + + switch ( theStatus ){ + + case TotalInternalReflection: //photon reflected inside fiber + Fiber = PreStepVolume->GetName(); + + if(strstr(Fiber.c_str(),S_fiber.c_str())){ //scintillating fibre step->GetTrack()->SetTrackStatus(fStopAndKill); - break; - - } //end of switch cases - - }//end of optical photon loop + } + + if( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fibre + copynumber = + step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(1); + int c_signal = cher_distribution(generator); + fEventAction->AddCherenkov(c_signal); + fEventAction->AddVectorCherPE(copynumber, c_signal); + step->GetTrack()->SetTrackStatus(fStopAndKill); + } + break; - }//end of FullOptic == false option + case Detection: + //To be used for SiPM simulation (optional) + // + break; + + default: + step->GetTrack()->SetTrackStatus(fStopAndKill); + break; + + } //end of switch cases + + }//end of optical photon loop - if (fFullOptic){ - if ( particlename == "opticalphoton" ){G4cout<<"fotone"< Date: Tue, 20 Jul 2021 10:41:40 +0200 Subject: [PATCH 4/6] adding FullOptic to OpticalPhysicsList --- DREMTubes.cc | 2 +- include/DREMTubesOpticalPhysics.hh | 4 +++- include/DREMTubesPhysicsList.hh | 7 ++++++- src/DREMTubesDetectorConstruction.cc | 6 +++--- src/DREMTubesOpticalPhysics.cc | 20 +++++++++++--------- src/DREMTubesPhysicsList.cc | 6 ++++-- 6 files changed, 28 insertions(+), 17 deletions(-) diff --git a/DREMTubes.cc b/DREMTubes.cc index 7eda123..4a7c0d1 100644 --- a/DREMTubes.cc +++ b/DREMTubes.cc @@ -101,7 +101,7 @@ int main(int argc, char** argv) { auto DetConstruction = new DREMTubesDetectorConstruction(); runManager->SetUserInitialization(DetConstruction); - runManager->SetUserInitialization(new DREMTubesPhysicsList(custom_pl)); + runManager->SetUserInitialization(new DREMTubesPhysicsList(custom_pl, FullOptic )); auto actionInitialization = new DREMTubesActionInitialization( FullOptic ); runManager->SetUserInitialization(actionInitialization); diff --git a/include/DREMTubesOpticalPhysics.hh b/include/DREMTubesOpticalPhysics.hh index 9ba5b7e..6d444b9 100644 --- a/include/DREMTubesOpticalPhysics.hh +++ b/include/DREMTubesOpticalPhysics.hh @@ -27,7 +27,7 @@ class DREMTubesOpticalPhysics : public G4VPhysicsConstructor { public: // Constructor // - DREMTubesOpticalPhysics(G4bool toggle=true); + DREMTubesOpticalPhysics( const G4bool FullOptic, G4bool toggle=true); // Deconstructor // virtual ~DREMTubesOpticalPhysics(); @@ -56,6 +56,8 @@ class DREMTubesOpticalPhysics : public G4VPhysicsConstructor { G4OpBoundaryProcess* theBoundaryProcess; G4bool AbsorptionOn; + + G4bool fFullOptic; }; diff --git a/include/DREMTubesPhysicsList.hh b/include/DREMTubesPhysicsList.hh index d124a63..6921433 100644 --- a/include/DREMTubesPhysicsList.hh +++ b/include/DREMTubesPhysicsList.hh @@ -23,7 +23,7 @@ class DREMTubesPhysicsList : public G4VModularPhysicsList{ public: //Constructor // - DREMTubesPhysicsList(G4String); + DREMTubesPhysicsList(G4String, const G4bool FullOptic ); //De-constructor // virtual ~DREMTubesPhysicsList(); @@ -31,6 +31,11 @@ class DREMTubesPhysicsList : public G4VModularPhysicsList{ DREMTubesOpticalPhysics* OpPhysics; G4bool AbsorptionOn; + + private: + + G4bool fFullOptic; + }; #endif diff --git a/src/DREMTubesDetectorConstruction.cc b/src/DREMTubesDetectorConstruction.cc index 70bf0a5..820e98e 100644 --- a/src/DREMTubesDetectorConstruction.cc +++ b/src/DREMTubesDetectorConstruction.cc @@ -317,7 +317,7 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { G4int NofFibersrow = 3*16; G4int NofFiberscolumn = 60; G4double moduleZ = (1000.)*mm; - double tolerance = 0.05*mm; + double tolerance = 0.0*mm; G4double moduleX = 3.*32.*mm+1.*mm+2.*tolerance*NofFibersrow; G4double moduleY = 59.*(1.733+2*tolerance)*mm+2.0*mm; @@ -453,8 +453,8 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { // Calorimeter placement (with rotation wrt beam axis) // G4RotationMatrix rotm = G4RotationMatrix(); - rotm.rotateY(1.0*deg); - rotm.rotateX(1.0*deg); + rotm.rotateY(0.0*deg); + rotm.rotateX(0.0*deg); G4ThreeVector position; position.setX(0.); position.setY(0.); diff --git a/src/DREMTubesOpticalPhysics.cc b/src/DREMTubesOpticalPhysics.cc index f33528d..e8e3665 100644 --- a/src/DREMTubesOpticalPhysics.cc +++ b/src/DREMTubesOpticalPhysics.cc @@ -18,8 +18,9 @@ // Constructor // -DREMTubesOpticalPhysics::DREMTubesOpticalPhysics(G4bool toggle) - : G4VPhysicsConstructor("Optical") { +DREMTubesOpticalPhysics::DREMTubesOpticalPhysics(const G4bool FullOptic, G4bool toggle) + : G4VPhysicsConstructor("Optical"), + fFullOptic( FullOptic ) { // Initialize private members // @@ -126,16 +127,17 @@ void DREMTubesOpticalPhysics::ConstructProcess() { pManager->SetProcessOrdering(theCerenkovProcess,idxPostStep); } // Add Scintillation process to each candidate + // Adding Scintillation process only if fFullOptic == true // if(theScintProcess->IsApplicable(*particle)){ - // If uncommented simulations gets reeeeeeeeeeally time consuming - // (your time-life could not be sufficient) - // - //pManager->AddProcess(theScintProcess); - //pManager->SetProcessOrderingToLast(theScintProcess,idxAtRest); - //pManager->SetProcessOrderingToLast(theScintProcess,idxPostStep); + if (fFullOptic) { + pManager->AddProcess(theScintProcess); + pManager->SetProcessOrderingToLast(theScintProcess,idxAtRest); + pManager->SetProcessOrderingToLast(theScintProcess,idxPostStep); + } + else {} } - } + }//end while } //************************************************** diff --git a/src/DREMTubesPhysicsList.cc b/src/DREMTubesPhysicsList.cc index c9e274c..18eba23 100644 --- a/src/DREMTubesPhysicsList.cc +++ b/src/DREMTubesPhysicsList.cc @@ -23,7 +23,9 @@ //Define constructor // -DREMTubesPhysicsList::DREMTubesPhysicsList(G4String physName):G4VModularPhysicsList(){ +DREMTubesPhysicsList::DREMTubesPhysicsList(G4String physName, const G4bool FullOptic ) + :G4VModularPhysicsList(), + fFullOptic( FullOptic ) { //Define physics list factor and modular physics list // @@ -58,7 +60,7 @@ DREMTubesPhysicsList::DREMTubesPhysicsList(G4String physName):G4VModularPhysicsL // Turn on and off the absorption of optical photons in materials // AbsorptionOn = true; - RegisterPhysics( OpPhysics = new DREMTubesOpticalPhysics(AbsorptionOn) ); + RegisterPhysics( OpPhysics = new DREMTubesOpticalPhysics(fFullOptic, AbsorptionOn) ); } From 6677a8c68b6df42d02fe887acf8e3ed81f561227 Mon Sep 17 00:00:00 2001 From: Lorenzo Pezzotti Date: Tue, 20 Jul 2021 10:59:02 +0200 Subject: [PATCH 5/6] adding PDGID --- include/DREMTubesEventAction.hh | 8 ++++---- src/DREMTubesEventAction.cc | 3 ++- src/DREMTubesRunAction.cc | 2 +- src/DREMTubesSteppingAction.cc | 6 ++++-- 4 files changed, 11 insertions(+), 8 deletions(-) diff --git a/include/DREMTubesEventAction.hh b/include/DREMTubesEventAction.hh index 72039b9..607434c 100644 --- a/include/DREMTubesEventAction.hh +++ b/include/DREMTubesEventAction.hh @@ -41,7 +41,7 @@ class DREMTubesEventAction : public G4UserEventAction { void Addenergy(G4double de);//Add all energy deposited //void AddEnergyfibre(G4double de, G4int number);//Add energy in fiber cpn //void AddSignalfibre(G4int number); - void SavePrimaryParticle(G4String name); + void SavePrimaryPDGID(G4int pdgid); void SaveAbsorberMaterial(G4String AbsorberMaterialName); void SavePrimaryEnergy(G4double primaryparticleenergy); void AddEscapedEnergy(G4double escapedenergy); @@ -66,7 +66,7 @@ class DREMTubesEventAction : public G4UserEventAction { G4double EnergyTot;//Total energy deposited (does not count invisibile energy) //G4double Signalfibre[64]; ////Signal in 64 single module fibers, to be used with AddEnergyfibre - G4String PrimaryParticleName; //Name of primary particle + G4int PrimaryPDGID; //PDGID of primary particle G4String AbsorberMaterial; //Name of absorber material G4double PrimaryParticleEnergy;//Primary particle energy G4double EscapedEnergy; @@ -84,8 +84,8 @@ inline void DREMTubesEventAction::AddEscapedEnergy(G4double escapedenergy){ EscapedEnergy += escapedenergy; } -inline void DREMTubesEventAction::SavePrimaryParticle(G4String name){ - PrimaryParticleName = name; +inline void DREMTubesEventAction::SavePrimaryPDGID(G4int pdgid){ + PrimaryPDGID = pdgid; } inline void DREMTubesEventAction::SaveAbsorberMaterial(G4String AbsorberMaterialName){ diff --git a/src/DREMTubesEventAction.cc b/src/DREMTubesEventAction.cc index e1bacf4..d989b24 100644 --- a/src/DREMTubesEventAction.cc +++ b/src/DREMTubesEventAction.cc @@ -34,6 +34,7 @@ DREMTubesEventAction::DREMTubesEventAction() NofCherenkovDetected(0), //NofScintillationDetected(0), EnergyTot(0.), + PrimaryPDGID(0), PrimaryParticleEnergy(0.), EscapedEnergy(0.), VectorSignals(0.), @@ -96,7 +97,7 @@ void DREMTubesEventAction::EndOfEventAction(const G4Event* ) { analysisManager->FillNtupleDColumn(3, NofCherenkovDetected); analysisManager->FillNtupleDColumn(4, EnergyTot); analysisManager->FillNtupleDColumn(5, PrimaryParticleEnergy); - analysisManager->FillNtupleSColumn(6, PrimaryParticleName); + analysisManager->FillNtupleIColumn(6, PrimaryPDGID); //analysisManager->FillNtupleSColumn(7, AbsorberMaterial); analysisManager->FillNtupleDColumn(7, EscapedEnergy); //analysisManager->FillNtupleDColumn(9, Energyem2); diff --git a/src/DREMTubesRunAction.cc b/src/DREMTubesRunAction.cc index 23ab9d4..5587c3b 100644 --- a/src/DREMTubesRunAction.cc +++ b/src/DREMTubesRunAction.cc @@ -51,7 +51,7 @@ DREMTubesRunAction::DREMTubesRunAction( DREMTubesEventAction* eventAction ) analysisManager->CreateNtupleDColumn("NofCherenkovDetected"); analysisManager->CreateNtupleDColumn("EnergyTot"); analysisManager->CreateNtupleDColumn("PrimaryParticleEnergy"); - analysisManager->CreateNtupleSColumn("PrimaryParticleName"); + analysisManager->CreateNtupleIColumn("PrimaryPDGID"); analysisManager->CreateNtupleDColumn("EscapedEnergy"); analysisManager->CreateNtupleDColumn ("VectorSignals", fEventAction->GetVectorSignals()); diff --git a/src/DREMTubesSteppingAction.cc b/src/DREMTubesSteppingAction.cc index e6ccfdb..04c930e 100644 --- a/src/DREMTubesSteppingAction.cc +++ b/src/DREMTubesSteppingAction.cc @@ -61,6 +61,8 @@ void DREMTubesSteppingAction::UserSteppingAction( const G4Step* step ) { //Define SlowSteppingAction() method // void DREMTubesSteppingAction::SlowSteppingAction( const G4Step* step ){ + + //step->GetTrack()->SetTrackStatus(fStopAndKill); } @@ -76,7 +78,7 @@ void DREMTubesSteppingAction::AuxSteppingAction( const G4Step* step ) { //= step->GetPostStepPoint()->GetTouchableHandle()->GetVolume(); G4double energydeposited = step->GetTotalEnergyDeposit(); G4String particlename = step->GetTrack()->GetDefinition()->GetParticleName(); - + G4int particlepdg = step->GetTrack()->GetDefinition()->GetPDGEncoding(); //-------------------------------------------------- //Store auxiliary information from event steps //-------------------------------------------------- @@ -102,7 +104,7 @@ void DREMTubesSteppingAction::AuxSteppingAction( const G4Step* step ) { if ( step->GetTrack()->GetTrackID() == 1 && step->GetTrack()->GetCurrentStepNumber() == 1){ //Save primary particle energy and name - fEventAction->SavePrimaryParticle(particlename); + fEventAction->SavePrimaryPDGID(particlepdg); fEventAction->SavePrimaryEnergy(step->GetTrack()->GetKineticEnergy()); } From 23bd465821836a771e093500d81e5a53b5373fad Mon Sep 17 00:00:00 2001 From: Lorenzo Pezzotti Date: Tue, 20 Jul 2021 19:11:36 +0200 Subject: [PATCH 6/6] first SlowSteppingAction implementation --- include/DREMTubesDetectorConstruction.hh | 2 +- include/DREMTubesRunAction.hh | 2 +- src/DREMTubesDetectorConstruction.cc | 44 +++++++----- src/DREMTubesEventAction.cc | 12 ++++ src/DREMTubesSteppingAction.cc | 90 +++++++++++++++++++++++- 5 files changed, 128 insertions(+), 22 deletions(-) diff --git a/include/DREMTubesDetectorConstruction.hh b/include/DREMTubesDetectorConstruction.hh index 26607ed..9fe0e8f 100644 --- a/include/DREMTubesDetectorConstruction.hh +++ b/include/DREMTubesDetectorConstruction.hh @@ -8,7 +8,7 @@ //Prevent including header multplie times // #ifndef DREMTubesDetectorConstruction_h -#define DREMTUbesDetectorConstruction_h 1 +#define DREMTubesDetectorConstruction_h 1 //Includers from Geant4 // diff --git a/include/DREMTubesRunAction.hh b/include/DREMTubesRunAction.hh index 67f1747..fdf4838 100644 --- a/include/DREMTubesRunAction.hh +++ b/include/DREMTubesRunAction.hh @@ -1,6 +1,6 @@ //************************************************** // \file DREMTubesRunAction.hh -// \brief: Definition of DREMTUbesBRunAction class +// \brief: Definition of DREMTubesBRunAction class // \author: Lorenzo Pezzotti (CERN EP-SFT-sim) @lopezzot // \start date: 7 July 2021 //************************************************** diff --git a/src/DREMTubesDetectorConstruction.cc b/src/DREMTubesDetectorConstruction.cc index 820e98e..465af60 100644 --- a/src/DREMTubesDetectorConstruction.cc +++ b/src/DREMTubesDetectorConstruction.cc @@ -496,15 +496,25 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { OpSurfaceGlassSi -> SetType(dielectric_metal); OpSurfaceGlassSi -> SetModel(glisur); OpSurfaceGlassSi -> SetFinish(polished); - G4double efficiencyOpSurfaceGlassSi[ENTRIES] = // detection efficiency - { 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4, - 0.4, 0.4, 0.4, 0.4}; + G4double efficiencyOpSurfaceGlassSi[ENTRIES] = //100% detection efficiency + { 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1, + 1, 1, 1, 1 }; + + /*G4double efficiencyOpSurfaceGlassSi[ENTRIES] = //0% detection efficiency + { 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0, + 0, 0, 0, 0 };*/ G4double reflectivityOpSurfaceGlassSi[ENTRIES] = // 0% reflection { 0., 0., 0., 0., @@ -559,8 +569,8 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { // Logical Skin Surface placement around the silicon of the SiPM // - /*G4LogicalSkinSurface* OpsurfaceSi = new G4LogicalSkinSurface("OpsurfaceSi", - SiLV, OpSurfaceGlassSi);*/ + /*G4LogicalSkinSurface* OpsurfaceSi =*/ new G4LogicalSkinSurface("OpsurfaceSi", + SiLV, OpSurfaceGlassSi); // Optical Surface properties between the scintillating fibers // and the default material @@ -598,8 +608,8 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { std::string S_name; std::string SiPM_name; S_name = "S_row" + S_fiber_row.str() + "_column_" + S_fiber_column.str(); - SiPM_name = - "SiPMS_row" + S_fiber_row.str() + "_column_" + S_fiber_column.str(); + SiPM_name = "S_SiPM"; + //SiPM_name = "SiPMS_row" + S_fiber_row.str() + "_column_" + S_fiber_column.str(); G4double S_x, S_y; G4ThreeVector vec_S_fiber; @@ -648,7 +658,7 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { SiPM_name, moduleequippedLV, false, - 0); + copynumber); //same copynumber of fibers /*logic_OpSurface_defaultAir[NofFibersrow][NofFiberscolumn] = new G4LogicalBorderSurface("logic_OpSurface_defaultAir", @@ -676,8 +686,8 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { std::string C_name; std::string SiPM_name; C_name = "C_row" + C_fiber_row.str() + "_column_" + C_fiber_column.str(); - SiPM_name = - "SiPMC_row" + C_fiber_row.str() + "_column_" + C_fiber_column.str(); + SiPM_name = "C_SiPM"; + //SiPM_name = "SiPMC_row" + C_fiber_row.str() + "_column_" + C_fiber_column.str(); G4double C_x, C_y; G4ThreeVector vec_C_fiber; @@ -723,7 +733,7 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { SiPM_name, moduleequippedLV, false, - 0); + copynumber); //same copynumber of fiber /*logic_OpSurface_defaultAir[NofFibersrow][NofFiberscolumn] = new G4LogicalBorderSurface("logic_OpSurface_defaultAir", diff --git a/src/DREMTubesEventAction.cc b/src/DREMTubesEventAction.cc index d989b24..9cfe321 100644 --- a/src/DREMTubesEventAction.cc +++ b/src/DREMTubesEventAction.cc @@ -103,6 +103,18 @@ void DREMTubesEventAction::EndOfEventAction(const G4Event* ) { //analysisManager->FillNtupleDColumn(9, Energyem2); analysisManager->AddNtupleRow(); + G4int tot_S = 0; + G4int tot_C = 0; + for(unsigned int i=0; iGetTrack()->SetTrackStatus(fStopAndKill); + + //Random seed and random number generator + // + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::default_random_engine generator(seed); + + std::poisson_distribution scin_distribution(0.008); + std::poisson_distribution cher_distribution(0.25); + + G4double energydeposited = step->GetTotalEnergyDeposit(); + G4VPhysicalVolume* PreStepVolume + = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume(); + + std::string Fiber; + std::string S_fiber = "S_fiber"; + std::string C_fiber = "C_fiber"; + Fiber = PreStepVolume->GetName(); //name of current step fiber + + if ( strstr( Fiber.c_str(), S_fiber.c_str() ) ) { //scintillating fiber + fEventAction->AddScin(energydeposited); + } + + if ( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fiber + fEventAction->AddCher(energydeposited); + } + + G4String particlename = step->GetTrack()->GetDefinition()->GetParticleName(); + + G4OpBoundaryProcessStatus theStatus = Undefined; + + G4ProcessManager* OpManager = + G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); + + if (OpManager) { + G4int MAXofPostStepLoops = + OpManager->GetPostStepProcessVector()->entries(); + G4ProcessVector* fPostStepDoItVector = + OpManager->GetPostStepProcessVector(typeDoIt); + + for ( G4int i=0; i(fCurrentProcess); + if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;} + } + } + + if( particlename == "opticalphoton" ) { //optical photons + switch ( theStatus ){ + + case TotalInternalReflection: + if(strstr(Fiber.c_str(),S_fiber.c_str())){ //scintillating fibre + if (scin_distribution(generator)==0 && step->GetTrack()->GetCurrentStepNumber() == 1){ + step->GetTrack()->SetTrackStatus(fStopAndKill); + } + } + else if( strstr( Fiber.c_str(), C_fiber.c_str() ) ) { //Cherenkov fibre + if (cher_distribution(generator)==0 && step->GetTrack()->GetCurrentStepNumber() == 1){ + step->GetTrack()->SetTrackStatus(fStopAndKill); + } + } + + case Detection: + + if (step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() == "C_SiPM"){ + fEventAction->AddCherenkov(1); + fEventAction->AddVectorCherPE(step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(), 1); + } + if (step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() == "S_SiPM"){ + fEventAction->AddVectorScinEnergy(1, step->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber()); + } + + //Print info on Detection position and volumes + // + /*G4cout<<"Detection "<GetPreStepPoint()->GetMaterial()->GetName()<< + " "<GetPostStepPoint()->GetMaterial()->GetName()<< + " "<GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName()<< + " "<GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber()<GetTotalEnergyDeposit(); G4String particlename = step->GetTrack()->GetDefinition()->GetParticleName(); G4int particlepdg = step->GetTrack()->GetDefinition()->GetPDGEncoding(); + //-------------------------------------------------- //Store auxiliary information from event steps //-------------------------------------------------- @@ -182,7 +266,7 @@ void DREMTubesSteppingAction::FastSteppingAction( const G4Step* step ) { for ( G4int i=0; i(fCurrentProcess); - if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break;} + if (fOpProcess) { theStatus = fOpProcess->GetStatus(); break; } } }