From 7410049a722adafe18c233864b3c4a81fb43b0cb Mon Sep 17 00:00:00 2001 From: Jan Kieseler Date: Wed, 9 Oct 2024 13:05:18 +0200 Subject: [PATCH] added multi particles --- bind/G4Calo.py | 96 ++++++++++++++--------- exampleB4a.cc | 3 +- include/ActionInitialization.hh | 16 +++- include/G4System.hh | 3 +- include/PrimaryGeneratorAction.hh | 57 +++++++------- src/G4System.cc | 8 +- src/PrimaryGeneratorAction.cc | 126 ++++++++++++++++++------------ 7 files changed, 184 insertions(+), 125 deletions(-) diff --git a/bind/G4Calo.py b/bind/G4Calo.py index afafd23..bb808e8 100644 --- a/bind/G4Calo.py +++ b/bind/G4Calo.py @@ -32,7 +32,7 @@ class __G4System(_G4System): def run_batch( self, nEvents: int, - particleSpec: str, + particleSpec , minEnergy_GeV: float, maxEnergy_GeV: float = -1.0, filename: str = "", @@ -41,6 +41,17 @@ class __G4System(_G4System): if maxEnergy_GeV < 0: maxEnergy_GeV = minEnergy_GeV + + #check if particleSpec is a list of strings or a single string and assert + if not isinstance(particleSpec, str): + assert isinstance(particleSpec, list), "particleSpec must be a string or a list of strings" + for p in particleSpec: + assert isinstance(p, str), "particleSpec must be a string or a list of strings" + + + if isinstance(particleSpec, str): + particleSpec = [particleSpec] + save_file = len(filename) > 0 # filename without file ending(!) filename = "._" + str(time.perf_counter_ns()) + ".root" @@ -141,7 +152,7 @@ def index_out_of_bounds_workaround(tbranch): def _run_mini_batch( cw : GeometryDescriptor, nEvents: int, - particleSpec: str, + particleSpec , minEnergy_GeV: float, maxEnergy_GeV: float = -1.0, batch_seed : int = 0, @@ -173,7 +184,7 @@ def init_g4calo_threading_lock(): def run_batch( gd : GeometryDescriptor, nEvents: int, - particleSpec: str, + particleSpec , minEnergy_GeV: float, maxEnergy_GeV: float = -1.0, filename: str = "", @@ -184,6 +195,11 @@ def run_batch( assert nEvents > 0 assert minEnergy_GeV > 0 + if not isinstance(particleSpec, str): + assert isinstance(particleSpec, list), "particleSpec must be a string or a list of strings" + for p in particleSpec: + assert isinstance(p, str), "particleSpec must be a string or a list of strings" + nCores = multiprocessing.cpu_count() #make sure to adjust cores such that at least 80 events are run per core nCores = min(nCores, nEvents // 80 + 1) @@ -230,7 +246,7 @@ def run_batch( return alldf def _fill_event(gd : GeometryDescriptor, - particleSpec: str, + particleSpec, energy: float, seed: int = -1): @@ -241,27 +257,34 @@ def _fill_event(gd : GeometryDescriptor, return gd def display_event(gd : GeometryDescriptor, - particleSpec: str, + particleSpec, energy: float, - logE = False, renderer=None, seed = -1): + logE = False, renderer=None, seed = -1, outfile : str = ""): #run _fill_event in forked mode using 1-core multiprocessing to avoid G4 singletons to interfere with multiprocessing.Pool(1) as pool: gd = pool.apply(_fill_event, (gd, particleSpec, energy)) + print('creating plot') #use gd to plot # for loop over all layers to_plot = [] material_dict = {} z0=0 - # sum up total deposited energy + # sum up total deposited energy, this is just a normalization factor total_dep_energy = 0 + max_dep_energy = 0 for layer in gd.getLayers(): for sensor in layer.sensors: total_dep_energy += sensor.getEnergy() + max_dep_energy = max(max_dep_energy, sensor.getEnergy()) if total_dep_energy == 0: print("No energy deposited in calorimeter!") total_dep_energy = 10**-8 # to avoid division by zero + if max_dep_energy == 0: + print("No energy deposited in any sensor!") + max_dep_energy = 10**-8 + # loop over layers, invert order for layer in gd.getLayers()[::-1]: @@ -282,7 +305,7 @@ def display_event(gd : GeometryDescriptor, 'color': col_dict[layer_material], 'showlegend': False, 'flatshading': True, - 'opacity': 0.2 * 20. / len(gd.getLayers())} + 'opacity': min(0.2 * 20. / len(gd.getLayers()) + 1e-3, 0.2)} # add legend entry to_plot.append(go.Mesh3d(x=[None], y=[None], z=[None], i=[0], j=[0], k=[0], color=material_dict[layer_material]['color'], @@ -310,7 +333,7 @@ def display_event(gd : GeometryDescriptor, y_center = sensor.getY() - corr hwidth = sensor.getdx() /2. energy = sensor.getEnergy() - use_energy = float(energy / total_dep_energy) + use_energy = float(energy / max_dep_energy) if logE: raise NotImplementedError use_energy = np.log(use_energy+1.) # - np.log(total_dep_energy) @@ -343,7 +366,7 @@ def display_event(gd : GeometryDescriptor, cmin=0, cmax=1, colorbar=dict( - title='Fraction of total deposited Energy', + title='Energy/max(Energy)', tickvals=[0, 1], ticktext=['0', '1'], ticks='outside', @@ -366,26 +389,27 @@ def display_event(gd : GeometryDescriptor, name='Incoming particle', showlegend=True, ) - # Calculate the direction vector for the arrow - direction_vector = [(end_point[0] - start_point[0]), (end_point[1] - start_point[1]), (end_point[2] - start_point[2])] - # Create the arrowhead at the start point with the opposite direction - arrowhead_trace = go.Cone( - x=[end_point[0]], - z=[end_point[1]], - y=[end_point[2]], - u=[direction_vector[0]], - w=[direction_vector[1]], - v=[direction_vector[2]], - sizemode='scaled', - sizeref=0.8, - showscale=False, - colorscale='Reds', - opacity=1.0, - anchor='tail', - ) - # Create the 3D scatter plot with both traces - to_plot.append(line_trace) - to_plot.append(arrowhead_trace) + if (not isinstance(particleSpec, list)) or len(particleSpec) == 1: + # Calculate the direction vector for the arrow + direction_vector = [(end_point[0] - start_point[0]), (end_point[1] - start_point[1]), (end_point[2] - start_point[2])] + # Create the arrowhead at the start point with the opposite direction + arrowhead_trace = go.Cone( + x=[end_point[0]], + z=[end_point[1]], + y=[end_point[2]], + u=[direction_vector[0]], + w=[direction_vector[1]], + v=[direction_vector[2]], + sizemode='scaled', + sizeref=0.8, + showscale=False, + colorscale='Reds', + opacity=1.0, + anchor='tail', + ) + # Create the 3D scatter plot with both traces + to_plot.append(line_trace) + to_plot.append(arrowhead_trace) # # finally show plot # @@ -403,12 +427,14 @@ def display_event(gd : GeometryDescriptor, #name the axes in the HEP way, so y and z switch names fig.update_layout(scene=dict(xaxis_title='x [mm]', yaxis_title='z [mm]', zaxis_title='y [mm]')) - #fig.write_html("temp.html") - #exit() - if renderer is not None: - fig.show(renderer=renderer) + if len(outfile) > 0: + fig.write_html(outfile) + else: - fig.show() + if renderer is not None: + fig.show(renderer=renderer) + else: + fig.show() diff --git a/exampleB4a.cc b/exampleB4a.cc index c9f4cf8..55351f9 100644 --- a/exampleB4a.cc +++ b/exampleB4a.cc @@ -273,7 +273,8 @@ int main(){ G4System builder(true); builder.init(cw); //builder.run_gui(); - builder.run_batch(10000, ((std::string)"e-").data(), 1, 100); + std::vector parts = {"e-"}; + builder.run_batch(10000,parts , 1, 100); //run(cw, 1000, "e-", 1, 100, true); diff --git a/include/ActionInitialization.hh b/include/ActionInitialization.hh index 0812553..5600786 100644 --- a/include/ActionInitialization.hh +++ b/include/ActionInitialization.hh @@ -53,15 +53,25 @@ class ActionInitialization : public G4VUserActionInitialization void BuildForMaster() const override; void Build() const override; - void setGeneratorProperties(G4double minEnergy, G4double maxEnergy, G4String partSpecies){ + void setGeneratorProperties(G4double minEnergy, G4double maxEnergy, const std::vector& partSpecies){ if(gen==nullptr){ throw std::runtime_error("PrimaryGeneratorAction not found"); } gen->setMinPartEnergy(minEnergy); gen->setMaxPartEnergy(maxEnergy); gen->setPartSpecies(partSpecies); - G4cout << "ActionInitialization: Setting particle species to " << partSpecies << G4endl; + G4cout << "ActionInitialization: Setting particle species to " ; + for(const auto& p : partSpecies){ + G4cout << p << " "; + } + G4cout << G4endl; } + void setGeneratorProperties(G4double minEnergy, G4double maxEnergy, const std::string& partSpecies){ + setGeneratorProperties(minEnergy, maxEnergy, std::vector{partSpecies}); + } + + + void setDetectorConstruction(B4::DetectorConstruction* det){ fDetConstruction = det; @@ -80,7 +90,7 @@ class ActionInitialization : public G4VUserActionInitialization mutable B4::RunAction * runact = nullptr; G4double minEnergy=1000; G4double maxEnergy=1000; - G4String partSpecies="e-"; + std::vector partSpecies; G4String filename=""; }; diff --git a/include/G4System.hh b/include/G4System.hh index 2b305fe..c807aba 100644 --- a/include/G4System.hh +++ b/include/G4System.hh @@ -17,6 +17,7 @@ #include "FTFP_BERT.hh" #include "Randomize.hh" +#include class G4System{ friend class GeometryDescriptor; @@ -39,7 +40,7 @@ void init(GeometryDescriptor &cw, int seed=0); void run_visualize(const std::string& partSpecies, double minEnergy_GeV, double maxEnergy_GeV); //runs the whole gui if available void run_gui(); -void run_batch(int nEvents, const std::string& partSpecies, double minEnergy_GeV, +void run_batch(int nEvents, const std::vector& partSpecies, double minEnergy_GeV, double maxEnergy_GeV, std::string filename="_1234567890_Hits.root"); void applyUICommand(const std::string& command){ diff --git a/include/PrimaryGeneratorAction.hh b/include/PrimaryGeneratorAction.hh index 9c279ae..0366a7f 100644 --- a/include/PrimaryGeneratorAction.hh +++ b/include/PrimaryGeneratorAction.hh @@ -59,32 +59,35 @@ public: G4double getPartEnergy() const{return partEnergy;} void setMinPartEnergy(G4double value){minPartEnergy = value;} void setMaxPartEnergy(G4double value){maxPartEnergy = value;} - void setPartSpecies(G4String value){ - partSpecies = value; - if(partSpecies == "e-"){ - G4cout << "Particle species set to electron" << G4endl; - } - else if(partSpecies == "e+"){ - G4cout << "Particle species set to positron" << G4endl; - } - else if(partSpecies == "gamma"){ - G4cout << "Particle species set to gamma" << G4endl; - } - else if(partSpecies == "proton"){ - G4cout << "Particle species set to proton" << G4endl; - } - else if(partSpecies == "neutron"){ - G4cout << "Particle species set to neutron" << G4endl; - } - else if(partSpecies == "pi+"){ - G4cout << "Particle species set to charged pion" << G4endl; - } - else if(partSpecies == "pi-"){ - G4cout << "Particle species set to charged pion" << G4endl; - } - else{ - G4cout << "Particle species "<< partSpecies << " not recognized" << G4endl; - throw std::invalid_argument("Particle species not recognized"); + void setPartSpecies(const std::vector& values){ + partSpecies = values; + for(const auto value : values){ + G4String val = value; + if(val == "e-"){ + G4cout << "Particle species set to electron" << G4endl; + } + else if(val == "e+"){ + G4cout << "Particle species set to positron" << G4endl; + } + else if(val == "gamma"){ + G4cout << "Particle species set to gamma" << G4endl; + } + else if(val == "proton"){ + G4cout << "Particle species set to proton" << G4endl; + } + else if(val == "neutron"){ + G4cout << "Particle species set to neutron" << G4endl; + } + else if(val == "pi+"){ + G4cout << "Particle species set to charged pion" << G4endl; + } + else if(val == "pi-"){ + G4cout << "Particle species set to charged pion" << G4endl; + } + else{ + G4cout << "Particle species "<< val << " not recognized" << G4endl; + throw std::invalid_argument("Particle species not recognized"); + } } } private: @@ -92,7 +95,7 @@ private: G4double partEnergy; G4double minPartEnergy=100; //MeV G4double maxPartEnergy=100; - G4String partSpecies="e-"; + std::vector partSpecies; }; diff --git a/src/G4System.cc b/src/G4System.cc index 72ac38c..9da8f63 100644 --- a/src/G4System.cc +++ b/src/G4System.cc @@ -162,16 +162,12 @@ void G4System::run_gui(){ } } -void G4System::run_batch(int nEvents,const std::string& partSpecies, double minEnergy, +void G4System::run_batch(int nEvents, const std::vector& partSpecies, double minEnergy, double maxEnergy, std::string filename){ check(); - - - G4String partSpec = partSpecies; - G4cout << "running with particle species " << partSpec << G4endl; UImanager->ApplyCommand("/run/initialize"); UImanager->ApplyCommand("/vis/disable"); - actionInitialization->setGeneratorProperties(minEnergy, maxEnergy, partSpec); + actionInitialization->setGeneratorProperties(minEnergy, maxEnergy, partSpecies); actionInitialization->setFileName(filename); runManager->BeamOn(nEvents); } diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc index 463d115..e6832d7 100644 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -55,7 +55,7 @@ PrimaryGeneratorAction::PrimaryGeneratorAction() G4int nofParticles = 1; fParticleGun = new G4ParticleGun(nofParticles); - partSpecies = "e-"; + partSpecies = std::vector{"e-"}; minPartEnergy = 1000; maxPartEnergy = 1000; @@ -84,62 +84,84 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) rand = G4INCL::Random::shoot(); partEnergy = rand*(maxPartEnergy-minPartEnergy)+minPartEnergy; - //G4cout << "shooting with " << partEnergy << " GeV" << G4endl; - - auto particleDefinition - = G4ParticleTable::GetParticleTable()->FindParticle(partSpecies); - if(particleDefinition == nullptr){ - G4cout << "Particle " << partSpecies << " not found" << G4endl; - throw std::runtime_error("Particle not found"); - return; + //create fractions for each entry in the partSpecies vector such that the partEnergy can be distributed + //according to the fractions + std::vector part_energies; + G4double sum = 0; + for(const auto& part : partSpecies){ + G4double randfrac = 0; + while(randfrac > 1 || randfrac <= 0) + randfrac = G4INCL::Random::shoot(); + part_energies.push_back(randfrac); + sum += randfrac; } - fParticleGun->SetParticleDefinition(particleDefinition); - fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); - fParticleGun->SetParticleEnergy(partEnergy*GeV); - - // In order to avoid dependence of PrimaryGeneratorAction - // on DetectorConstruction class we get world volume - // from G4LogicalVolumeStore - // - G4double worldZHalfLength = 0.; - auto worldLV = G4LogicalVolumeStore::GetInstance()->GetVolume("World"); - - // Check that the world volume has box shape - G4Box* worldBox = nullptr; - if ( worldLV ) { - worldBox = dynamic_cast(worldLV->GetSolid()); + //normalise w.r.t total + for(auto& en : part_energies){ + en /= sum; + en *= partEnergy; } - if ( worldBox ) { - worldZHalfLength = worldBox->GetZHalfLength(); - } - else { - G4ExceptionDescription msg; - msg << "World volume of box shape not found." << G4endl; - msg << "Perhaps you have changed geometry." << G4endl; - msg << "The gun will be place in the center."; - G4Exception("PrimaryGeneratorAction::GeneratePrimaries()", - "MyCode0002", JustWarning, msg); - } - if(false){ - //draw particle position from a unit distribution in x-y with +-3 cm around the center - G4double x = -1000; - G4double y = -1000; - G4double x_width = 6; - G4double y_width = 6; - while(x < -x_width/2. || x > x_width/2. || y < -y_width/2. || y > y_width/2.){ - x = G4INCL::Random::shoot()*x_width - x_width/2; - y = G4INCL::Random::shoot()*y_width - y_width/2; + // now create a vertex for every particle + + for (size_t i = 0; i < partSpecies.size(); i++){ + G4String single_partType = partSpecies[i]; + G4double single_partEnergy = part_energies[i]; + //create a vertex for the particle + auto particleDefinition + = G4ParticleTable::GetParticleTable()->FindParticle(single_partType); + if(particleDefinition == nullptr){ + G4cout << "Particle " << single_partType << " not found" << G4endl; + throw std::runtime_error("Particle not found"); + return; } - fParticleGun->SetParticlePosition(G4ThreeVector(x*cm, y*cm, -worldZHalfLength)); + fParticleGun->SetParticleDefinition(particleDefinition); + fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); + fParticleGun->SetParticleEnergy(single_partEnergy*GeV); + + // In order to avoid dependence of PrimaryGeneratorAction + // on DetectorConstruction class we get world volume + // from G4LogicalVolumeStore + // + G4double worldZHalfLength = 0.; + auto worldLV = G4LogicalVolumeStore::GetInstance()->GetVolume("World"); + + // Check that the world volume has box shape + G4Box* worldBox = nullptr; + if ( worldLV ) { + worldBox = dynamic_cast(worldLV->GetSolid()); + } + + if ( worldBox ) { + worldZHalfLength = worldBox->GetZHalfLength(); + } + else { + G4ExceptionDescription msg; + msg << "World volume of box shape not found." << G4endl; + msg << "Perhaps you have changed geometry." << G4endl; + msg << "The gun will be place in the center."; + G4Exception("PrimaryGeneratorAction::GeneratePrimaries()", + "MyCode0002", JustWarning, msg); + } + if(partSpecies.size()){ //only change impact point if there are multiple particles + //draw particle position from a unit distribution in x-y with +-5 cm around the center + G4double x = -1000; + G4double y = -1000; + G4double x_width = 10; + G4double y_width = 10; + while(x < -x_width/2. || x > x_width/2. || y < -y_width/2. || y > y_width/2.){ + x = G4INCL::Random::shoot()*x_width - x_width/2; + y = G4INCL::Random::shoot()*y_width - y_width/2; + } + fParticleGun->SetParticlePosition(G4ThreeVector(x*cm, y*cm, -worldZHalfLength)); + + } + else{ + fParticleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, -worldZHalfLength)); + } + //printout for checking + + fParticleGun->GeneratePrimaryVertex(anEvent); } - else{ - fParticleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, -worldZHalfLength)); - } - //printout for checking - //G4cout << "shooting at " << x << " " << y << " " << -worldZHalfLength << G4endl; - - fParticleGun->GeneratePrimaryVertex(anEvent); }