output of energy of muons x angle

include/B1Hits.hh

@@ -12,13 +12,24 @@ public:
 
   void print();
   void set_partdef (const G4String particle_name);
-
+  void set_energy (const G4double particle_en);
+  void set_position (const G4ThreeVector particle_pos);
   const G4String getParticleInTarget() {
     return fParticleInTarget;
   }
+  const G4double getParticleEnergy() {
+    return fparticle_energy;
+  }
+  const G4ThreeVector getParticlePos() {
+    return fparticle_position;;
+  }
+  
+  
 private:
 
  G4String fParticleInTarget;
+ G4double fparticle_energy;
+ G4ThreeVector fparticle_position;
 
 };
 

include/B1PrimaryGeneratorAction.hh

@@ -55,9 +55,11 @@ class B1PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     // method to access particle gun
     const G4ParticleGun* GetParticleGun() const { return fParticleGun; }
     
+  
   private:
     G4ParticleGun*  fParticleGun; // pointer a to G4 gun class
     G4Box* fEnvelopeBox;
+    double angle;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

include/B1RunAction.hh

@@ -54,10 +54,14 @@ class B1RunAction : public G4UserRunAction
 
     void AddEdep1 (G4double edep); 
     void AddEdep2 (G4double edep);
+    void AddE_mum (G4double edep); 
+    void AddE_mup (G4double edep);
 
   private:
     G4Accumulable<G4double> fEdep1;
     G4Accumulable<G4double> fEdep2;
+    G4Accumulable<G4double> fE_mum;
+    G4Accumulable<G4double> fE_mup;
 };
 
 #endif

src/B1EventAction.cc

@@ -31,7 +31,7 @@
 #include "B1RunAction.hh"
 #include "B1Hits.hh"
 #include "B1DetectorConstruction.hh"
-
+#include "G4SystemOfUnits.hh"
 #include "G4HCofThisEvent.hh"
 #include "G4VHitsCollection.hh"
 #include "G4SDManager.hh"
@@ -58,7 +58,7 @@ B1EventAction::~B1EventAction()
 void B1EventAction::BeginOfEventAction(const G4Event*)
 {    
   fEdep1 = 0.;
-  fEdep2 = 0;
+  fEdep2 = 0.;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -93,12 +93,42 @@ if (CHCID < 0) {
      G4cout << "My detector has " << n_hit << "hits" << G4endl;
      B1Hits* hit = new B1Hits;
      std::map<const G4String, int> fparticles;
+
+     double n_mu_p = 0.0;
+     double n_mu_m = 0.0;
+     G4double total_energy_mu_p = 0;
+     G4double total_energy_mu_m = 0;
+
+
      for(int i1 = 0; i1 < n_hit; i1++) {
       B1Hits* hit = (*HitsCol)[i1];
       const G4String name = hit->getParticleInTarget();
       G4cout << name << G4endl;
+      if (name == "mu+" || name=="mu-") {
+	 G4double energy = hit->getParticleEnergy();
+         G4ThreeVector position = hit->getParticlePos();
+	if (name == "mu+") {
+	n_mu_p++;
+ 	total_energy_mu_p += energy;
+	} else if (name=="mu-") {
+        total_energy_mu_m += energy;
+	n_mu_m++;
+
+	}
+    }
+      
 }
 
+ if(n_mu_p != 0.0) {
+       G4double value = total_energy_mu_p/n_mu_p;
+       G4cout << value/(GeV) << G4endl;
+       fRunAction->AddE_mup(value);
+    } else {fRunAction->AddE_mup(total_energy_mu_p);}
+ 
+ if(n_mu_m != 0.0) {
+       G4double value = total_energy_mu_m/n_mu_m;
+       fRunAction->AddE_mup(value);
+    } else {fRunAction->AddE_mum(total_energy_mu_m);}
 
 } 
 

src/B1Hits.cc

@@ -24,3 +24,12 @@ B1Hits::~B1Hits() {
       fParticleInTarget = particle_name;
 }
 
+ void B1Hits::set_energy(const G4double particle_energy) {
+      fparticle_energy = particle_energy;
+
+}
+
+void B1Hits::set_position(const G4ThreeVector particle_position) {
+      fparticle_position = particle_position;
+
+}

src/B1PrimaryGeneratorAction.cc

@@ -59,10 +59,10 @@ B1PrimaryGeneratorAction::B1PrimaryGeneratorAction()
   //ângulo de incidência
   #define pi 3.14159265
 
-  double angle = 0; // em graus
+  angle = 0; // em graus
   double p_y = tan(angle*pi/180.0);
   
-  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0,0.,1.));
+  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0,p_y,1.));
   fParticleGun->SetParticleEnergy(1.*GeV);
 }
 

src/B1RunAction.cc

@@ -31,7 +31,8 @@
 #include "B1PrimaryGeneratorAction.hh"
 #include "B1DetectorConstruction.hh"
 // #include "B1Run.hh"
-
+#include <stdio.h>      /* printf */
+#include <math.h> 
 #include "G4RunManager.hh"
 #include "G4Run.hh"
 #include "G4AccumulableManager.hh"
