lars/minicalosim
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

snapshot with output

Browse Source
This commit is contained in:
Jan Kieseler
2023-09-08 14:33:33 +02:00
parent 0adb64f763
commit 90975a7c24
14 changed files with 347 additions and 143 deletions
Download Patch File Download Diff File Expand all files Collapse all files
exampleB4a.cc
+11 -8
View File
@@ -39,6 +39,8 @@
#include "FTFP_BERT.hh"
#include "Randomize.hh"
#include "ConstructionWrapper.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace {
@@ -112,18 +114,18 @@ int main(int argc,char** argv)
}
#endif
ConstructionWrapper cw;
cw.addLayer(1, "G4_Pb", false, 2, 2);
cw.addLayer(10, "G4_Si", true, 2, 2);
cw.addLayer(1, "G4_Pb", false, 2, 2);
cw.addLayer(20, "G4_Si", true, 2, 2);
cw.addLayer(4, "G4_Pb", false, 2, 2);
cw.addLayer(3, "G4_Si", true, 2, 2);
ConstructionWrapper * cw = new ConstructionWrapper();
cw->addLayer(1, "G4_Pb", false);
cw->addLayer(10, "G4_Si", true, 9);
cw->addLayer(1, "G4_Pb", false);
cw->addLayer(20, "G4_Si", true, 15);
cw->addLayer(4, "G4_Pb", false);
cw->addLayer(3, "G4_Si", true, 23);
// Set mandatory initialization classes
//
auto detConstruction = new B4::DetectorConstruction(cw);
// FOR LATER: assign constructionwrapper here (will be passed from other function)
runManager->SetUserInitialization(detConstruction);
@@ -167,6 +169,7 @@ int main(int argc,char** argv)
delete visManager;
delete runManager;
delete cw;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
include/ConstructionWrapper.hh
+54 -3
View File
@@ -1,7 +1,34 @@
#ifndef CONSTRUCTIONWRAPPER_HH
#define CONSTRUCTIONWRAPPER_HH
#include "G4ThreeVector.hh"
#include <vector>
#include <string>
#include "G4VPhysicalVolume.hh"
//#include "G4VPhysicalVolume.hh"
class G4VPhysicalVolume;
class Sensor{
public:
Sensor(){};
~Sensor(){};
const G4double getEnergy()const{
return energy;
}
G4ThreeVector getPos()const{
return position;
}
G4ThreeVector getSize()const{
return size;
}
G4ThreeVector position;
G4ThreeVector size;
mutable G4double energy;
};
class Layer{
@@ -25,6 +52,8 @@ public:
bool isActive;
G4VPhysicalVolume* physicalVolume;
std::string name;
std::vector<Sensor> sensors;
};
@@ -34,7 +63,7 @@ public:
ConstructionWrapper(double xy_width=50);
~ConstructionWrapper() {};
void addLayer(double thickness_cm, std::string material, bool isActive=true, int nx=1, int ny=1);
void addLayer(double thickness_cm, std::string material, bool isActive=true, int nx=1, int ny=-1);
std::vector<Layer>& getLayers() ;
const std::vector<Layer>& getLayers() const;
@@ -42,7 +71,29 @@ public:
double getXYWidth() const{
return xywidth;
}
void resetSensorEnergies()const;//energies are mutable
int getNSensors()const{
int n_sensors = 0;
for(const auto& layer: layers){
n_sensors += layer.sensors.size();
}
return n_sensors;
}
Layer* getLayerByVolume(G4VPhysicalVolume* volume);
const Layer* getLayerByVolume(G4VPhysicalVolume* volume)const;
Sensor* getSensorByVolume(G4VPhysicalVolume* volume);
const Sensor* getSensorByVolume(G4VPhysicalVolume* volume)const;
void printSensorEnergies()const;
private:
double xywidth;
std::vector<Layer> layers;
};
};
#endif
include/DetectorConstruction.hh
+10 -21
View File
@@ -32,10 +32,10 @@
#include "G4VUserDetectorConstruction.hh"
#include "globals.hh"
#include "ConstructionWrapper.hh"
class G4VPhysicalVolume;
class G4GlobalMagFieldMessenger;
class ConstructionWrapper;
namespace B4
{
@@ -59,19 +59,20 @@ namespace B4
class DetectorConstruction : public G4VUserDetectorConstruction
{
public:
DetectorConstruction(const ConstructionWrapper& cw){
this->cw = cw;
}
DetectorConstruction( ConstructionWrapper* Cw): G4VUserDetectorConstruction(), cw(Cw){}
~DetectorConstruction() override = default;
public:
G4VPhysicalVolume* Construct() override;
void ConstructSDandField() override;
// get methods
//
const G4VPhysicalVolume* GetAbsorberPV() const;
const G4VPhysicalVolume* GetGapPV() const;
ConstructionWrapper* getConstructionWrapper(){
return cw;
}
const ConstructionWrapper* getConstructionWrapper()const{
return cw;
}
private:
// methods
@@ -84,23 +85,11 @@ class DetectorConstruction : public G4VUserDetectorConstruction
static G4ThreadLocal G4GlobalMagFieldMessenger* fMagFieldMessenger;
// magnetic field messenger
G4VPhysicalVolume* fAbsorberPV = nullptr; // the absorber physical volume
G4VPhysicalVolume* fGapPV = nullptr; // the gap physical volume
G4bool fCheckOverlaps = true; // option to activate checking of volumes overlaps
ConstructionWrapper cw;
ConstructionWrapper* cw=nullptr;
};
// inline functions
inline const G4VPhysicalVolume* DetectorConstruction::GetAbsorberPV() const {
return fAbsorberPV;
}
inline const G4VPhysicalVolume* DetectorConstruction::GetGapPV() const {
return fGapPV;
}
}
include/EventAction.hh
+16 -17
View File
@@ -31,7 +31,10 @@
#define B4aEventAction_h 1
#include "G4UserEventAction.hh"
#include "RunAction.hh"
#include "PrimaryGeneratorAction.hh"
#include "globals.hh"
class ConstructionWrapper;
namespace B4a
{
@@ -51,30 +54,26 @@ class EventAction : public G4UserEventAction
~EventAction() override = default;
void BeginOfEventAction(const G4Event* event) override;
void EndOfEventAction(const G4Event* event) override;
void EndOfEventAction(const G4Event* event) override;
void AddAbs(G4double de, G4double dl);
void AddGap(G4double de, G4double dl);
void setConstructionWrapper(const ConstructionWrapper * Cw){
this->cw = Cw;
}
void setPrimaryGeneratorAction(const B4::PrimaryGeneratorAction * Gen){
this->gen = Gen;
}
void setRunAction(const B4::RunAction * Run){
this->run = Run;
}
private:
G4double fEnergyAbs = 0.;
G4double fEnergyGap = 0.;
G4double fTrackLAbs = 0.;
G4double fTrackLGap = 0.;
const ConstructionWrapper * cw=nullptr;
const B4::PrimaryGeneratorAction * gen=nullptr;
const B4::RunAction * run=nullptr;
};
// inline functions
inline void EventAction::AddAbs(G4double de, G4double dl) {
fEnergyAbs += de;
fTrackLAbs += dl;
}
inline void EventAction::AddGap(G4double de, G4double dl) {
fEnergyGap += de;
fTrackLGap += dl;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
include/PrimaryGeneratorAction.hh
+9 -1
View File
@@ -56,9 +56,17 @@ public:
// set methods
void SetRandomFlag(G4bool value);
G4double getPartEnergy() const{return partEnergy;}
void setMinPartEnergy(G4double value){minPartEnergy = value;}
void setMaxPartEnergy(G4double value){maxPartEnergy = value;}
void setPartSpecies(G4String value){partSpecies = value;}
private:
G4ParticleGun* fParticleGun = nullptr; // G4 particle gun
G4double partEnergy;
G4double minPartEnergy; //MeV
G4double maxPartEnergy;
G4String partSpecies;
};
}
include/RunAction.hh
+14
View File
@@ -64,6 +64,20 @@ class RunAction : public G4UserRunAction
void BeginOfRunAction(const G4Run*) override;
void EndOfRunAction(const G4Run*) override;
//now this is dumb but the way this is setup it has to be done like this
//yeah mutable fine whatever
mutable std::vector<double> hitEnergy;
mutable std::vector<double> hitX;
mutable std::vector<double> hitY;
mutable std::vector<double> hitZ;
mutable std::vector<double> hitDX;
mutable std::vector<double> hitDY;
mutable std::vector<double> hitDZ;
mutable std::vector<int> hitLayer;
mutable std::vector<int> hitCopyNumber;
};
}
include/SteppingAction.hh
+2
View File
@@ -31,6 +31,7 @@
#define B4aSteppingAction_h 1
#include "G4UserSteppingAction.hh"
#include "ConstructionWrapper.hh"
namespace B4
{
@@ -57,6 +58,7 @@ public:
void UserSteppingAction(const G4Step* step) override;
private:
const B4::DetectorConstruction* fDetConstruction = nullptr;
EventAction* fEventAction = nullptr;
src/ActionInitialization.cc
+10 -3
View File
@@ -52,15 +52,22 @@ void ActionInitialization::BuildForMaster() const
SetUserAction(new RunAction);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ActionInitialization::Build() const
{
SetUserAction(new PrimaryGeneratorAction);
SetUserAction(new RunAction);
auto gen = new PrimaryGeneratorAction;
SetUserAction(gen);
auto runact = new RunAction;
SetUserAction(runact);
auto eventAction = new EventAction;
eventAction->setConstructionWrapper(fDetConstruction->getConstructionWrapper());
eventAction->setRunAction(runact);
eventAction->setPrimaryGeneratorAction(gen);
SetUserAction(eventAction);
SetUserAction(new SteppingAction(fDetConstruction,eventAction));
auto steppingAction = new SteppingAction(fDetConstruction,eventAction);
SetUserAction(steppingAction);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
src/ConstructionWrapper.cc
+67
View File
@@ -1,4 +1,7 @@
#include "ConstructionWrapper.hh"
#include "G4VPhysicalVolume.hh"
Layer::Layer(){
thickness = 0;
@@ -38,6 +41,10 @@ ConstructionWrapper::ConstructionWrapper(double xy_width){
}
void ConstructionWrapper::addLayer(double thickness, std::string material, bool isActive, int nx, int ny){
if(ny<0){
ny = nx;
}
Layer layer;
layer.setThickness(thickness);
layer.setMaterial(material);
@@ -60,3 +67,63 @@ const std::vector<Layer> & ConstructionWrapper::getLayers() const{
std::vector<Layer> & ConstructionWrapper::getLayers(){
return layers;
}
void ConstructionWrapper::resetSensorEnergies()const {
for(const auto & layer : layers){
for(const auto & sensor : layer.sensors){
sensor.energy = 0;
}
}
}
Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume){
for(auto & layer : layers){
if(layer.physicalVolume == volume){
return &layer;
}
}
return nullptr;
}
const Layer* ConstructionWrapper::getLayerByVolume(G4VPhysicalVolume* volume)const{
for(auto & layer : layers){
if(layer.physicalVolume == volume){
return &layer;
}
}
return nullptr;
}
Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume){
int copyNo = volume->GetCopyNo();
auto layer = getLayerByVolume(volume);
if(layer == nullptr){
return nullptr;
}
if(copyNo >= layer->sensors.size()){
return nullptr;
}
return &layer->sensors[copyNo];
}
const Sensor* ConstructionWrapper::getSensorByVolume(G4VPhysicalVolume* volume)const{
int copyNo = volume->GetCopyNo();
auto layer = getLayerByVolume(volume);
if(layer == nullptr){
return nullptr;
}
if(copyNo >= layer->sensors.size()){
return nullptr;
}
return &layer->sensors[copyNo];
}
void ConstructionWrapper::printSensorEnergies()const{
for(const auto& layer: layers){
for(const auto& sensor: layer.sensors){
G4cout << sensor.energy << " ";
}
G4cout << G4endl;
}
}
src/DetectorConstruction.cc
+36 -5
View File
@@ -28,6 +28,7 @@
/// \brief Implementation of the B4::DetectorConstruction class
#include "DetectorConstruction.hh"
#include "ConstructionWrapper.hh"
#include "G4Material.hh"
#include "G4NistManager.hh"
@@ -100,6 +101,7 @@ G4GlobalMagFieldMessenger* DetectorConstruction::fMagFieldMessenger = nullptr;
G4VPhysicalVolume* DetectorConstruction::Construct()
{
fCheckOverlaps=true;
// Define materials
DefineMaterials();
@@ -114,7 +116,7 @@ void DetectorConstruction::DefineMaterials()
// Lead material defined using NIST Manager
auto nistManager = G4NistManager::Instance();
nistManager->FindOrBuildMaterial("G4_AIR");
auto cwLayers = cw.getLayers();
auto cwLayers = cw->getLayers();
for(auto layer : cwLayers){
nistManager->FindOrBuildMaterial(layer.material);
}
@@ -128,14 +130,14 @@ nistManager->FindOrBuildMaterial("G4_AIR");
G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
{
// Geometry parameters
auto & cwLayers = cw.getLayers();
auto & cwLayers = cw->getLayers();
G4int nofLayers = cwLayers.size();
G4double caloLength = 0;
for(auto layer : cwLayers){
caloLength += layer.thickness * cm;
}
G4double calorSizeXY = cw.getXYWidth() * cm;
G4double calorSizeXY = cw->getXYWidth() * cm;
auto worldSizeXY = 1.2 * calorSizeXY;
auto worldSizeZ = 1.2 * caloLength;
@@ -210,6 +212,21 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
defaultMaterial, //G4Material::GetMaterial(layer.material), // its material
layer.name); // its name
//set layerLV to be visible and green if active and transparent otherwise
if(layer.isActive){
auto layerVisAtt = new G4VisAttributes(G4Colour(0.0,1.0,0.0));
layerVisAtt->SetVisibility(true);
layerVisAtt->SetForceSolid(true);
layerLV->SetVisAttributes(layerVisAtt);
} //else set the absorbers to be transparent
else{
auto layerVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.1));
layerVisAtt->SetVisibility(true);
layerVisAtt->SetForceSolid(true);
layerLV->SetVisAttributes(layerVisAtt);
}
auto sensorLV = new G4LogicalVolume(
new G4Box("sensor", layer.sens_xwidth/2*cm, layer.sens_ywidth/2*cm, layer.thickness/2*cm),
G4Material::GetMaterial(layer.material),
@@ -231,7 +248,20 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
0, // copy number
fCheckOverlaps); // checking overlaps
layer.assignPhysicalVolume(pv);
layer.assignPhysicalVolume(ppv); //maybe this needs to be ppv - check
// assign sensors to layer by copyNumber; access the copyNumber of the volumes
if(layer.isActive){
for(int i = 0; i < layer.nx*layer.ny; i++){
Sensor sensor;
sensor.position = G4ThreeVector(layer.sens_xwidth * (i % layer.nx) - layer.sens_xwidth * ((float)layer.nx-1) / 2.*cm,
layer.sens_ywidth * (i / layer.nx) - layer.sens_ywidth * ((float)layer.ny-1) / 2.*cm, position);
sensor.size = G4ThreeVector(layer.sens_xwidth, layer.sens_ywidth, layer.thickness);
sensor.energy = 0;
layer.sensors.push_back(sensor);//this should now be aligned with copy number
}
}
position += layer.thickness / 2 *cm; //assign the physical volume to the layer
}
//
@@ -239,10 +269,11 @@ G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
//
worldLV->SetVisAttributes (G4VisAttributes::GetInvisible());
auto simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0));
auto simpleBoxVisAtt= new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.));
simpleBoxVisAtt->SetVisibility(true);
calorLV->SetVisAttributes(simpleBoxVisAtt);
//
// Always return the physical World
//
src/EventAction.cc
+61 -31
View File
@@ -27,6 +27,8 @@
/// \file B4/B4a/src/EventAction.cc
/// \brief Implementation of the B4a::EventAction class
#include "ConstructionWrapper.hh"
#include "EventAction.hh"
#include "RunAction.hh"
@@ -34,6 +36,7 @@
#include "G4RunManager.hh"
#include "G4Event.hh"
#include "G4UnitsTable.hh"
#include "PrimaryGeneratorAction.hh"
#include "Randomize.hh"
#include <iomanip>
@@ -46,10 +49,28 @@ namespace B4a
void EventAction::BeginOfEventAction(const G4Event* /*event*/)
{
// initialisation per event
fEnergyAbs = 0.;
fEnergyGap = 0.;
fTrackLAbs = 0.;
fTrackLGap = 0.;
if(cw == nullptr){
G4cout << "ConstructionWrapper not found" << G4endl;
throw std::runtime_error("ConstructionWrapper not found");
}
if(gen == nullptr){
G4cout << "PrimaryGeneratorAction not found" << G4endl;
throw std::runtime_error("PrimaryGeneratorAction not found");
}
cw->resetSensorEnergies();
if(run->hitEnergy.size() != cw->getNSensors()){
run->hitEnergy.resize(cw->getNSensors(),0);
run->hitX.resize(cw->getNSensors(),0);
run->hitY.resize(cw->getNSensors(),0);
run->hitZ.resize(cw->getNSensors(),0);
run->hitDX.resize(cw->getNSensors(),0);
run->hitDY.resize(cw->getNSensors(),0);
run->hitDZ.resize(cw->getNSensors(),0);
run->hitLayer.resize(cw->getNSensors(),0);
run->hitCopyNumber.resize(cw->getNSensors(),0);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -62,38 +83,47 @@ void EventAction::EndOfEventAction(const G4Event* event)
// get analysis manager
auto analysisManager = G4AnalysisManager::Instance();
// fill histograms
analysisManager->FillH1(0, fEnergyAbs);
analysisManager->FillH1(1, fEnergyGap);
analysisManager->FillH1(2, fTrackLAbs);
analysisManager->FillH1(3, fTrackLGap);
//cw->printSensorEnergies();//DEBUG
std::vector<double> sensorEnergies;
double tot_energy=0;
int sensorNumber=0;
int layerNumber=0;
for(const auto& layer: cw->getLayers()){
for(const auto& sensor: layer.sensors){
sensorEnergies.push_back(sensor.getEnergy());
tot_energy+= sensor.getEnergy();
//fill the vectors
run->hitEnergy[sensorNumber] = sensor.getEnergy();
run->hitX[sensorNumber] = sensor.getPos().x();
run->hitY[sensorNumber] = sensor.getPos().y();
run->hitZ[sensorNumber] = sensor.getPos().z();
run->hitDX[sensorNumber] = sensor.getSize().x();
run->hitDY[sensorNumber] = sensor.getSize().y();
run->hitDZ[sensorNumber] = sensor.getSize().z();
run->hitLayer[sensorNumber] = layerNumber;
run->hitCopyNumber[sensorNumber] = sensorNumber;
sensorNumber++;
}
layerNumber++;
}
//assign the vectors
// fill ntuple
analysisManager->FillNtupleDColumn(0, fEnergyAbs);
analysisManager->FillNtupleDColumn(1, fEnergyGap);
analysisManager->FillNtupleDColumn(2, fTrackLAbs);
analysisManager->FillNtupleDColumn(3, fTrackLGap);
analysisManager->FillNtupleDColumn(0, gen->getPartEnergy()); //needs to be true energy DEBUG
analysisManager->FillNtupleDColumn(1, tot_energy);
analysisManager->FillNtupleIColumn(2, (int)cw->getLayers().size());
analysisManager->FillNtupleIColumn(3, cw->getNSensors());
//vectors are stored automatically
analysisManager->AddNtupleRow();
// Print per event (modulo n)
//
auto eventID = event->GetEventID();
auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
if ( ( printModulo > 0 ) && ( eventID % printModulo == 0 ) ) {
G4cout << "---> End of event: " << eventID << G4endl;
G4cout
<< " Absorber: total energy: " << std::setw(7)
<< G4BestUnit(fEnergyAbs,"Energy")
<< " total track length: " << std::setw(7)
<< G4BestUnit(fTrackLAbs,"Length")
<< G4endl
<< " Gap: total energy: " << std::setw(7)
<< G4BestUnit(fEnergyGap,"Energy")
<< " total track length: " << std::setw(7)
<< G4BestUnit(fTrackLGap,"Length")
<< G4endl;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
src/PrimaryGeneratorAction.cc
+21 -5
View File
@@ -39,6 +39,8 @@
#include "G4ParticleDefinition.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "G4INCLRandom.hh"
#include <G4INCLGeant4Random.hh>
namespace B4
{
@@ -50,13 +52,15 @@ PrimaryGeneratorAction::PrimaryGeneratorAction()
G4int nofParticles = 1;
fParticleGun = new G4ParticleGun(nofParticles);
partSpecies = "e-";
minPartEnergy = 1000;
maxPartEnergy = 1000;
G4INCL::Random::setGenerator( new G4INCL::Geant4RandomGenerator());
// default particle kinematic
//
auto particleDefinition
= G4ParticleTable::GetParticleTable()->FindParticle("e-");
fParticleGun->SetParticleDefinition(particleDefinition);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
fParticleGun->SetParticleEnergy(50.*MeV);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -70,8 +74,19 @@ PrimaryGeneratorAction::~PrimaryGeneratorAction()
void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
// This function is called at the begining of event
G4double rand = G4INCL::Random::shoot();
partEnergy = rand*(maxPartEnergy-minPartEnergy)+minPartEnergy;
auto particleDefinition
= G4ParticleTable::GetParticleTable()->FindParticle(partSpecies);
fParticleGun->SetParticleDefinition(particleDefinition);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
fParticleGun->SetParticleEnergy(partEnergy*MeV);
// In order to avoid dependence of PrimaryGeneratorAction
// on DetectorConstruction class we get world volume
// from G4LogicalVolumeStore
@@ -102,6 +117,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
->SetParticlePosition(G4ThreeVector(0., 0., -worldZHalfLength));
fParticleGun->GeneratePrimaryVertex(anEvent);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
src/RunAction.cc
+19 -40
View File
@@ -60,19 +60,26 @@ RunAction::RunAction()
// Book histograms, ntuple
//
// Creating histograms
analysisManager->CreateH1("Eabs","Edep in absorber", 100, 0., 800*MeV);
analysisManager->CreateH1("Egap","Edep in gap", 100, 0., 100*MeV);
analysisManager->CreateH1("Labs","trackL in absorber", 100, 0., 1*m);
analysisManager->CreateH1("Lgap","trackL in gap", 100, 0., 50*cm);
// Creating ntuple
//
analysisManager->CreateNtuple("B4", "Edep and TrackL");
analysisManager->CreateNtupleDColumn("Eabs");
analysisManager->CreateNtupleDColumn("Egap");
analysisManager->CreateNtupleDColumn("Labs");
analysisManager->CreateNtupleDColumn("Lgap");
analysisManager->CreateNtuple("Hits", "Hits");
analysisManager->CreateNtupleDColumn("true_energy");
analysisManager->CreateNtupleDColumn("total_dep_energy");
analysisManager->CreateNtupleIColumn("N_layers");
analysisManager->CreateNtupleIColumn("N_sensors");
analysisManager->CreateNtupleDColumn("sensor_energy",hitEnergy);
analysisManager->CreateNtupleDColumn("sensor_x",hitX);
analysisManager->CreateNtupleDColumn("sensor_y",hitY);
analysisManager->CreateNtupleDColumn("sensor_z",hitZ);
analysisManager->CreateNtupleDColumn("sensor_dx",hitDX);
analysisManager->CreateNtupleDColumn("sensor_dy",hitDY);
analysisManager->CreateNtupleDColumn("sensor_dz",hitDZ);
analysisManager->CreateNtupleIColumn("sensor_layer",hitLayer);
analysisManager->CreateNtupleIColumn("sensor_copy_number",hitCopyNumber);
analysisManager->FinishNtuple();
}
@@ -104,35 +111,7 @@ void RunAction::EndOfRunAction(const G4Run* /*run*/)
// print histogram statistics
//
auto analysisManager = G4AnalysisManager::Instance();
if ( analysisManager->GetH1(1) ) {
G4cout << G4endl << " ----> print histograms statistic ";
if(isMaster) {
G4cout << "for the entire run " << G4endl << G4endl;
}
else {
G4cout << "for the local thread " << G4endl << G4endl;
}
G4cout << " EAbs : mean = "
<< G4BestUnit(analysisManager->GetH1(0)->mean(), "Energy")
<< " rms = "
<< G4BestUnit(analysisManager->GetH1(0)->rms(), "Energy") << G4endl;
G4cout << " EGap : mean = "
<< G4BestUnit(analysisManager->GetH1(1)->mean(), "Energy")
<< " rms = "
<< G4BestUnit(analysisManager->GetH1(1)->rms(), "Energy") << G4endl;
G4cout << " LAbs : mean = "
<< G4BestUnit(analysisManager->GetH1(2)->mean(), "Length")
<< " rms = "
<< G4BestUnit(analysisManager->GetH1(2)->rms(), "Length") << G4endl;
G4cout << " LGap : mean = "
<< G4BestUnit(analysisManager->GetH1(3)->mean(), "Length")
<< " rms = "
<< G4BestUnit(analysisManager->GetH1(3)->rms(), "Length") << G4endl;
}
// save histograms & ntuple
//
src/SteppingAction.cc
+17 -9
View File
@@ -59,19 +59,27 @@ void SteppingAction::UserSteppingAction(const G4Step* step)
// energy deposit
auto edep = step->GetTotalEnergyDeposit();
// step length
G4double stepLength = 0.;
if ( step->GetTrack()->GetDefinition()->GetPDGCharge() != 0. ) {
stepLength = step->GetStepLength();
//return;
//find the layer through the constructionwrapper that corresponds to the volume
auto cw = fDetConstruction->getConstructionWrapper();
if(cw == nullptr){
G4cout << "ConstructionWrapper not found" << G4endl;
throw std::runtime_error("ConstructionWrapper not found");
}
if ( volume == fDetConstruction->GetAbsorberPV() ) {
fEventAction->AddAbs(edep,stepLength);
auto sensor = fDetConstruction->getConstructionWrapper()->getSensorByVolume(volume);
if(sensor == nullptr){
//can simply be not an active volume, so no need to throw an error
//G4cout << "Sensor not found in" << volume->GetName() << G4endl; //DEBUG
return;
}
//G4cout << "Sensor found in " << volume->GetName() << " with copy number " << volume->GetCopyNo()<< G4endl; //DEBUG
sensor->energy += edep;
//find sensor by volume copy number
if ( volume == fDetConstruction->GetGapPV() ) {
fEventAction->AddGap(edep,stepLength);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......