From 85673e78b471c6bcebe4864fc0019cdb8ef22197 Mon Sep 17 00:00:00 2001
From: Elidiano Tronchin <elidiano.tronchin@gmail.com>
Date: Fri, 12 Oct 2018 10:11:59 +0200
Subject: [PATCH] * Chqnge to generate muon

---
 NPSimulation/Core/ParticleStack.cc               |  2 ++
 NPSimulation/Detectors/Dali/Dali.cc              |  1 -
 .../EventGenerator/EventGeneratorCosmic.cc       | 16 +++++++++-------
 .../EventGenerator/EventGeneratorIsotropic.cc    |  4 +++-
 .../EventGeneratorMultipleParticle.cc            |  3 ++-
 NPSimulation/Process/PhysicsList.cc              |  6 ++++++
 NPSimulation/Simulation.cc                       |  9 +++++++++
 7 files changed, 31 insertions(+), 10 deletions(-)

diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index c8d6641f9..769fcb2bc 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -240,6 +240,8 @@ string ParticleStack::ChangeNameToG4Standard(string OriginalName){
     else if (FinalName=="deuteron") FinalName="deuteron";
     else if (FinalName=="triton")   FinalName="triton";
     else if (FinalName=="alpha")    FinalName="alpha";
+    else if (FinalName=="mu+")   FinalName="mu+";
+    else if (FinalName=="mu-")   FinalName="mu-";
     else if (FinalName=="n")        FinalName="neutron";
     else if (FinalName=="neutron")  FinalName="neutron";
     return(FinalName);
diff --git a/NPSimulation/Detectors/Dali/Dali.cc b/NPSimulation/Detectors/Dali/Dali.cc
index af3aa10ac..fc7e58ba6 100644
--- a/NPSimulation/Detectors/Dali/Dali.cc
+++ b/NPSimulation/Detectors/Dali/Dali.cc
@@ -368,7 +368,6 @@ void Dali::ReadSensitive(const G4Event* ){
 
   
   unsigned int size = Scorer->GetMult();
-  // cout << "size " << size << endl;
   for(unsigned int i = 0 ; i < size ; i++){
     vector<unsigned int> level = Scorer->GetLevel(i); 
     double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Dali_NS::ResoEnergy);
diff --git a/NPSimulation/EventGenerator/EventGeneratorCosmic.cc b/NPSimulation/EventGenerator/EventGeneratorCosmic.cc
index e5e68318c..c5ce9d1a0 100644
--- a/NPSimulation/EventGenerator/EventGeneratorCosmic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorCosmic.cc
@@ -174,8 +174,8 @@ void EventGeneratorCosmic::GenerateEvent(G4Event*){
 
         /* //Putting a cylinder
         if(randomize>0){          //top
-        momentum_x = cos(angle)*CosmicAngle;
-        momentum_z = sin(angle)*CosmicAngle;
+        momentum_x = cos(angle)*sin(CosmicAngle);
+        momentum_z = sin(angle)*sin(CosmicAngle);
           
         x0 = cos(angle2*2)*R*(randomize*2);
         z0 = sin(angle2*2)*R*(randomize*2);
@@ -186,17 +186,19 @@ void EventGeneratorCosmic::GenerateEvent(G4Event*){
           z0 = sin(angle)*R;
           par.m_y0 = (.5-RandFlat::shoot())*H;   ///!!!!!!!
 
-          momentum_x = -cos(angle+angle2)*CosmicAngle;
-          momentum_z = -sin(angle+angle2)*CosmicAngle;
+          momentum_x = -cos(angle+angle2)*sin(CosmicAngle);
+          momentum_z = -sin(angle+angle2)*sin(CosmicAngle);
         }
         */
 
         // Begin Constrain to pass in a circle with radius 3R
-        momentum_x = cos(angle)*CosmicAngle;
-        momentum_z = sin(angle)*CosmicAngle;
+
+        momentum_y = cos(CosmicAngle);
+        momentum_x = cos(angle)*sin(CosmicAngle);
+        momentum_z = sin(angle)*sin(CosmicAngle);
         x0 = cos(angle2*2)*R*(0.5-randomize)*3-momentum_x*( H/2 / momentum_y);// *( H/2 / momentum_y) this is to have the origin always with par.m_y0 = H/2;
         z0 = sin(angle2*2)*R*(0.5-randomize)*3-momentum_z*( H/2 / momentum_y);
-        par.m_y0 = H/2;
+        par.m_y0 = H/2; // momentum_y*( H/2 / momentum_y);
         // End Constrain to pass in a circle with radius 3R
 
 
diff --git a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
index db497821f..1647dc6cb 100644
--- a/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorIsotropic.cc
@@ -91,6 +91,8 @@ void EventGeneratorIsotropic::ReadConfiguration(NPL::InputParser parser){
         else if(particleName[j]=="3He" || particleName[j]=="He3") { it->m_particleName.push_back("3He") ; }
         else if(particleName[j]=="alpha") { it->m_particleName.push_back("4He") ; }
         else if(particleName[j]=="gamma") { it->m_particleName.push_back("gamma") ;}
+        else if(particleName[j]=="mu+") { it->m_particleName.push_back("mu+") ;}
+        else if(particleName[j]=="mu-") { it->m_particleName.push_back("mu-") ;}
         else if(particleName[j]=="neutron") {it->m_particleName.push_back("neutron") ;}
         else it->m_particleName.push_back(particleName[j]);
       }
@@ -128,7 +130,7 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event*){
         par.m_particle=NULL;
         if(par.m_particle==NULL){
 
-          if(par.m_particleName[p]=="gamma" || par.m_particleName[p]=="neutron" ||  par.m_particleName[p]=="opticalphoton"  ||  par.m_particleName[p]=="mu+"){
+          if(par.m_particleName[p]=="gamma" || par.m_particleName[p]=="neutron" ||  par.m_particleName[p]=="opticalphoton"  ||  par.m_particleName[p]=="mu+" ||  par.m_particleName[p]=="mu-"){
             par.m_particle =  G4ParticleTable::GetParticleTable()->FindParticle(par.m_particleName[p].c_str());
           }
           else{
diff --git a/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc b/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc
index f9f08578b..a4d56722a 100644
--- a/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorMultipleParticle.cc
@@ -118,6 +118,7 @@ void EventGeneratorMultipleParticle::ReadConfiguration(NPL::InputParser parser){
                     else if(sParticle=="3He" || sParticle=="He3"){vParticle.push_back("3He");}
                     else if(sParticle=="alpha"){vParticle.push_back("4He");}
                     else if(sParticle=="gamma") { vParticle.push_back("gamma") ;}
+                    else if(sParticle=="mu+") { vParticle.push_back("mu+") ;}
                     else if(sParticle=="neutron") {vParticle.push_back("neutron") ;}
                     else vParticle.push_back(sParticle);
                     
@@ -154,7 +155,7 @@ void EventGeneratorMultipleParticle::GenerateEvent(G4Event* evt){
     for(int i=0; i<m_Multiplicity[evtID]; i++){
         m_particle=NULL;
         if(m_particle==NULL){
-            if(m_particleName[evtID][i]=="gamma" || m_particleName[evtID][i]=="neutron" ||  m_particleName[evtID][i]=="opticalphoton"){
+            if(m_particleName[evtID][i]=="gamma" || m_particleName[evtID][i]=="neutron" ||  m_particleName[evtID][i]=="opticalphoton" ||  m_particleName[evtID][i]=="mu+"){
                 m_particle =  G4ParticleTable::GetParticleTable()->FindParticle(m_particleName[evtID][i].c_str());
             }
             else{
diff --git a/NPSimulation/Process/PhysicsList.cc b/NPSimulation/Process/PhysicsList.cc
index 3ab889c95..ee7b76530 100644
--- a/NPSimulation/Process/PhysicsList.cc
+++ b/NPSimulation/Process/PhysicsList.cc
@@ -43,6 +43,7 @@
 #include "G4IonFluctuations.hh"
 #include "G4IonParametrisedLossModel.hh"
 #include "G4UniversalFluctuation.hh"
+#include "G4UImanager.hh"
 
 #include "G4PAIModel.hh"
 #include "G4PAIPhotModel.hh"
@@ -50,6 +51,11 @@
 
 /////////////////////////////////////////////////////////////////////////////
 PhysicsList::PhysicsList() : G4VUserPhysicsList(){
+
+  // UI commands to activate one or more High Energy Processes
+  //  G4UImanager* UI = G4UImanager::GetUIpointer();
+  //  UI->ApplyCommand("/physics_list/em/MuonNuclear true");
+    
     m_EmList = "Option4";
     defaultCutValue = 1*mm;//0.2*mm;
     opticalPhysicsList = NULL;
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index 03a09d5c5..176481979 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -6,6 +6,8 @@
 #include "G4PhysListFactory.hh"
 // UI
 #include "G4UImanager.hh"
+#include "QBBC.hh"
+
 #include "G4UIterminal.hh"
 #include "G4UItcsh.hh"
 
@@ -72,7 +74,14 @@ int main(int argc, char** argv){
     runManager->SetUserInitialization(detector);
     
     PhysicsList* physicsList   = new PhysicsList();
+    //G4VModularPhysicsList* physicsList   = new QBBC;
+    physicsList->SetVerboseLevel(0);
     runManager->SetUserInitialization(physicsList);
+        
+
+    
+
+    
     PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector);
 
     // Initialize Geant4 kernel
-- 
GitLab