added multi particles

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
+    if (not isinstance(particleSpec, list)) or len(particleSpec) == 1:
-    direction_vector = [(end_point[0] - start_point[0]), (end_point[1] - start_point[1]), (end_point[2] - start_point[2])]
+        # Calculate the direction vector for the arrow
-    # Create the arrowhead at the start point with the opposite direction
+        direction_vector = [(end_point[0] - start_point[0]), (end_point[1] - start_point[1]), (end_point[2] - start_point[2])]
-    arrowhead_trace = go.Cone(
+        # Create the arrowhead at the start point with the opposite direction
-        x=[end_point[0]],
+        arrowhead_trace = go.Cone(
-        z=[end_point[1]],
+            x=[end_point[0]],
-        y=[end_point[2]],
+            z=[end_point[1]],
-        u=[direction_vector[0]],
+            y=[end_point[2]],
-        w=[direction_vector[1]],
+            u=[direction_vector[0]],
-        v=[direction_vector[2]],
+            w=[direction_vector[1]],
-        sizemode='scaled',
+            v=[direction_vector[2]],
-        sizeref=0.8,
+            sizemode='scaled',
-        showscale=False,
+            sizeref=0.8,
-        colorscale='Reds',
+            showscale=False,
-        opacity=1.0,
+            colorscale='Reds',
-        anchor='tail',
+            opacity=1.0,
-    )
+            anchor='tail',
-    # Create the 3D scatter plot with both traces
+        )
-    to_plot.append(line_trace)
+        # Create the 3D scatter plot with both traces
-    to_plot.append(arrowhead_trace)
+        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")
+    if len(outfile) > 0:
-    #exit()
+        fig.write_html(outfile)
-    if renderer is not None:
+    
-        fig.show(renderer=renderer) 
     else:
-        fig.show() 
+        if renderer is not None:
+            fig.show(renderer=renderer) 
+        else:
+            fig.show() 
 
 
 

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<std::string> parts = {"e-"};
+  builder.run_batch(10000,parts , 1, 100);
 
   //run(cw, 1000, "e-", 1, 100, true);
 

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<std::string>& 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<std::string>{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<std::string> partSpecies;
     G4String filename="";
 };
 

include/G4System.hh

@@ -17,6 +17,7 @@
 
 #include "FTFP_BERT.hh"
 #include "Randomize.hh"
+#include <vector>
 
 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<std::string>& partSpecies, double minEnergy_GeV, 
            double maxEnergy_GeV, std::string filename="_1234567890_Hits.root");
 
 void applyUICommand(const std::string& command){

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){
+  void setPartSpecies(const std::vector<std::string>& values){
-    partSpecies = value;
+    partSpecies = values;
-    if(partSpecies == "e-"){
+    for(const auto value : values){
-      G4cout << "Particle species set to electron" << G4endl;
+      G4String val = value;
-    }
+      if(val == "e-"){
-    else if(partSpecies == "e+"){
+        G4cout << "Particle species set to electron" << G4endl;
-      G4cout << "Particle species set to positron" << G4endl;
+      }
-    }
+      else if(val == "e+"){
-    else if(partSpecies == "gamma"){
+        G4cout << "Particle species set to positron" << G4endl;
-      G4cout << "Particle species set to gamma" << G4endl;
+      }
-    }
+      else if(val == "gamma"){
-    else if(partSpecies == "proton"){
+        G4cout << "Particle species set to gamma" << G4endl;
-      G4cout << "Particle species set to proton" << G4endl;
+      }
-    }
+      else if(val == "proton"){
-    else if(partSpecies == "neutron"){
+        G4cout << "Particle species set to proton" << G4endl;
-      G4cout << "Particle species set to neutron" << G4endl;
+      }
-    }
+      else if(val == "neutron"){
-    else if(partSpecies == "pi+"){
+        G4cout << "Particle species set to neutron" << G4endl;
-      G4cout << "Particle species set to charged pion" << G4endl;
+      }
-    }
+      else if(val == "pi+"){
-    else if(partSpecies == "pi-"){
+        G4cout << "Particle species set to charged pion" << G4endl;
-      G4cout << "Particle species set to charged pion" << G4endl;
+      }
-    }
+      else if(val == "pi-"){
-    else{
+        G4cout << "Particle species set to charged pion" << G4endl;
-      G4cout << "Particle species "<< partSpecies << " not recognized" << G4endl;
+      }
-      throw std::invalid_argument("Particle species not recognized");
+      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<std::string> partSpecies;
 
 };
 

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<std::string>& 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);
 }

src/PrimaryGeneratorAction.cc

@@ -55,7 +55,7 @@ PrimaryGeneratorAction::PrimaryGeneratorAction()
   G4int nofParticles = 1;
   fParticleGun = new G4ParticleGun(nofParticles);
 
-  partSpecies = "e-";
+  partSpecies = std::vector<std::string>{"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;
+  //create fractions for each entry in the partSpecies vector such that the partEnergy can be distributed 
-
+  //according to the fractions
-  auto particleDefinition
+  std::vector<G4double> part_energies;
-    = G4ParticleTable::GetParticleTable()->FindParticle(partSpecies);
+  G4double sum = 0;
-  if(particleDefinition == nullptr){
+  for(const auto& part : partSpecies){
-    G4cout << "Particle " << partSpecies << " not found" << G4endl;
+    G4double randfrac = 0;
-    throw std::runtime_error("Particle not found");
+    while(randfrac > 1 || randfrac <= 0)
-    return;
+        randfrac = G4INCL::Random::shoot();
+    part_energies.push_back(randfrac);
+    sum += randfrac;
   }
-  fParticleGun->SetParticleDefinition(particleDefinition);
+  //normalise w.r.t total 
-  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
+  for(auto& en : part_energies){
-  fParticleGun->SetParticleEnergy(partEnergy*GeV);
+    en /= sum;
-
+    en *= partEnergy;
-  // 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<G4Box*>(worldLV->GetSolid());
   }
 
-  if ( worldBox ) {
+  // now create a vertex for every particle
-    worldZHalfLength = worldBox->GetZHalfLength();
+
-  }
+  for (size_t i = 0; i < partSpecies.size(); i++){
-  else  {
+    G4String single_partType = partSpecies[i];
-    G4ExceptionDescription msg;
+    G4double single_partEnergy = part_energies[i];
-    msg << "World volume of box shape not found." << G4endl;
+    //create a vertex for the particle
-    msg << "Perhaps you have changed geometry." << G4endl;
+    auto particleDefinition
-    msg << "The gun will be place in the center.";
+      = G4ParticleTable::GetParticleTable()->FindParticle(single_partType);
-    G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
+    if(particleDefinition == nullptr){
-      "MyCode0002", JustWarning, msg);
+      G4cout << "Particle " << single_partType << " not found" << G4endl;
-  }
+      throw std::runtime_error("Particle not found");
-  if(false){
+      return;
-    //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;
     }
-    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<G4Box*>(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);
 
 }