position graph

2dplot.p

+set terminal png size 500,500
+set parametric
+set trange[0:2*pi]
+set size square
+set output 'hits_detector.png'
+
+z_0 = 0.9
+carbon_pDz = 1.125
+concrete_pDz = 0.76
+pos_after_detec = 3
+initial_angle = 2
+final_angle = 10
+
+initial_radius = (z_0 + 2*carbon_pDz + 2*concrete_pDz + pos_after_detec)*tan(initial_angle*pi/180.00)
+final_radius = (z_0 + 2*carbon_pDz + 2*concrete_pDz + pos_after_detec)*tan(final_angle*pi/180.00)
+
+plot initial_radius*cos(t), initial_radius*sin(t) title " ", final_radius*cos(t) , final_radius*sin(t) title " ", "position_mu_p.csv" u 1:2 title "mu+", "position_mu_m.csv" u 1:2 title "mu-"
+

CMakeLists.txt

@@ -49,6 +49,9 @@ set(EXAMPLEB1_SCRIPTS
   run1.mac
   run2.mac
   vis.mac
+  events.in
+  angle_var.in
+  2dplot.p
   )
 
 foreach(_script ${EXAMPLEB1_SCRIPTS})

angle_var.in

+# Macro file for absorver changind 1 degress starting from 2deg to 10deg
+
+/run/initialize
+
+/gun/energy 100 GeV
+/gun/particle mu+
+
+#2deg
+/gun/direction 0 0.036 1
+/run/beamOn 300
+
+#3deg
+/gun/direction 0 0.052 1
+/run/beamOn 300
+
+#4deg
+/gun/direction 0 0.069 1
+/run/beamOn 300
+
+#5deg
+/gun/direction 0 0.087 1
+/run/beamOn 300
+
+#6deg
+/gun/direction 0 0.105 1
+/run/beamOn 300
+
+#7deg
+/gun/direction 0 0.122 1
+/run/beamOn 300 
+
+#8deg
+/gun/direction 0 0.14 1
+/run/beamOn 300
+
+#9deg
+/gun/direction 0 0.158 1
+
+#10deg
+/gun/direction 0 0.176 1
+/run/beamOn 100
+
+
+
+/run/beamOn 100

events.in

+# Macro file for running a lot of events in a certain angle
+
+/run/initialize
+
+
+/gun/energy 100 GeV
+#5deg
+/gun/direction 0 0.087 1
+
+
+/gun/particle mu+
+/run/beamOn 300
+
+/gun/particle mu-
+/run/beamOn 300

include/B1EventAction.hh

@@ -29,7 +29,7 @@
 
 #ifndef B1EventAction_h
 #define B1EventAction_h 1
-
+#include "B1PrimaryGeneratorAction.hh"
 #include "G4UserEventAction.hh"
 #include "globals.hh"
 #include <iostream>
@@ -41,7 +41,7 @@ class B1RunAction;
 class B1DetectorConstruction;
 /// Event action class
 ///
-
+class G4Run;
 using namespace std;
 class B1EventAction : public G4UserEventAction
 {

src/B1EventAction.cc

@@ -80,7 +80,7 @@ if (CHCID < 0) {
  
  G4SDManager * SDman = G4SDManager::GetSDMpointer();
  G4HCofThisEvent* HCE = event->GetHCofThisEvent();
-
+ 
  B1HitsCollection* HitsCol = 0;
 
 
@@ -99,6 +99,19 @@ if (CHCID < 0) {
      G4double total_energy_mu_p = 0;
      G4double total_energy_mu_m = 0;
 
+    const B1PrimaryGeneratorAction* generatorAction
+   = static_cast<const B1PrimaryGeneratorAction*>
+     (G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
+ 
+
+   const G4ParticleGun* particleGun = generatorAction->GetParticleGun();
+
+
+   G4double particleEnergy = particleGun->GetParticleEnergy();
+     
+     std::ofstream mu_p_pos("position_mu_p.csv",std::ios_base::app);
+ 
+  std::ofstream mu_m_pos("position_mu_m.csv",std::ios_base::app);
 
      for(int i1 = 0; i1 < n_hit; i1++) {
       B1Hits* hit = (*HitsCol)[i1];
@@ -109,9 +122,11 @@ if (CHCID < 0) {
          G4ThreeVector position = hit->getParticlePos();
 	if (name == "mu+") {
 	n_mu_p++;
+        mu_p_pos << position.x()/(m) << "  " << position.y()/(m) << "\n";
  	total_energy_mu_p += energy;
 	} else if (name=="mu-") {
         total_energy_mu_m += energy;
+        mu_m_pos << position.x()/(m) << "  " << position.y()/(m) << "\n";
 	n_mu_m++;
 
 	}
@@ -119,6 +134,10 @@ if (CHCID < 0) {
       
 }
 
+ mu_p_pos.close();
+mu_m_pos.close();
+
+
  if(n_mu_p != 0.0) {
        G4double value = total_energy_mu_p/n_mu_p;
        G4cout << value/(GeV) << G4endl;