@@ -45,7 +46,9 @@
 B1RunAction::B1RunAction()
 : G4UserRunAction(),
   fEdep1(0.),
-  fEdep2(0.)
+  fEdep2(0.),
+  fE_mup(0.),
+  fE_mum(0.)
 { 
   // add new units for dose
   // 
@@ -63,6 +66,8 @@ B1RunAction::B1RunAction()
   G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance();
   accumulableManager->RegisterAccumulable(fEdep1);
   accumulableManager->RegisterAccumulable(fEdep2);
+  accumulableManager->RegisterAccumulable(fE_mup);
+  accumulableManager->RegisterAccumulable(fE_mum);
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -100,7 +105,8 @@ void B1RunAction::EndOfRunAction(const G4Run* run)
   //
   G4double edep  = fEdep1.GetValue() + fEdep2.GetValue();
   //G4double edep2 = fEdep2.GetValue();
-  
+  G4double mu_p_energy = fE_mup.GetValue();
+  G4double mu_m_energy = fE_mum.GetValue();
   //G4double rms = edep2 - edep*edep/nofEvents;
  // if (rms > 0.) rms = std::sqrt(rms); else rms = 0.;  
 
@@ -109,6 +115,8 @@ void B1RunAction::EndOfRunAction(const G4Run* run)
      (G4RunManager::GetRunManager()->GetUserDetectorConstruction());
   //G4double mass = detectorConstruction->GetScoringVolume()->GetMass();
   G4double dose = edep/nofEvents;
+  G4double mu_p_en = mu_p_energy/nofEvents;
+  G4double mu_m_en = mu_m_energy/nofEvents;
   //G4double rmsDose = rms/mass;
 
   // Run conditions
@@ -124,9 +132,29 @@ void B1RunAction::EndOfRunAction(const G4Run* run)
     runCondition += particleGun->GetParticleDefinition()->GetParticleName();
     runCondition += " of ";
     G4double particleEnergy = particleGun->GetParticleEnergy();
+    //G4ParticleMomentum direction = particleGun->GetParticleMomentumDirection();
+    G4double y = particleGun->GetParticleMomentumDirection().y();
+    G4double z = particleGun->GetParticleMomentumDirection().z();
+    
     runCondition += G4BestUnit(particleEnergy,"Energy");
+   
+   #define PI 3.14159265
+    G4double angle = atan(y/z)*180/PI;
+
+  std::ofstream mu_p_energy("energy_mu_p" + std::to_string(particleEnergy/GeV) + ".csv",std::ios_base::app);
+ 
+  std::ofstream mu_m_energy("energy_mu_m" + std::to_string(particleEnergy/GeV) + ".csv",std::ios_base::app);
+        
+   mu_p_energy << mu_p_en/GeV << " " << angle << "\n";
+   mu_m_energy << mu_m_en/GeV << " " << angle << "\n";
+
+  mu_p_energy.close();
+  mu_m_energy.close();
   }
  
+  //write on files
+ 
+  
   // Print
   //  
   if (IsMaster()) {
@@ -145,7 +173,7 @@ void B1RunAction::EndOfRunAction(const G4Run* run)
      << " The run consists of " << nofEvents << " "<< runCondition
      << G4endl
      << " Average Cumulated energy in volumes : " 
-     << G4BestUnit(dose,"Dose") << G4endl
+     << G4BestUnit(dose,"Energy") << G4endl
      << "------------------------------------------------------------"
      << G4endl
      << G4endl;
@@ -163,6 +191,16 @@ void B1RunAction::AddEdep2(G4double edep)
 {
   fEdep2  += edep;
   
+}
+void B1RunAction::AddE_mup(G4double edep)
+{
+  fE_mup  += edep;
+  
+}
+void B1RunAction::AddE_mum(G4double edep)
+{
+  fE_mum  += edep;
+  
 }
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 

src/B1SD.cc

@@ -26,11 +26,18 @@ G4bool B1SD::ProcessHits(G4Step* step, G4TouchableHistory* ROhist) {
   const G4String particle_name = step->GetTrack()->GetDynamicParticle()->GetParticleDefinition()->GetParticleName();
   //pega o Track relacionado à ela
   G4int track = step->GetTrack()->GetTrackID();
+//pega energia da particula
+ const G4double particle_energy = step->GetPreStepPoint()->GetTotalEnergy();
+//pega posição da particula
+const G4ThreeVector particle_position = step->GetPreStepPoint()->GetPosition();
+ 
 //Checa se a partícula é repetida
   if (track == track_id) {
     if (track == 0) {
       B1Hits* hit = new B1Hits();
       hit->set_partdef(particle_name);
+      hit->set_energy(particle_energy);
+      hit->set_position(particle_position);
       //coloca o valor na hitCollecion
       hitCollection->insert(hit);
       track_id = track;
@@ -40,6 +47,8 @@ G4bool B1SD::ProcessHits(G4Step* step, G4TouchableHistory* ROhist) {
   } else {
     B1Hits* hit = new B1Hits();
     hit->set_partdef(particle_name);
+    hit->set_energy(particle_energy);
+    hit->set_position(particle_position);
     hitCollection->insert(hit);
     track_id = track;
     return true;