Skip to content
Snippets Groups Projects
Commit fbd6e4c3 authored by MANZANILLAS Luis's avatar MANZANILLAS Luis
Browse files

first commit

parent f25c39e9
Branches
No related tags found
No related merge requests found
Ge2SOS.cc 0 → 100644
#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4SteppingVerbose.hh"
#include "Randomize.hh"
#include "G4MTRunManager.hh"
#include "DetectorConstruction.hh"
#include "PhysicsList.hh"
#include "ActionInitialization.hh"
#include "G4VisExecutive.hh"
#include "G4UIExecutive.hh"
#include "unistd.h"
#include "time.h"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace {
void PrintUsage() {
G4cerr << " Usage: " << G4endl;
G4cerr << " OpNovice [-m macro ] [-u UIsession] [-t nThreads] [-r seed] "
<< G4endl;
G4cerr << " note: -t option is available only for multi-threaded mode."
<< G4endl;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
int main(int argc,char** argv)
{
G4cout << " " << G4endl
<< "======================================================="<< G4endl
<< " "<< G4endl
<< " ###### ###### ###### ###### "<< G4endl
<< " # # # # # "<< G4endl
<< " # # # # # "<< G4endl
<< " # ###### #### ####### # # ####### "<< G4endl
<< " ## # # # # # # # "<< G4endl
<< " # # # # # # # # "<< G4endl
<< " ######## #### ####### ###### ####### "<< G4endl
<< " "<< G4endl
<< "======================================================="<< G4endl
<< " Germanium detectors at SOLEIL "<< G4endl
<< " Paco Iguaz Gutierrez, Luis Manzanillas Velez "<< G4endl
<< "======================================================="<< G4endl
<< G4endl << G4endl;
// Evaluate arguments
//
G4UIExecutive* ui = nullptr;
if (argc == 1) ui = new G4UIExecutive(argc,argv);
G4String macro;
G4String session;
//G4long myseed = 345354;
G4long myseed = time((time_t *)NULL);
for ( G4int i=1; i<argc; i=i+2 ) {
if ( G4String(argv[i]) == "-m" ) macro = argv[i+1];
else if ( G4String(argv[i]) == "-u" ) session = argv[i+1];
else if ( G4String(argv[i]) == "-r" ) myseed = atoi(argv[i+1]);
else {
PrintUsage();
return 1;
}
}
// Choose the Random engine
//
G4Random::setTheEngine(new CLHEP::RanecuEngine);
// Construct the default run manager
//
G4RunManager * runManager = new G4RunManager;
// Seed the random number generator manually
G4Random::setTheSeed(myseed);
// Initialize G4 kernel
//
// Set mandatory initialization classes
//
// Detector construction
DetectorConstruction* det = new DetectorConstruction;
runManager->SetUserInitialization(det);
//
// Physics list
runManager-> SetUserInitialization(new PhysicsList());
// User action initialization
runManager->SetUserInitialization(new ActionInitialization(det));
runManager->Initialize();
//initialize visualization
G4VisManager* visManager = nullptr;
//get the pointer to the User Interface manager
G4UImanager* UImanager = G4UImanager::GetUIpointer();
if (ui) {
//interactive mode
visManager = new G4VisExecutive;
visManager->Initialize();
UImanager->ApplyCommand("/control/execute vis.mac");
ui->SessionStart();
delete ui;
} else {
//batch mode
G4String command = "/control/execute ";
UImanager->ApplyCommand(command+macro);
}
//job termination
delete visManager;
delete runManager;
return 0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
gui.mac 0 → 100644
/process/em/fluo true
/process/em/auger true
/process/em/augerCascade true
/process/em/pixe true
#/Ge2SOS/det/setOutputDirectory /home/iwsatlas1/manzanil/Documents/Varios/SOLEIL/
/Ge2SOS/det/setDetectorType 1
/Ge2SOS/det/setNTargetSamples 1
/Ge2SOS/det/setGeDetectorThickness 2. mm
#Set the position of the collimator, a values between -37 and 37 mm (size of samples)
#Choise of source type: 0 Cs-137; 1 Bi-207; 2 Sr-90; 3 e-
/run/initialize
/Ge2SOS/gun/sourceType 0
/run/beamOn 1
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// $Id: ActionInitialization.cc 98759 2016-08-09 14:03:09Z gcosmo $
//
/// \file ActionInitialization.cc
/// \brief Implementation of the ActionInitialization class
#include "ActionInitialization.hh"
#include "DetectorConstruction.hh"
#include "PrimaryGeneratorAction.hh"
#include "RunAction.hh"
#include "EventAction.hh"
#include "SteppingAction.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ActionInitialization::ActionInitialization(DetectorConstruction* det)
//ActionInitialization::ActionInitialization()
: G4VUserActionInitialization()
{
fDetector = det;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ActionInitialization::~ActionInitialization()
{
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ActionInitialization::BuildForMaster() const
{
//PrimaryGeneratorAction* prim = new PrimaryGeneratorAction(fDetector);
//SetUserAction(new RunAction(fDetector, prim));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ActionInitialization::Build() const
{
PrimaryGeneratorAction* prim = new PrimaryGeneratorAction(fDetector);
SetUserAction(prim);
//G4cout << "Prim"<<G4endl;
RunAction* run = new RunAction(fDetector,prim);
SetUserAction(run);
//G4cout << "Run"<<G4endl;
EventAction* event = new EventAction(fDetector,run);
SetUserAction(event);
//G4cout<<"Event"<<G4endl;
SetUserAction(new SteppingAction(event));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
This diff is collapsed.
#include "DetectorConstruction.hh"
#include "DetectorMessenger.hh"
#include "SiliconPlateConstruction.hh"
#include "CollimatorConstruction.hh"
#include "PrimaryGeneratorAction.hh"
#include "G4MaterialTable.hh"
#include "SoleilMaterials.hh"
#include "G4Material.hh"
#include "G4Element.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4OpticalSurface.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4Sphere.hh"
#include "G4Cons.hh"
#include "G4Torus.hh"
#include "G4Hype.hh"
#include "G4Transform3D.hh"
#include "G4LogicalVolume.hh"
#include "G4ThreeVector.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4UnitsTable.hh"
#include "G4NistManager.hh"
#include <G4UserLimits.hh>
#include "G4MultiFunctionalDetector.hh"
#include "G4SDManager.hh"
#include "G4PSEnergyDeposit.hh"
#include <G4VPrimitiveScorer.hh>
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4UnionSolid.hh"
#include "G4VoxelLimits.hh"
#include "G4MTRunManager.hh"
#include "G4PhysicalConstants.hh"
#include "G4GeometryManager.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4SolidStore.hh"
#include "G4Navigator.hh"
#include "G4TransportationManager.hh"
#include "G4GDMLParser.hh"
#include <G4VisAttributes.hh>
#include <iostream>
#include <fstream>
#include <iterator>
using namespace std;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
/*
Constructs DetectorConstruction, defines default values.
*/
DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(), geDetectorTube(nullptr),
logicCollimatorLEAPS(nullptr),
logicSample(nullptr),
geLogicVolume(nullptr)
{
fDetectorMessenger = new DetectorMessenger(this);
halfSizeDarkBoxX = halfSizeDarkBoxY = halfSizeDarkBoxZ = 1.*m;
fSetupName = "Polar";
fDataType = "csv";
detectorInnerRadius = 200.0*mm;
detectorThickness = 10.0*mm;
detectorHeight = 11.5*mm;
sampleInnerRadius = 0.0*mm;
sampleThickness = 1.5*mm;
sampleHeight = 5.0*mm;
beWindowRadius = 8.0*mm;
CollimatorThickness = 20.0*mm;
distanceCollimatorDetector = 5.0*mm;
ContactThickness = 300.0*nm;
nSamples = 1;
fDetectorType = 0;
fDetectorName = "Polarimeter";
fVolName = "World";
materialConstruction = new SoleilMaterials;
DefineMaterials();
data_output_directory = "./";
fCollimatorMaterial = G4Material::GetMaterial("titanium_grade1");
fTargetMaterial = G4Material::GetMaterial("G4_Ge");
fSampleMaterial = G4Material::GetMaterial("G4_POLYETHYLENE");
fWorldMaterial = G4Material::GetMaterial("G4_Galactic");
materialGeContainer = G4Material::GetMaterial("Aluminium");
materialGeContainerCoating = G4Material::GetMaterial("Aluminium");
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
DetectorConstruction::~DetectorConstruction(){
//delete physicWorldBox;
delete fDetectorMessenger;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4VPhysicalVolume* DetectorConstruction::GetPhysicalVolumeByName(const G4String &name)
{
// access the store of physical volumes
G4PhysicalVolumeStore * pvs= G4PhysicalVolumeStore::GetInstance();
G4VPhysicalVolume* pv;
G4int npv= pvs->size();
G4int ipv;
for (ipv=0; ipv<npv; ipv++) {
pv= (*pvs)[ipv];
if (!pv)
break;
//G4cout<<" pv->GetName() "<<pv->GetName()<<G4endl;
if (pv->GetName() == name)
return pv;
}
return NULL;
}
/*
Sets which detector geometry is used.
*/
void DetectorConstruction::SetDetectorType(G4int value){
fDetectorType=value;
G4RunManager::GetRunManager()->ReinitializeGeometry();
//G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
void DetectorConstruction::SetNumberOfTargetSamples(G4int value){
nSamples = value;
G4RunManager::GetRunManager()->ReinitializeGeometry();
//UpdateGeometry();
}
//Sets dimmensions of target, thickness corresponds to the Z coordinate, Length to x.
void DetectorConstruction::SetGeDetectorLength(G4double value){
detectorInnerRadius = (value/2.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
void DetectorConstruction::SetGeDetectorThickness(G4double value){
detectorThickness = (value/2.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
void DetectorConstruction::SetGeDetectorWidth(G4double value){
detectorHeight = (value/2.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
void DetectorConstruction::SetContactThickness(G4double value){
ContactThickness = (value/1.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
//Set sample dimmensions
void DetectorConstruction::SetSampleLength(G4double value){
sampleInnerRadius = (value/2.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
void DetectorConstruction::SetSampleThickness(G4double value){
sampleThickness = (value/2.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
void DetectorConstruction::SetSampleWidth(G4double value){
sampleHeight = (value/2.)*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
void DetectorConstruction::SetBeWindowRadius(G4double value){
beWindowRadius = value*mm;
//UpdateGeometry();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
//Collimator thickness
void DetectorConstruction::SetCollimatorThickness(G4double value){
CollimatorThickness = value*mm;
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
//Collimator to detector distance
void DetectorConstruction::SetDistanceCollimatorDetector(G4double value){
distanceCollimatorDetector = value*mm;
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
//Sets material of target.
void DetectorConstruction::SetTargetMaterial(G4String materialChoice)
{
// search the material by its name
G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
if (pttoMaterial) {
fTargetMaterial = pttoMaterial;
if(geLogicVolume)geLogicVolume->SetMaterial(fTargetMaterial);
G4cout<<" material: "<<fTargetMaterial->GetName()<<G4endl;
} else {
G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found" << G4endl;
}
G4RunManager::GetRunManager()->ReinitializeGeometry();
G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
//Sets material of Ge Container
void DetectorConstruction::SetGeContainerMaterial(G4String materialChoice)
{
// search the material by its name
G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
if (pttoMaterial) {
materialGeContainer = pttoMaterial;
if(logicGeContainer)logicGeContainer->SetMaterial(materialGeContainer);
G4cout<<" material container: "<<materialGeContainer->GetName()<<G4endl;
} else {
G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found" << G4endl;
}
G4RunManager::GetRunManager()->ReinitializeGeometry();
G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
//Sets material of Ge Container coating
void DetectorConstruction::SetGeContainerMaterialCoating(G4String materialChoice)
{
// search the material by its name
G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
if (pttoMaterial) {
materialGeContainerCoating = pttoMaterial;
if(logicGeContainerCoating)logicGeContainerCoating->SetMaterial(materialGeContainerCoating);
G4cout<<" material container coating: "<<materialGeContainerCoating->GetName()<<G4endl;
} else {
G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found" << G4endl;
}
G4RunManager::GetRunManager()->ReinitializeGeometry();
G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
void DetectorConstruction::SetDetectorName(G4String detectorNAME)
{
fDetectorName = detectorNAME;
}
void DetectorConstruction::SetSetupName(G4String setupNAME)
{
fSetupName = setupNAME;
}
void DetectorConstruction::SetDataType(G4String dataType)
{
fDataType = dataType;
}
/*
Sets material of sample.
*/
void DetectorConstruction::SetSampleMaterial(G4String materialChoice)
{
// search the material by its name
G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
if (pttoMaterial) {
fSampleMaterial = pttoMaterial;
if(logicSample)logicSample->SetMaterial(fSampleMaterial);
G4cout<<" material "<<fSampleMaterial->GetName()<<G4endl;
} else {
G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found" << G4endl;
}
G4RunManager::GetRunManager()->ReinitializeGeometry();
G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
void DetectorConstruction::SetOutputDirectory(G4String output_directory)
{
data_output_directory = output_directory;
}
/*
Sets material of world volume.
*/
void DetectorConstruction::SetWorldMaterial(G4String materialChoice)
{
// search the material by its name
G4Material* pttoMaterial =
G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);
if (pttoMaterial) {
fWorldMaterial = pttoMaterial;
if ( logicWorldBox ) { logicWorldBox->SetMaterial(fWorldMaterial); }
} else {
G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found" << G4endl;
}
G4RunManager::GetRunManager()->ReinitializeGeometry();
//G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
void DetectorConstruction::SetCollimatorMaterial(G4String materialChoice)
{
// search the material by its name
G4Material* pttoMaterial =
G4NistManager::Instance()->FindOrBuildMaterial(materialChoice);
if (pttoMaterial) {
fCollimatorMaterial = pttoMaterial;
if ( logicCollimatorLEAPS ) { logicCollimatorLEAPS->SetMaterial(fCollimatorMaterial); }
} else {
G4cout << "\n--> warning from DetectorConstruction::SetMaterial : "
<< materialChoice << " not found" << G4endl;
}
G4RunManager::GetRunManager()->ReinitializeGeometry();
//G4MTRunManager::GetRunManager()->PhysicsHasBeenModified();
}
/*
Defines materials used in simulation. Sets material properties for PEN and other optical components.
*/
void DetectorConstruction::DefineMaterials(){
// ============================================================= Materials =============================================================
//materialConstruction = new PenMaterials;
materialConstruction-> Construct();
materialAir = G4Material::GetMaterial("Air");
materialBialkali = G4Material::GetMaterial("Bialkali");
fGlass = G4Material::GetMaterial("BorosilicateGlass");
PenMaterial = G4Material::GetMaterial("PEN");
PVTMaterial = G4Material::GetMaterial("PVT_structure");
materialSi = G4Material::GetMaterial("G4_Si");
materialGe = G4Material::GetMaterial("G4_Ge");
materialTriggerFoilEJ212 = G4Material::GetMaterial("EJ212");
Pstyrene = G4Material::GetMaterial("Polystyrene");
materialPMMA = G4Material::GetMaterial("PMMA");
fVacuum = G4Material::GetMaterial("G4_Galactic");
materialGreaseEJ550 = G4Material::GetMaterial("Grease");
materialTeflon = G4Material::GetMaterial("G4_TEFLON");
materialVikuiti = G4Material::GetMaterial("Vikuiti");
materialTitanium = G4Material::GetMaterial("titanium_grade1");
materialCollimatorCoating = G4Material::GetMaterial("Aluminium");
materialSample = G4Material::GetMaterial("Iron");
materialBeWindow = G4Material::GetMaterial("G4_Be");
materialMetalTube = G4Material::GetMaterial("G4_STAINLESS-STEEL");
materialSupportPins = G4Material::GetMaterial("PEEK");
materialContacts = G4Material::GetMaterial("Aluminium");
materialPolyethylene = G4Material::GetMaterial("G4_POLYETHYLENE");
G4cout<<" materials imported succesfully "<<G4endl;
}
void DetectorConstruction::SetVolName(G4ThreeVector thePoint){
G4Navigator* theNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
G4VPhysicalVolume* myVolume= theNavigator->LocateGlobalPointAndSetup(thePoint);
fVolName = myVolume->GetName();
}
void DetectorConstruction::UpdateGeometry(){
//define new one
G4GeometryManager::GetInstance()->OpenGeometry();
G4PhysicalVolumeStore::Clean();
G4LogicalVolumeStore::Clean();
G4SolidStore::Clean();
G4RunManager::GetRunManager()->ReinitializeGeometry();
}
/*
Clears stored geometry, then constructs all volumes that can be used in the simulation.
Builds and places volumes in world.
Defines detector sensitivities and properties.
*/
G4VPhysicalVolume* DetectorConstruction::Construct()
{
// ============================================================= Define Volumes =============================================================
G4GeometryManager::GetInstance()->OpenGeometry();
G4PhysicalVolumeStore::GetInstance()->Clean();
G4LogicalVolumeStore::GetInstance()->Clean();
G4SolidStore::GetInstance()->Clean();
G4SolidStore::GetInstance()->Clean();
G4LogicalSkinSurface::CleanSurfaceTable();
G4LogicalBorderSurface::CleanSurfaceTable();
//The experimental Dark Box
fWorldBox = new G4Box("World",halfSizeDarkBoxX,halfSizeDarkBoxY,halfSizeDarkBoxZ);
logicWorldBox = new G4LogicalVolume(fWorldBox,materialAir,"World",0,0,0);
physicWorldBox = new G4PVPlacement(0,G4ThreeVector(),logicWorldBox,"World",0,false,0);
G4double detectorExternalRadius = detectorInnerRadius + detectorThickness;
geDetectorTube = new G4Tubs("target", detectorInnerRadius, detectorExternalRadius, detectorHeight, 0.*deg, 360.*deg);
geLogicVolume = new G4LogicalVolume(geDetectorTube,fTargetMaterial, "target",0,0,0);
//solidGeContainer = new G4Tubs("ge_container", geContainerInnerRadius, geContainerOuterRadius, geContainerHalfThickness, 0.*deg,360.*deg);
G4double sampleExternalRadius = sampleInnerRadius + sampleThickness;
sampleTube = new G4Tubs("sample", sampleInnerRadius, sampleExternalRadius, sampleHeight, 0.*deg, 360.*deg);
logicSample = new G4LogicalVolume(sampleTube,fSampleMaterial, "sample",0,0,0);
//Detector
G4VisAttributes* visualAttributesGeDetector = new G4VisAttributes(G4Colour::Blue());
visualAttributesGeDetector->SetVisibility(true);
visualAttributesGeDetector->SetForceWireframe(true);
visualAttributesGeDetector->SetForceAuxEdgeVisible(true);
visualAttributesGeDetector->SetForceSolid(true);
geLogicVolume->SetVisAttributes(visualAttributesGeDetector);
G4VisAttributes* visualAttributesSample = new G4VisAttributes(G4Colour::Black());
visualAttributesSample->SetVisibility(true);
visualAttributesSample->SetForceSolid(true);
logicSample->SetVisAttributes(visualAttributesSample);
// ============================================================= Place volumes =============================================================
// Place main detector always at center of world volume
switch(fDetectorType){
//Detector alone for "direct" beam radiation
case 0:
new G4PVPlacement(0,
G4ThreeVector(0,0,0),
geLogicVolume,
"target_"+std::to_string(1),
logicWorldBox,false,1,false);
new G4PVPlacement(0,
G4ThreeVector(0,0,0),
logicSample,
"sample"+std::to_string(1),
logicWorldBox,false,1,false);
break;
//Detector + collimator for "direct" beam radiation
}
//// Definition of simulation steps.
//logicWorldBox->SetUserLimits(new G4UserLimits(0.1*mm));
//geLogicVolume->SetUserLimits(new G4UserLimits(1.0*um));
//logicBeWindow->SetUserLimits(new G4UserLimits(50.0*um));
//logicContactVolume->SetUserLimits(new G4UserLimits(5.0*nm));
return physicWorldBox;
}
#include "DetectorMessenger.hh"
#include "DetectorConstruction.hh"
#include "G4UIdirectory.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWithADouble.hh"
#include "G4UIcmdWithoutParameter.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4UIcmdWithABool.hh"
#include "G4RunManager.hh"
DetectorMessenger::DetectorMessenger(DetectorConstruction * Det)
:G4UImessenger(),
fDetector(Det),
fGe2SOSDir(0),
fDetDir(0),
commandSetWorldMaterial(0),
commandSetDetectorType(0),
commandSetNumberOfTargetSamples(0),
commandSetGeDetectorLength(0),
commandSetGeDetectorThickness(0),
commandSetGeDetectorWidth(0),
commandSetContactThickness(0),
commandSetBeWindowRadius(0),
commandSetCollimatorThickness(0),
commandSetDistanceCollimatorDetector(0),
commandSetCollimatorMaterial(0),
commandSetTargetMaterial(0),
commandSetGeContainerMaterial(0),
commandSetGeContainerMaterialCoating(0),
commandSetDetectorName(0),
commandSetSetupName(0),
commandSetDataType(0),
commandSetOutputDirectory(0)
{
fDetDir = new G4UIdirectory("/Ge2SOS/det/");
fDetDir->SetGuidance("detector construction commands");
commandSetCollimatorMaterial = new G4UIcmdWithAString("/Ge2SOS/det/setCollimatorMat",this);
commandSetCollimatorMaterial->SetGuidance("Select material of the collimator.");
commandSetCollimatorMaterial->SetParameterName("choice",false);
commandSetCollimatorMaterial->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetCollimatorMaterial->SetToBeBroadcasted(false);
commandSetTargetMaterial = new G4UIcmdWithAString("/Ge2SOS/det/setGeDetectorMat",this);
commandSetTargetMaterial->SetGuidance("Select material of the target.");
commandSetTargetMaterial->SetParameterName("choice",false);
commandSetTargetMaterial->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetTargetMaterial->SetToBeBroadcasted(false);
commandSetGeContainerMaterial = new G4UIcmdWithAString("/Ge2SOS/det/setGeContainerMat",this);
commandSetGeContainerMaterial->SetGuidance("Select material of the ge container.");
commandSetGeContainerMaterial->SetParameterName("choice",false);
commandSetGeContainerMaterial->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetGeContainerMaterial->SetToBeBroadcasted(false);
commandSetGeContainerMaterialCoating = new G4UIcmdWithAString("/Ge2SOS/det/setGeContainerMatCoating",this);
commandSetGeContainerMaterialCoating->SetGuidance("Select material of the ge container coating.");
commandSetGeContainerMaterialCoating->SetParameterName("choice",false);
commandSetGeContainerMaterialCoating->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetGeContainerMaterialCoating->SetToBeBroadcasted(false);
commandSetDetectorName = new G4UIcmdWithAString("/Ge2SOS/det/setGeDetectorName",this);
commandSetDetectorName->SetGuidance("Select name of detector.");
commandSetDetectorName->SetParameterName("choice",false);
commandSetDetectorName->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetDetectorName->SetToBeBroadcasted(false);
commandSetSetupName = new G4UIcmdWithAString("/Ge2SOS/det/setSetupName",this);
commandSetSetupName->SetGuidance("Select name of setup.");
commandSetSetupName->SetParameterName("choice",false);
commandSetSetupName->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetSetupName->SetToBeBroadcasted(false);
commandSetDataType = new G4UIcmdWithAString("/Ge2SOS/det/setDataType",this);
commandSetDataType->SetGuidance("Select format of data: csv, hdf5, root.");
commandSetDataType->SetParameterName("choice",false);
commandSetDataType->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetDataType->SetToBeBroadcasted(false);
commandSetOutputDirectory = new G4UIcmdWithAString("/Ge2SOS/det/setOutputDirectory",this);
commandSetOutputDirectory->SetGuidance("Set output directory");
commandSetOutputDirectory->SetParameterName("choice",false);
commandSetOutputDirectory->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetOutputDirectory->SetToBeBroadcasted(false);
commandSetWorldMaterial = new G4UIcmdWithAString("/Ge2SOS/det/setWorldMat",this);
commandSetWorldMaterial->SetGuidance("Select material of the world.");
commandSetWorldMaterial->SetParameterName("choice",false);
commandSetWorldMaterial->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetWorldMaterial->SetToBeBroadcasted(false);
commandSetGeDetectorLength = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setGeDetectorLength",this);
commandSetGeDetectorLength->SetGuidance("Set length of target samples");
commandSetGeDetectorLength->SetParameterName("SampleLength",false);
commandSetGeDetectorLength->SetRange("SampleLength>0.");
commandSetGeDetectorLength->SetUnitCategory("Length");
commandSetGeDetectorLength->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetGeDetectorLength->SetToBeBroadcasted(false);
commandSetGeDetectorThickness = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setGeDetectorThickness",this);
commandSetGeDetectorThickness->SetGuidance("Set thickness of target samples");
commandSetGeDetectorThickness->SetParameterName("SampleThickness",false);
commandSetGeDetectorThickness->SetRange("SampleThickness>0.");
commandSetGeDetectorThickness->SetUnitCategory("Length");
commandSetGeDetectorThickness->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetGeDetectorThickness->SetToBeBroadcasted(false);
commandSetGeDetectorWidth = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setGeDetectorWidth",this);
commandSetGeDetectorWidth->SetGuidance("Set width of target samples");
commandSetGeDetectorWidth->SetParameterName("SampleWidth",false);
commandSetGeDetectorWidth->SetRange("SampleWidth>0.");
commandSetGeDetectorWidth->SetUnitCategory("Length");
commandSetGeDetectorWidth->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetGeDetectorWidth->SetToBeBroadcasted(false);
commandSetContactThickness = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setContactThickness",this);
commandSetContactThickness->SetGuidance("Set thickness of Al contact");
commandSetContactThickness->SetParameterName("ContactThickness",false);
commandSetContactThickness->SetRange("ContactThickness>0.");
commandSetContactThickness->SetUnitCategory("Length");
commandSetContactThickness->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetContactThickness->SetToBeBroadcasted(false);
commandSetBeWindowRadius = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setBeWindowRadius",this);
commandSetBeWindowRadius->SetGuidance("Set radius of Be Window");
commandSetBeWindowRadius->SetParameterName("WindowRadius",false);
commandSetBeWindowRadius->SetRange("WindowRadius>0.");
commandSetBeWindowRadius->SetUnitCategory("Length");
commandSetBeWindowRadius->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetBeWindowRadius->SetToBeBroadcasted(false);
commandSetCollimatorThickness = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setCollimatorThickness",this);
commandSetCollimatorThickness->SetGuidance("Set thickness of target samples");
commandSetCollimatorThickness->SetParameterName("CollimatorThickness",false);
commandSetCollimatorThickness->SetRange("CollimatorThickness>0.");
commandSetCollimatorThickness->SetUnitCategory("Length");
commandSetCollimatorThickness->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetCollimatorThickness->SetToBeBroadcasted(false);
commandSetDistanceCollimatorDetector = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/det/setDistanceCollimatorDetector",this);
commandSetDistanceCollimatorDetector->SetGuidance("Set distance collimator detector");
commandSetDistanceCollimatorDetector->SetParameterName("CollimatorDistance",false);
commandSetDistanceCollimatorDetector->SetRange("CollimatorDistance>0.");
commandSetDistanceCollimatorDetector->SetUnitCategory("Length");
commandSetDistanceCollimatorDetector->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetDistanceCollimatorDetector->SetToBeBroadcasted(false);
commandSetDetectorType = new G4UIcmdWithAnInteger("/Ge2SOS/det/setDetectorType",this);
commandSetDetectorType->SetGuidance("Set detector type");
commandSetDetectorType->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetDetectorType->SetToBeBroadcasted(false);
commandSetNumberOfTargetSamples = new G4UIcmdWithAnInteger("/Ge2SOS/det/setNTargetSamples",this);
commandSetNumberOfTargetSamples->SetGuidance("Set number of target Samples");
commandSetNumberOfTargetSamples->AvailableForStates(G4State_PreInit,G4State_Idle);
commandSetNumberOfTargetSamples->SetToBeBroadcasted(false);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
DetectorMessenger::~DetectorMessenger()
{
delete fDetDir;
delete fGe2SOSDir;
delete commandSetWorldMaterial;
delete commandSetDetectorType;
delete commandSetNumberOfTargetSamples;
delete commandSetGeDetectorLength;
delete commandSetGeDetectorThickness;
delete commandSetGeDetectorWidth;
delete commandSetContactThickness;
delete commandSetBeWindowRadius;
delete commandSetCollimatorThickness;
delete commandSetDistanceCollimatorDetector;
delete commandSetCollimatorMaterial;
delete commandSetTargetMaterial;
delete commandSetGeContainerMaterial;
delete commandSetGeContainerMaterialCoating;
delete commandSetDetectorName;
delete commandSetSetupName;
delete commandSetDataType;
delete commandSetOutputDirectory;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
{
if( command == commandSetTargetMaterial ){
fDetector->SetTargetMaterial(newValue);
}
if( command == commandSetGeContainerMaterial ){
fDetector->SetGeContainerMaterial(newValue);
}
if( command == commandSetGeContainerMaterialCoating ){
fDetector->SetGeContainerMaterialCoating(newValue);
}
if( command == commandSetDetectorName ){
fDetector->SetDetectorName(newValue);
}
if( command == commandSetSetupName ){
fDetector->SetSetupName(newValue);
}
if( command == commandSetDataType ){
fDetector->SetDataType(newValue);
}
if( command == commandSetCollimatorMaterial ){
fDetector->SetCollimatorMaterial(newValue);
}
if( command == commandSetOutputDirectory ){
fDetector->SetOutputDirectory(newValue);
}
if( command == commandSetWorldMaterial ){
fDetector->SetWorldMaterial(newValue);
}
if( command == commandSetGeDetectorLength ){
fDetector->SetGeDetectorLength(commandSetGeDetectorLength->GetNewDoubleValue(newValue));
}
if( command == commandSetGeDetectorThickness ){
fDetector->SetGeDetectorThickness(commandSetGeDetectorThickness->GetNewDoubleValue(newValue));
}
if( command == commandSetContactThickness ){
fDetector->SetContactThickness(commandSetContactThickness->GetNewDoubleValue(newValue));
}
if( command == commandSetGeDetectorWidth ){
fDetector->SetGeDetectorWidth(commandSetGeDetectorWidth->GetNewDoubleValue(newValue));
}
if( command == commandSetBeWindowRadius ){
fDetector->SetBeWindowRadius(commandSetBeWindowRadius->GetNewDoubleValue(newValue));
}
if( command == commandSetCollimatorThickness ){
fDetector->SetCollimatorThickness(commandSetCollimatorThickness->GetNewDoubleValue(newValue));
}
if( command == commandSetDistanceCollimatorDetector ){
fDetector->SetDistanceCollimatorDetector(commandSetDistanceCollimatorDetector->GetNewDoubleValue(newValue));
}
if( command == commandSetDetectorType ){
fDetector->SetDetectorType(commandSetDetectorType->GetNewIntValue(newValue));
}
if( command == commandSetNumberOfTargetSamples ){
fDetector->SetNumberOfTargetSamples(commandSetNumberOfTargetSamples->GetNewIntValue(newValue));
}
}
#include "EventAction.hh"
#include "RunAction.hh"
#include "G4SystemOfUnits.hh"
#include "G4PhysicalConstants.hh"
#include "DetectorConstruction.hh"
#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4LogicalVolume.hh"
#include "G4Box.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventAction::EventAction(DetectorConstruction* det, RunAction* runAction)
: G4UserEventAction(),fDetector(det),
fRunAction(runAction)
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventAction::~EventAction()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventAction::BeginOfEventAction(const G4Event* myEvent)
{
depositedEnergyPENStackedSample1 = 0.;
eventNumber = myEvent->GetEventID();
xFirstPen = 0.;
yFirstPen = 0.;
zFirstPen = 0.;
if (myEvent->GetEventID() % 10000 == 0)
G4cout << "starting event no.: " << myEvent->GetEventID() << G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventAction::AddInfo(G4double xF, G4double yF, G4double zF, G4double depE){
if (depE > 1.0){
auto analysisManager = G4AnalysisManager::Instance();
analysisManager->FillNtupleIColumn(0, eventNumber);
analysisManager->FillNtupleDColumn(1, xF);
analysisManager->FillNtupleDColumn(2, yF);
analysisManager->FillNtupleDColumn(3, zF);
analysisManager->FillNtupleDColumn(4, depE);
analysisManager->AddNtupleRow(0);
}
}
void EventAction::EndOfEventAction(const G4Event* myEvent)
{
if (myEvent->GetEventID() % 10000 == 0)
G4cout << "finish event no.: " << myEvent->GetEventID() << G4endl;
}
#include "globals.hh"
#include "PhysicsList.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4BosonConstructor.hh"
#include "G4LeptonConstructor.hh"
#include "G4MesonConstructor.hh"
#include "G4BaryonConstructor.hh"
#include "G4IonConstructor.hh"
#include "G4ShortLivedConstructor.hh"
#include "G4EmParameters.hh"
#include "G4ProcessManager.hh"
#include "G4Cerenkov.hh"
#include "G4Scintillation.hh"
#include "G4OpAbsorption.hh"
#include "G4OpRayleigh.hh"
#include "G4OpMieHG.hh"
#include "G4OpBoundaryProcess.hh"
#include "G4LossTableManager.hh"
#include "G4EmSaturation.hh"
#include "G4UnitsTable.hh"
#include "G4ParticleTypes.hh"
#include "G4PhysicsListHelper.hh"
#include "G4RadioactiveDecay.hh"
#include "G4SystemOfUnits.hh"
#include "G4NuclideTable.hh"
#include "G4LossTableManager.hh"
#include "G4UAtomicDeexcitation.hh"
#include "G4NuclearLevelData.hh"
#include "G4DeexPrecoParameters.hh"
G4ThreadLocal G4int PhysicsList::fVerboseLevel = 0;
G4ThreadLocal G4int PhysicsList::fMaxNumPhotonStep = 20;
G4ThreadLocal G4Cerenkov* PhysicsList::fCerenkovProcess = 0;
G4ThreadLocal G4Scintillation* PhysicsList::fScintillationProcess = 0;
G4ThreadLocal G4OpAbsorption* PhysicsList::fAbsorptionProcess = 0;
G4ThreadLocal G4OpRayleigh* PhysicsList::fRayleighScatteringProcess = 0;
G4ThreadLocal G4OpMieHG* PhysicsList::fMieHGScatteringProcess = 0;
G4ThreadLocal G4OpBoundaryProcess* PhysicsList::fBoundaryProcess = 0;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PhysicsList::PhysicsList()
: G4VUserPhysicsList() {
//modified by Luis
defaultCutValue = 1.0*micrometer; //
cutForGamma = defaultCutValue;
cutForElectron = 1.0*micrometer;
cutForPositron = defaultCutValue;
VerboseLevel = 0;
OpVerbLevel = 0;
SetVerboseLevel(VerboseLevel);
//end modifs
//read new PhotonEvaporation data set
// mandatory for G4NuclideTable
//
G4NuclideTable::GetInstance()->SetThresholdOfHalfLife(0.1*picosecond);
G4NuclideTable::GetInstance()->SetLevelTolerance(1.0*eV);
//
G4DeexPrecoParameters* deex =
G4NuclearLevelData::GetInstance()->GetParameters();
deex->SetCorrelatedGamma(false);
deex->SetStoreAllLevels(true);
deex->SetIsomerProduction(true);
deex->SetMaxLifeTime(G4NuclideTable::GetInstance()->GetThresholdOfHalfLife()
/std::log(2.));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PhysicsList::~PhysicsList()
{;}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::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 bConstructor;
bConstructor.ConstructParticle();
G4LeptonConstructor lConstructor;
lConstructor.ConstructParticle();
G4MesonConstructor mConstructor;
mConstructor.ConstructParticle();
G4BaryonConstructor rConstructor;
rConstructor.ConstructParticle();
G4IonConstructor iConstructor;
iConstructor.ConstructParticle();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructProcess()
{
AddTransportation();
ConstructDecay();
ConstructEM();
ConstructOp();
G4EmParameters* param = G4EmParameters::Instance();
param->SetFluo(true);
param->SetAuger(true);
param->SetPixe(true);
SetVerbose(0);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4Decay.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructDecay()
{
// Add Decay Process
G4Decay* theDecayProcess = new G4Decay();
auto particleIterator=GetParticleIterator();
particleIterator->reset();
while( (*particleIterator)() ){
G4ParticleDefinition* particle = particleIterator->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);
}
}
G4RadioactiveDecay* radioactiveDecay = new G4RadioactiveDecay();
//radioactiveDecay->SetARM(false); //Atomic Rearangement
radioactiveDecay->SetARM(true); //Atomic Rearangement
G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
ph->RegisterProcess(radioactiveDecay, G4GenericIon::GenericIon());
// Need to initialize atomic deexcitation outside of radioactive decay
G4LossTableManager* theManager = G4LossTableManager::Instance();
G4VAtomDeexcitation* p = theManager->AtomDeexcitation();
if (!p) {
G4UAtomicDeexcitation* atomDeex = new G4UAtomicDeexcitation();
theManager->SetAtomDeexcitation(atomDeex);
atomDeex->InitialiseAtomicDeexcitation();
}
//
// mandatory for G4NuclideTable
//
G4NuclideTable::GetInstance()->SetThresholdOfHalfLife(0.1*picosecond);
G4NuclideTable::GetInstance()->SetLevelTolerance(1.0*eV);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4EmLivermorePolarizedPhysics.hh"
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"
#include "G4eMultipleScattering.hh"
#include "G4MuMultipleScattering.hh"
#include "G4hMultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4LivermoreBremsstrahlungModel.hh"
#include "G4LivermoreIonisationModel.hh"
#include "G4eplusAnnihilation.hh"
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4hIonisation.hh"
//////
//
#include "G4PolarizedAnnihilation.hh"
#include "G4PolarizedBremsstrahlung.hh"
#include "G4PolarizedCompton.hh"
#include "G4PolarizedGammaConversion.hh"
#include "G4PolarizedIonisation.hh"
#include "G4PolarizedPhotoElectric.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructEM()
{
emPhysicsList = new G4EmLivermorePolarizedPhysics;
emPhysicsList->ConstructProcess();
/*
auto particleIterator=GetParticleIterator();
particleIterator->reset();
while( (*particleIterator)() ){
G4ParticleDefinition* particle = particleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (particleName == "gamma") {
// gamma
pmanager->AddDiscreteProcess(new G4PolarizedPhotoElectric);
pmanager->AddDiscreteProcess(new G4PolarizedCompton);
pmanager->AddDiscreteProcess(new G4PolarizedGammaConversion);
// Construct processes for gamma
//pmanager->AddDiscreteProcess(new G4GammaConversion());
//pmanager->AddDiscreteProcess(new G4ComptonScattering());
//pmanager->AddDiscreteProcess(new G4PhotoElectricEffect());
} else if (particleName == "e-") {
//electron
//Modified by Luis
G4eMultipleScattering* msc = new G4eMultipleScattering();
msc->SetStepLimitType(fUseDistanceToBoundary);
pmanager->AddProcess(msc,-1, 1, -1);
// Ionisation
G4eIonisation* eIonisation = new G4eIonisation();
eIonisation->SetEmModel(new G4LivermoreIonisationModel());
eIonisation->SetStepFunction(0.2, 100*um); //improved precision in tracking
pmanager->AddProcess(eIonisation,-1, 2, 2);
// Bremsstrahlung
G4eBremsstrahlung* eBremsstrahlung = new G4eBremsstrahlung();
eBremsstrahlung->SetEmModel(new G4LivermoreBremsstrahlungModel());
pmanager->AddProcess(eBremsstrahlung, -1,-3, 3);
//End of modifs by Luis
//Construct processes for 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
// Construct processes for 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
// Construct processes for 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 ((particle->GetPDGCharge() != 0.0) &&
(particle->GetParticleName() != "chargedgeantino") &&
!particle->IsShortLived()) {
// 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 "G4Threading.hh"
void PhysicsList::ConstructOp()
{
fCerenkovProcess = new G4Cerenkov("Cerenkov");
fCerenkovProcess->SetMaxNumPhotonsPerStep(fMaxNumPhotonStep);
fCerenkovProcess->SetMaxBetaChangePerStep(10.0);
fCerenkovProcess->SetTrackSecondariesFirst(true);
fScintillationProcess = new G4Scintillation("Scintillation");
//fScintillationProcess->SetScintillationYieldFactor(1.0);
//fScintillationProcess->SetScintillationExcitationRatio(0.0);//Line added by Luis
fScintillationProcess->SetTrackSecondariesFirst(true);
fScintillationProcess->SetVerboseLevel(OpVerbLevel);
// scintillation process for alpha:
G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
// theScintProcessNuc->DumpPhysicsTable();
theScintProcessAlpha->SetTrackSecondariesFirst(true);
//theScintProcessAlpha->SetScintillationYieldFactor(1.1);
//theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
theScintProcessAlpha->SetVerboseLevel(OpVerbLevel);
// scintillation process for heavy nuclei
G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
// theScintProcessNuc->DumpPhysicsTable();
theScintProcessNuc->SetTrackSecondariesFirst(true);
//theScintProcessNuc->SetScintillationYieldFactor(0.2);
//theScintProcessNuc->SetScintillationExcitationRatio(1.0);
theScintProcessNuc->SetVerboseLevel(OpVerbLevel);
fAbsorptionProcess = new G4OpAbsorption();
fRayleighScatteringProcess = new G4OpRayleigh();
fMieHGScatteringProcess = new G4OpMieHG();
fBoundaryProcess = new G4OpBoundaryProcess();
fCerenkovProcess->SetVerboseLevel(fVerboseLevel);
fScintillationProcess->SetVerboseLevel(OpVerbLevel);
fAbsorptionProcess->SetVerboseLevel(OpVerbLevel);
fRayleighScatteringProcess->SetVerboseLevel(fVerboseLevel);
fMieHGScatteringProcess->SetVerboseLevel(fVerboseLevel);
fBoundaryProcess->SetVerboseLevel(fVerboseLevel);
// Use Birks Correction in the Scintillation process
if(G4Threading::IsMasterThread())
{
G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
fScintillationProcess->AddSaturation(emSaturation);
}
auto particleIterator=GetParticleIterator();
particleIterator->reset();
while( (*particleIterator)() ){
G4ParticleDefinition* particle = particleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (fCerenkovProcess->IsApplicable(*particle)) {
pmanager->AddProcess(fCerenkovProcess);
pmanager->SetProcessOrdering(fCerenkovProcess,idxPostStep);
}
//if (fScintillationProcess->IsApplicable(*particle) && particle->GetParticleType() != "nucleus") {
if (fScintillationProcess->IsApplicable(*particle)) {
if(particle->GetParticleName() == "GenericIon") {
pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete
pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest);
pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep);
}
else if(particle->GetParticleName() == "alpha") {
pmanager->AddProcess(theScintProcessAlpha);
pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest);
pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep);
}
else {
pmanager->AddProcess(fScintillationProcess);
pmanager->SetProcessOrderingToLast(fScintillationProcess, idxAtRest);
pmanager->SetProcessOrderingToLast(fScintillationProcess, idxPostStep);
}
}
if (particleName == "opticalphoton") {
G4cout << " AddDiscreteProcess to OpticalPhoton " << G4endl;
pmanager->AddDiscreteProcess(fAbsorptionProcess);
pmanager->AddDiscreteProcess(fRayleighScatteringProcess);
pmanager->AddDiscreteProcess(fMieHGScatteringProcess);
pmanager->AddDiscreteProcess(fBoundaryProcess);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::SetVerbose(G4int verbose)
{
fVerboseLevel = verbose;
fCerenkovProcess->SetVerboseLevel(fVerboseLevel);
fScintillationProcess->SetVerboseLevel(OpVerbLevel);
fAbsorptionProcess->SetVerboseLevel(fVerboseLevel);
fRayleighScatteringProcess->SetVerboseLevel(fVerboseLevel);
fMieHGScatteringProcess->SetVerboseLevel(fVerboseLevel);
fBoundaryProcess->SetVerboseLevel(fVerboseLevel);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::SetNbOfPhotonsCerenkov(G4int MaxNumber)
{
fMaxNumPhotonStep = MaxNumber;
fCerenkovProcess->SetMaxNumPhotonsPerStep(fMaxNumPhotonStep);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::SetCuts()
{
// " G4VUserPhysicsList::SetCutsWithDefault" method sets
// the default cut value for all particle types
//
//SetCutsWithDefault();
//special for low energy physics
G4double lowlimit=250*eV;
G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit,100.*GeV);
// 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(cutForGamma, "gamma");
SetCutValue(cutForElectron, "e-");
SetCutValue(cutForPositron, "e+");
// if (verboseLevel>0) DumpCutValuesTable();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "PrimaryGeneratorAction.hh"
#include "PrimaryGeneratorMessenger.hh"
#include "G4MTRunManager.hh"
#include "Randomize.hh"
#include "DetectorConstruction.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4SystemOfUnits.hh"
#include "G4IonTable.hh"
#include "G4DecayTable.hh"
#include "G4OpticalPhoton.hh"
#include "G4VPhysicalVolume.hh"
#include "G4PhysicalConstants.hh"
#include "G4Navigator.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* det)
: G4VUserPrimaryGeneratorAction()
{
fDetector = det;
fPrimaryMessenger = new PrimaryGeneratorMessenger(this);
G4int n_particle = 1;
fParticleGun = new G4ParticleGun(n_particle);
fPositionX = 0.0 *mm;
fPositionY = 0.0 *mm;
fPositionZ = 30.*mm;
fDiameter = 0.1*mm;
fSourceType = 0;
fSourceDirectionType = 0;
fSourceEnergy = 10*keV;
fTheta_polar = 0; //in degrees, so the default is horizontal polarization
fPolarization_degree = 0.95; // values from 0 to 1, so the default is 95% polarization
fPhotonWavelength = 0;
fParticleName = "void";
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
delete fParticleGun;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Randomises placement and momentum vectors for required sources.
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4double ionCharge = 0.*eplus;
G4double excitEnergy = 0.*eV;
G4int Z=0, A=0;
G4ParticleDefinition* ion;
//Define a more realistic source geometry
G4double r = fDiameter * std::sqrt(G4UniformRand()) / 2.;
G4double theta = 2. * pi * G4UniformRand()*rad;
G4double x_sourceframe = fPositionX + r*cos(theta)*CLHEP::mm;
G4double y_sourceframe = fPositionY + r*sin(theta)*CLHEP::mm;
G4double z_sourceframe = r*sin(theta)*CLHEP::mm;
//G4cout<<" r "<<r<<" theta: "<<theta<<" x "<<x_sourceframe<<" y "<<y_sourceframe<<G4endl;
G4ThreeVector position = G4ThreeVector(x_sourceframe, y_sourceframe, fPositionZ);
//For random direction in order to optimize we generate random point in a sphere, but we are interested only in the values close to phi = 180 degrees, which is why we use random()/20 -1
theta = 2. * pi * G4UniformRand()*rad;
G4double phi = acos(G4UniformRand()/20. -1.);
G4double x_direction = sin(phi)*cos(theta);
G4double y_direction = sin(phi)*sin(theta);
G4double z_direction = cos(phi);
G4ThreeVector direction = G4ThreeVector(0,0,-1);
//Polarazation of gammas
G4double theta_polar_radians = fTheta_polar * pi * rad / 180;
G4double stoke_vector_e1 = cos(theta_polar_radians);
G4double stoke_vector_e2 = sin(theta_polar_radians);
G4ThreeVector stokes_vector = G4ThreeVector(stoke_vector_e1,stoke_vector_e2,0);
if( fPolarization_degree < G4UniformRand()){
stokes_vector = G4ThreeVector(0,0,0);
}
if(fSourceDirectionType == 1){
direction = G4ThreeVector(x_direction,y_direction,z_direction);
}
if(fDetector->GetDetectorType() == 2 || fDetector->GetDetectorType() == 3){
position = G4ThreeVector(fPositionZ, y_sourceframe/4., fDetector->GetSamplePositionZ() + z_sourceframe/4.);
direction = G4ThreeVector(-1,0,0);
G4cout<<" sample position "<<fDetector->GetSamplePositionZ()<<G4endl;
}
switch (fSourceType) {
case 0:
fParticleGun->SetParticleDefinition(particleTable->FindParticle("gamma"));
fParticleGun->SetParticleEnergy(fSourceEnergy);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(direction);
fParticleGun->SetParticlePolarization(stokes_vector);
break;
case 1:
Z = 26;
A = 55;
ion = G4IonTable::GetIonTable()->GetIon(Z,A,excitEnergy);
fParticleGun->SetParticleEnergy(0.*eV);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,0.));
fParticleGun->SetParticleDefinition(ion);
fParticleGun->SetParticleCharge(ionCharge);
break;
case 2:
Z = 55;
A = 137;
ion = G4IonTable::GetIonTable()->GetIon(Z,A,excitEnergy);
fParticleGun->SetParticleEnergy(0.*eV);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,0.));
fParticleGun->SetParticleDefinition(ion);
fParticleGun->SetParticleCharge(ionCharge);
break;
case 3:
Z = 83;
A = 207;
ion = G4IonTable::GetIonTable()->GetIon(Z,A,excitEnergy);
fParticleGun->SetParticleEnergy(0.*eV);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,0.));
fParticleGun->SetParticleDefinition(ion);
fParticleGun->SetParticleCharge(ionCharge);
break;
case 4:
Z = 38;
A = 90;
ion = G4IonTable::GetIonTable()->GetIon(Z,A,excitEnergy);
fParticleGun->SetParticleEnergy(0.*eV);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,0.));
fParticleGun->SetParticleDefinition(ion);
fParticleGun->SetParticleCharge(ionCharge);
break;
case 5:
//Z = 95;
//A = 241;
//ion = G4IonTable::GetIonTable()->GetIon(Z,A,excitEnergy);
fParticleGun->SetParticleDefinition(particleTable->FindParticle("gamma"));
fParticleGun->SetParticleEnergy(59.5*keV);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(x_direction,y_direction,z_direction));
//fParticleGun->SetParticleDefinition(ion);
//fParticleGun->SetParticleCharge(ionCharge);
break;
case 6:
fParticleGun->SetParticleDefinition(particleTable->FindParticle("e-"));
fParticleGun->SetParticleEnergy(fSourceEnergy);
fParticleGun->SetParticlePosition(position);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0,0,-1));
break;
case 7:
fParticleGun->SetParticleDefinition(particleTable->FindParticle("opticalphoton"));
fParticleGun->SetParticleEnergy(fSourceEnergy);
fParticleGun->SetParticlePosition(G4ThreeVector(fPositionX,fPositionY,fPositionZ));
fParticleGun->SetParticleMomentumDirection(direction);
break;
}
fParticleGun->GeneratePrimaryVertex(anEvent);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::SetParticleName(G4int Z, G4int A, G4double excitEnergy)
{
fParticleName = G4IonTable::GetIonTable()->GetIonName(Z,A,excitEnergy);
}
void PrimaryGeneratorAction::SetSourcePositionX(G4double newPosition){
fPositionX = newPosition;
}
void PrimaryGeneratorAction::SetSourcePositionY(G4double newPosition){
fPositionY = newPosition;
}
void PrimaryGeneratorAction::SetSourcePositionZ(G4double newPosition){
fPositionZ = newPosition;
}
void PrimaryGeneratorAction::SetSourceDiameter(G4double newDiameter){
fDiameter = newDiameter;
}
void PrimaryGeneratorAction::SetSourcePolarizationAngle(G4double newAngle){
fTheta_polar = newAngle;
}
void PrimaryGeneratorAction::SetSourcePolarizationDegree(G4double newPolarDegree){
fPolarization_degree = newPolarDegree;
}
void PrimaryGeneratorAction::SetSourceType(G4int newType)
{
if (newType <= 5 && newType >= 0)
{
fSourceType = newType;
}
else
{
G4cerr << "The option is out of the possible values (0-6)!" << G4endl;
G4cerr << "The default option (0) is set!" << G4endl;
fSourceType = 0;
}
}
void PrimaryGeneratorAction::SetSourceDirectionType(G4int newType)
{
if (newType <= 1 && newType >= 0)
{
fSourceDirectionType = newType;
}
else
{
G4cerr << "The option is out of the possible values (0-1)!" << G4endl;
G4cerr << "The default option (0) is set!" << G4endl;
fSourceType = 0;
}
}
void PrimaryGeneratorAction::SetSourceGeometry(G4int newType)
{
if (newType <= 2 && newType >= 0){
fSourceGeometry = newType;}
else{
G4cerr<<"Possible values are 0 for square and 1 for circunference"<<G4endl;
G4cerr<<"Setting the default value to square "<<G4endl;
fSourceGeometry = 0;
}
}
void PrimaryGeneratorAction::SetSourceEnergy(G4double newEnergy)
{
if (newEnergy>0)
{
fSourceEnergy = newEnergy;
}
else{
G4cerr << "New energy is < 0." << G4endl;
G4cerr << "The default option 10 keV is set!" << G4endl;
fSourceEnergy = 10.*keV;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "PrimaryGeneratorMessenger.hh"
#include "PrimaryGeneratorAction.hh"
#include "G4UIdirectory.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWithADouble.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4SystemOfUnits.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorMessenger::PrimaryGeneratorMessenger(PrimaryGeneratorAction* Gun)
: G4UImessenger(),
fAction(Gun),fGunDir(0)
{
fGunDir = new G4UIdirectory("/Ge2SOS/gun/");
fGunDir->SetGuidance("PrimaryGenerator control");
fSourceType = new G4UIcmdWithAnInteger("/Ge2SOS/gun/sourceType",this);
fSourceType->SetGuidance("Choose the type of source");
fSourceType->SetParameterName("sourceType",true);
fSourceType->SetDefaultValue(0);
fSourceType->AvailableForStates(G4State_Idle);
fSourceDirectionType = new G4UIcmdWithAnInteger("/Ge2SOS/gun/sourceDirectionType",this);
fSourceDirectionType->SetGuidance("Choose the direction type");
fSourceDirectionType->SetParameterName("sourceDirectionType",true);
fSourceDirectionType->SetDefaultValue(0);
fSourceDirectionType->AvailableForStates(G4State_Idle);
fSourceGeometry = new G4UIcmdWithAnInteger("/Ge2SOS/gun/sourceGeometry",this);
fSourceGeometry->SetGuidance("Choose the type of geometry for the source");
fSourceGeometry->SetParameterName("sourceGeometry",true);
fSourceGeometry->SetDefaultValue(0);
fSourceGeometry->AvailableForStates(G4State_Idle);
fSourceEnergy = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/gun/sourceEnergy", this);
fSourceEnergy->SetGuidance("Choose source energy");
fSourceEnergy->SetParameterName("sourceEnergy",true);
fSourceEnergy->SetDefaultValue(20.*keV);
fSourceEnergy->AvailableForStates(G4State_Idle);
fSourcePolarizationAngle = new G4UIcmdWithADouble("/Ge2SOS/gun/sourceGammaPolarizationAngle", this);
fSourcePolarizationAngle->SetGuidance("Choose angle of polarization in degrees");
fSourcePolarizationAngle->SetParameterName("sourcePolarAngle",true);
fSourcePolarizationAngle->SetDefaultValue(0.);
fSourcePolarizationAngle->AvailableForStates(G4State_Idle);
fSourcePolarizationDegree = new G4UIcmdWithADouble("/Ge2SOS/gun/sourceGammaPolarizationDegree", this);
fSourcePolarizationDegree->SetGuidance("Choose degree of polarization o to 1");
fSourcePolarizationDegree->SetParameterName("sourcePolarDegree",true);
fSourcePolarizationDegree->SetDefaultValue(1.);
fSourcePolarizationDegree->AvailableForStates(G4State_Idle);
fSourcePositionX = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/gun/sourcePositionX",this);
fSourcePositionX->SetGuidance("Set Source x position");
fSourcePositionX->SetParameterName("fPositionX",true);
fSourcePositionX->SetDefaultValue(0.*mm);
fSourcePositionX->AvailableForStates(G4State_Idle);
fSourcePositionY = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/gun/sourcePositionY",this);
fSourcePositionY->SetGuidance("Set Source y position");
fSourcePositionY->SetParameterName("fPositionY",true);
fSourcePositionY->SetDefaultValue(0.*mm);
fSourcePositionY->AvailableForStates(G4State_Idle);
fSourcePositionZ = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/gun/sourcePositionZ",this);
fSourcePositionZ->SetGuidance("Set Source z position");
fSourcePositionZ->SetParameterName("fPositionZ",true);
fSourcePositionZ->SetDefaultValue(30.*mm);
fSourcePositionZ->AvailableForStates(G4State_Idle);
fSourceDiameter = new G4UIcmdWithADoubleAndUnit("/Ge2SOS/gun/sourceDiameter",this);
fSourceDiameter->SetGuidance("Set Source diameter or side");
fSourceDiameter->SetParameterName("fDiameter",true);
fSourceDiameter->SetDefaultValue(5.*mm);
fSourceDiameter->AvailableForStates(G4State_Idle);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorMessenger::~PrimaryGeneratorMessenger()
{
delete fSourceType;
delete fSourceDirectionType;
delete fSourceGeometry;
delete fSourceEnergy;
delete fSourcePolarizationAngle;
delete fSourcePolarizationDegree;
delete fGunDir;
delete fSourcePositionX;
delete fSourcePositionY;
delete fSourcePositionZ;
delete fSourceDiameter;
//delete fAction;//debug
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorMessenger::SetNewValue(G4UIcommand* command, G4String newValue)
{
if(command == fSourceType) {
fAction->SetSourceType(fSourceType->GetNewIntValue(newValue));
}
if(command == fSourceDirectionType) {
fAction->SetSourceDirectionType(fSourceDirectionType->GetNewIntValue(newValue));
}
if(command == fSourceGeometry) {
fAction->SetSourceGeometry(fSourceGeometry->GetNewIntValue(newValue));
}
if(command == fSourceEnergy){
fAction->SetSourceEnergy(fSourceEnergy->GetNewDoubleValue(newValue));
}
if(command == fSourcePolarizationAngle){
fAction->SetSourcePolarizationAngle(fSourcePolarizationAngle->GetNewDoubleValue(newValue));
}
if(command == fSourcePolarizationDegree){
fAction->SetSourcePolarizationDegree(fSourcePolarizationDegree->GetNewDoubleValue(newValue));
}
if(command == fSourcePositionX){
fAction->SetSourcePositionX(fSourcePositionX->GetNewDoubleValue(newValue));
}
if(command == fSourcePositionY){
fAction->SetSourcePositionY(fSourcePositionY->GetNewDoubleValue(newValue));
}
if(command == fSourcePositionZ){
fAction->SetSourcePositionZ(fSourcePositionZ->GetNewDoubleValue(newValue));
}
if(command == fSourceDiameter){
fAction->SetSourceDiameter(fSourceDiameter->GetNewDoubleValue(newValue));
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Make this appear first!
#include "G4Timer.hh"
#include "RunAction.hh"
#include "RunActionMessenger.hh"
#include "G4Run.hh"
#include "G4Types.hh"
#include "G4Material.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "EventAction.hh"
#include <string>
#include <ctime>
#include <sys/stat.h>
#include "DetectorConstruction.hh"
#include "PrimaryGeneratorAction.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
: G4UserRunAction(),fDetector(det),fPrimary(kin),
fTimer(0)
{
fTimer = new G4Timer;
fMan = G4AnalysisManager::Instance();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
RunAction::~RunAction()
{
delete fTimer;
}
std::string datetime()
{
time_t rawtime;
struct tm * timeinfo;
char buffer[80];
time (&rawtime);
timeinfo = localtime(&rawtime);
strftime(buffer,80,"%d_%m_%Y_%H_%M",timeinfo);
return std::string(buffer);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void RunAction::BeginOfRunAction(const G4Run* aRun)
{
fMan = G4AnalysisManager::Instance();
G4int sourceName = fPrimary->GetSourceType();
G4String sourceString;
switch (sourceName){
case 0:
sourceString = "gamma";
break;
case 1:
sourceString = "Fe55";
break;
case 2:
sourceString = "Cs-137";
break;
case 3:
sourceString = "Bi-207";
break;
case 4:
sourceString = "Sr-90";
break;
case 5:
sourceString = "Am-241";
break;
case 6:
sourceString = "optical photon";
break;
}
G4cout << sourceString << G4endl;
//G4String positionZ = G4BestUnit((fDetector->GetGeDetectorThickness()-fPrimary->GetSourcePositionZ()),"Length");
G4String positionX ="_x_source_"+std::to_string((int)fPrimary->GetSourcePositionX())+"mm";
G4String sourceEnergy =std::to_string((int)(1000. * fPrimary->GetSourceEnergy()))+"_keV";
G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
fTimer->Start();
//Target material
G4String s_setup_Name =fDetector->GetSetupName();
//n Samples
//Thickness
//std::stringstream stream_thickness;
//stream_thickness << std::fixed << std::setprecision(2) <<fDetector->GetGeDetectorThickness();
G4String s_target_thickness ="_GeThickness_"+std::to_string((int)round(fDetector->GetGeDetectorThickness()))+"mm";
G4String s_collimator_material =fDetector->GetCollimatorMaterial()->GetName();
G4String s_collimatorThickness = "_collimator_"+s_collimator_material+"_"+std::to_string((int)fDetector->GetCollimatorThickness())+"mm";
G4String s_setup_type = "_setup_"+std::to_string((int)fDetector->GetDetectorType());
if(fDetector->GetDetectorType()==0 || fDetector->GetDetectorType()==2){
s_collimatorThickness = "_no_collimator";
}
G4String directorName = fDetector->GetDetectorOutputDataDirectory()+"/"+sourceString+"_"+sourceEnergy
//+datetime()
+positionX
+s_target_thickness
+s_setup_type
+s_collimatorThickness
+"/";
mkdir(directorName, 0777);
fFileName = directorName+fDetector->GetDetectorName()
+"_"+s_setup_Name+"."
+fDetector->GetDataType();
fMan->SetVerboseLevel(0);
fMan->OpenFile(fFileName);
fMan->CreateNtuple("GeSOS","Detailed MC Info");
fMan->CreateNtupleIColumn("EventID");
fMan->CreateNtupleDColumn("x");
fMan->CreateNtupleDColumn("y");
fMan->CreateNtupleDColumn("z");
fMan->CreateNtupleDColumn("E");
fMan->FinishNtuple();
}
void RunAction::SetFileName(G4String fileName)
{
fFileName = fileName;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void RunAction::EndOfRunAction(const G4Run* aRun)
{
fTimer->Stop();
G4cout << "number of event = " << aRun->GetNumberOfEvent()
<< " " << *fTimer << G4endl;
G4cout << "End Run" << G4endl;
fMan->Write();
fMan->CloseFile();
//delete fMan;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4UnionSolid.hh"
#include "G4SubtractionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4UnionSolid.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "SiliconPlateConstruction.hh"
G4VSolid* SiliconPlateConstruction::ConstructPlate(){
fSiliconPlate_h = 1.5*mm;
fHolderWidth=90.00*mm;
G4double TubsStartAngle = 0;
G4double TubsSpanningAngle = 360 * deg;
G4double trg_b = (fHolderWidth/sqrt(2.))*mm;
G4double trg_h = 45.*mm;
G4double rect_x = 2.*(trg_h + 9.5*mm);
G4double rect_y = fHolderWidth;
G4double rect_bc_x = 28.*mm;
G4double cut_x = 20.*mm;
G4double cut_y = 31.*mm;
G4RotationMatrix* rm = new G4RotationMatrix();
rm->rotateZ(45.*deg);
G4RotationMatrix* rm1 = new G4RotationMatrix();
rm1->rotateZ(-120.*deg);
G4RotationMatrix* rm2 = new G4RotationMatrix();
rm2->rotateZ(120.*deg);
char solidname[100];
G4VSolid* solid_LowerPlate_Rect = new G4Box(solidname,
rect_x/2,
rect_y/2,
0.5*fSiliconPlate_h);
// build rectangular cuts
G4VSolid* solid_LowerPlate_CutRect0 = new G4Box(solidname,
(0.5*rect_x- rect_bc_x+1*mm)*mm/2,
(rect_y+1*mm)/2,
0.5*(fSiliconPlate_h+.1*mm));
G4VSolid* solid_LowerPlate_CutRect1 = new G4Box(solidname,
trg_b/2,
trg_b/2,
0.5*(fSiliconPlate_h+.1*mm));
G4VSolid* solid_LowerPlate_CutRect2 = new G4Box(solidname,
cut_x/2,
cut_y/2,
0.5*(fSiliconPlate_h+.1*mm));
//build rectangular cut plates (to imitate rounded corners, important to fit into MS)
G4VSolid* solid_LowerPlate_CutRect3 = new G4Box(solidname,
9*mm/2,
20*mm/2,
0.5*(fSiliconPlate_h+.1*mm));
G4VSolid* solid_LowerPlate_CutRect4 = new G4Box(solidname,
6*mm/2,
20*mm/2,
0.5*(fSiliconPlate_h+.1*mm));
// cut final silicon plate
G4VSolid* solid_LowerPlate_Base0 = new G4SubtractionSolid(solidname,
solid_LowerPlate_Rect,
solid_LowerPlate_CutRect0,
0,
G4ThreeVector((-.5*rect_bc_x-.5*.5*rect_x-.5*1*mm),0.,0.));
G4VSolid* solid_LowerPlate_Base1 = new G4SubtractionSolid(solidname,
solid_LowerPlate_Base0,
solid_LowerPlate_CutRect1,
rm,
G4ThreeVector(rect_x/2,trg_h,0.));
G4VSolid* solid_LowerPlate_Base2 = new G4SubtractionSolid(solidname,
solid_LowerPlate_Base1,
solid_LowerPlate_CutRect1,
rm,
G4ThreeVector(rect_x/2,-trg_h,0.));
G4VSolid* solid_LowerPlate_Base3 = new G4SubtractionSolid(solidname,
solid_LowerPlate_Base2,
solid_LowerPlate_CutRect2,
0,
G4ThreeVector(0.,(0.5*cut_y-3.5*mm),0.));
G4VSolid* solid_LowerPlate_Base4 = new G4SubtractionSolid(solidname,
solid_LowerPlate_Base3,
solid_LowerPlate_CutRect3,
0,
G4ThreeVector(rect_x/2,0,0.));
G4VSolid* solid_LowerPlate_Base5 = new G4SubtractionSolid(solidname,
solid_LowerPlate_Base4,
solid_LowerPlate_CutRect4,
rm1,
G4ThreeVector(-rect_bc_x,trg_h,0.));
G4VSolid* solid_LowerPlate_Base = new G4SubtractionSolid(solidname,
solid_LowerPlate_Base5,
solid_LowerPlate_CutRect4,
rm2,
G4ThreeVector(-rect_bc_x,-trg_h,0.));
G4Tubs* solid_LowerPlate_hole = new G4Tubs(solidname,
0*mm,
(1.5*mm+0.1*mm),
0.5*fSiliconPlate_h,
TubsStartAngle, TubsSpanningAngle);
//if(fCrystalRadius>41.*mm){fMaximalAllowedCrystalRadius=46.*mm;}
//radius condition identifies GTFs
G4double tmpR = 46.*mm + 1*mm + 1.5*mm;
G4VSolid* current_solid_LowerPlate = solid_LowerPlate_Base;
for (int i = 0; i < 3; i++) {
const G4double VertBarAngle = ((G4double) i) * 120.*deg;
const G4ThreeVector VertBarTranslation (/* x */ tmpR * std::cos(VertBarAngle), /* y */ tmpR * std::sin(VertBarAngle), /* z */ 0.);
G4VSolid* solid_LowerPlate_nHoles = new G4SubtractionSolid(solidname,
current_solid_LowerPlate,
solid_LowerPlate_hole,
0,
VertBarTranslation);
current_solid_LowerPlate = solid_LowerPlate_nHoles;
}
G4VSolid* final_plate = current_solid_LowerPlate;
return final_plate;
}
This diff is collapsed.
#include "SteppingAction.hh"
#include "EventAction.hh"
#include "G4SteppingManager.hh"
#include "G4SDManager.hh"
#include "G4EventManager.hh"
#include "G4ProcessManager.hh"
#include "G4Track.hh"
#include "G4Step.hh"
#include "G4Event.hh"
#include "G4StepPoint.hh"
#include "G4TrackStatus.hh"
#include "G4VPhysicalVolume.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTypes.hh"
#include "G4SystemOfUnits.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
SteppingAction::SteppingAction(EventAction* eventAction)
: fOneStepPrimaries(false), fEventAction(eventAction)
{
fExpectedNextStatus = Undefined;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
SteppingAction::~SteppingAction() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void SteppingAction::UserSteppingAction(const G4Step * theStep)
{
//auto analysisManager = G4AnalysisManager::Instance();
//G4Track* theTrack = theStep->GetTrack();
//G4ParticleDefinition* particleType = theTrack->GetDefinition();
//G4cout << "Begin Stepping" << G4endl;
fExpectedNextStatus = Undefined;
G4StepPoint* thePrePoint = theStep->GetPreStepPoint();
G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
//G4TouchableHistory* theTouchable = (G4TouchableHistory*)(thePrePoint->GetTouchable());
//
//G4StepStatus PreStepStatus = thePrePoint->GetStepStatus();
G4double edepStep = theStep->GetTotalEnergyDeposit()/eV;
//double PosDirz = step->GetTrack()->GetPosition().x();
//step->GetPreStepPoint()->GetPosition()
//
if ( thePrePV->GetName()=="target_1"){
fEventAction->AddInfo(thePrePoint->GetPosition().x(),thePrePoint->GetPosition().y(),thePrePoint->GetPosition().z(),edepStep);
}
G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
if(!thePostPV)
{//out of world
fExpectedNextStatus=Undefined;
return;
}
}
test.mac 0 → 100644
#just in case, it shoud be true by default
/process/em/fluo true
/process/em/auger true
/process/em/augerCascade true
/process/em/pixe true
#/run/numberOfThreads 4
#/Ge2SOS/det/setOutputDirectory /nfs/tegile/work/experiences/detecteurs/manzanillas/LEAPS_INNOV/Geant4_output/
/Ge2SOS/det/setDetectorType 5
/Ge2SOS/det/setSetupName run_1
#select output format, options are: csv root hdf5
/Ge2SOS/det/setDataType csv
/Ge2SOS/det/setNTargetSamples 1
/Ge2SOS/det/setGeDetectorThickness 6.99 mm
/Ge2SOS/det/setGeDetectorLength 36. mm
/Ge2SOS/det/setGeDetectorWidth 36. mm
/Ge2SOS/det/setCollimatorThickness 1. mm
/Ge2SOS/det/setCollimatorMat titanium_grade1
/Ge2SOS/det/setContactThickness 300. nm
/Ge2SOS/det/setBeWindowRadius 26.0 mm
#Set the position of the collimator, a values between -37 and 37 mm (size of samples)
#Choise of source type: 0 gamma,1 Fe-55,2 Cs-137; 1 Bi-207; 2 Sr-90; 3 e-
/run/initialize
/Ge2SOS/gun/sourceType 1
/Ge2SOS/gun/sourceDiameter 0.1 mm
#/Ge2SOS/gun/sourceEnergy 60.0 keV
/run/beamOn 9000000
vis.mac 0 → 100644
# Use this open statement to create an OpenGL view:
/vis/open OGL 600x600-0+0
#
# Use this open statement to create a .prim file suitable for
# viewing in DAWN:
#/vis/open DAWNFILE
#
# Use this open statement to create a .heprep file suitable for
# viewing in HepRApp:
#/vis/open HepRepFile
#
# Use this open statement to create a .wrl file suitable for
# viewing in a VRML viewer:
#/vis/open VRML2FILE
#
/Ge2SOS/det/setDetectorType 0
/Ge2SOS/gun/sourceType 0
# Disable auto refresh and quieten vis messages whilst scene and
# trajectories are established:
/vis/viewer/set/autoRefresh false
#/vis/verbose errors
#
# Draw geometry:
/vis/drawVolume
#
/vis/viewer/set/background white
# Specify view angle:
#/vis/viewer/set/viewpointThetaPhi 45. 0.
#
# Specify zoom value:
/vis/viewer/zoom 1.0
#
# Specify style (surface or wireframe):
#/vis/viewer/set/style wireframe
#
# Draw coordinate axes:
#/vis/scene/add/axes 0 0 0 1 m
#
# Draw smooth trajectories at end of event, showing trajectory points
# as markers 2 pixels wide:
/vis/scene/add/trajectories smooth
/vis/modeling/trajectories/create/drawByCharge
/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true
/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 2
# (if too many tracks cause core dump => /tracking/storeTrajectory 0)
#
# Draw hits at end of event:
#/vis/scene/add/hits
#
# To draw only gammas:
#/vis/filtering/trajectories/create/particleFilter
#/vis/filtering/trajectories/particleFilter-0/add gamma
#
# To invert the above, drawing all particles except gammas,
# keep the above two lines but also add:
#/vis/filtering/trajectories/particleFilter-0/invert true
#
# Many other options are available with /vis/modeling and /vis/filtering.
# For example, to select colour by particle ID:
/vis/modeling/trajectories/create/drawByParticleID
/vis/modeling/trajectories/drawByParticleID-0/set gamma green
#
# To superimpose all of the events from a given run:
#/vis/scene/endOfEventAction accumulate
#
# Re-establish auto refreshing and verbosity:
/vis/viewer/set/autoRefresh true
/vis/verbose warnings
#
# For file-based drivers, use this to create an empty detector view:
#/vis/viewer/flush
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment