diff --git a/NPLib/Makefile b/NPLib/Makefile
index 7412348ca5b938e88ba33453b67a719ef223a75f..41cede0101f4b850935eae500475eca163cd2fad 100644
--- a/NPLib/Makefile
+++ b/NPLib/Makefile
@@ -1,5 +1,5 @@
include $(NPTOOL)/NPLib/Makefile.arch
-
+#include /usr/local/root-5.32.01/etc/Makefile.arch
#------------------------------------------------------------------------------
SHARELIB = SharedLib
FILLINC = FillIncludeDir
diff --git a/NPLib/Makefile.arch b/NPLib/Makefile.arch
index 2d8aec4e0bb51415ed0d2c1dd3932814d18548d0..543ca7973398e093c958af35f59c71654fb20786 100644
--- a/NPLib/Makefile.arch
+++ b/NPLib/Makefile.arch
@@ -1,6 +1,6 @@
-# -*- mode: makefile -*-
-#
-# Makefile containing platform dependencies for ROOT based projects.
+# Makefile for the ROOT test programs.
+# This Makefile shows nicely how to compile and link applications
+# using the ROOT libraries on all supported platforms.
#
# Copyright (c) 2000 Rene Brun and Fons Rademakers
#
@@ -43,44 +43,8 @@ HASTHREAD := $(shell $(ROOTCONFIG) --has-thread)
ROOTDICTTYPE := $(shell $(ROOTCONFIG) --dicttype)
ROOTCINT := rootcint
-ifeq ($(PLATFORM),macosx)
-MACOSX_MINOR := $(shell sw_vers | sed -n 's/ProductVersion://p' | cut -d . -f 2)
-ifeq ($(MACOSX_DEPLOYMENT_TARGET),)
-MACOSXTARGET := 10.$(MACOSX_MINOR)
-else
-MACOSXTARGET := $(MACOSX_DEPLOYMENT_TARGET)
-endif
-endif
-
-ifeq ($(ARCH),hpuxacc)
-# HP-UX 10.x with aCC
-CXX = aCC
-CXXFLAGS = $(OPT) +Z
-LD = aCC
-LDFLAGS = $(OPT) -z
-SOFLAGS = -b
-endif
-
-ifeq ($(ARCH),hpuxia64acc)
-# HP-UX 11i 1.5 (IA-64) with aCC
-CXX = aCC
-CXXFLAGS = +DD64 $(OPT) +Z
-LD = aCC
-LDFLAGS = +DD64 $(OPT) -z
-SOFLAGS = -b
-endif
-
-ifeq ($(ARCH),hpuxgcc)
-# HP-UX 10.x with g++
-CXXFLAGS = $(OPT) -fPIC
-CXX = g++
-LD = g++
-LDFLAGS = $(OPT)
-SOFLAGS = -fPIC -shared
-endif
-
-ifeq ($(ARCH),hurddeb)
-# GNU/Hurd
+ifeq ($(ARCH),linux)
+# Linux with egcs, gcc 2.9x, gcc 3.x
CXX = g++
CXXFLAGS = $(OPT2) -Wall -fPIC
LD = g++
@@ -88,134 +52,13 @@ LDFLAGS = $(OPT2)
SOFLAGS = -shared
endif
-ifeq ($(ARCH),aix)
-# IBM AIX xlC 4.x
-CXX = xlC
-CXXFLAGS = $(OPT)
-LD = xlC
-LDFLAGS = $(OPT)
-SOFLAGS =
-DllSuf = a
-endif
-
-ifeq ($(ARCH),aix5)
-# IBM AIX xlC 5.x
-CXX = xlC
-CXXFLAGS = $(OPT)
-LD = xlC
-LDFLAGS = $(OPT)
+ifeq ($(ARCH),linuxkcc)
+# Linux with the KAI compiler
+CXX = KCC --one_instantiation_per_object
+CXXFLAGS = $(OPT) -fPIC +K0
+LD = KCC
+LDFLAGS = $(OPT) $(shell $(ROOTCONFIG) --cflags)
SOFLAGS =
-DllSuf = a
-endif
-
-ifeq ($(ARCH),aixgcc)
-# IBM AIX with GCC
-CXX = g++
-CXXFLAGS = $(OPT)
-LD = g++
-LDFLAGS = $(OPT)
-SOFLAGS = -shared
-DllSuf = a
-EXPLLINKLIBS = $(ROOTLIBS) $(ROOTGLIBS)
-endif
-
-ifeq ($(ARCH),solaris)
-# Solaris CC
-CXX = /opt/SUNWspro/bin/CC
-CXXFLAGS = $(OPT) -KPIC
-LD = /opt/SUNWspro/bin/CC
-LDFLAGS = $(OPT)
-SOFLAGS = -G
-endif
-
-ifeq ($(ARCH),solarisCC5)
-# Solaris CC 5.0
-CXX = CC
-CXXFLAGS = $(OPT) -KPIC
-LD = CC
-LDFLAGS = $(OPT)
-SOFLAGS = -G
-endif
-
-ifeq ($(ARCH),solaris64CC5)
-# Solaris CC 5.0 64-bit
-CXX = CC
-CXXFLAGS = $(OPT) -KPIC
-LD = CC
-LDFLAGS = $(OPT)
-SOFLAGS = -G
-endif
-
-ifeq ($(ARCH),solarisgcc)
-# Solaris gcc
-CXX = g++
-CXXFLAGS = $(OPT) -fPIC
-LD = g++
-LDFLAGS = $(OPT)
-SOFLAGS = -shared
-endif
-
-ifeq ($(ARCH),sgicc)
-# SGI
-CXX = CC -n32 -I/usr/include/CC.sgi
-CXXFLAGS = $(OPT)
-LD = CC -n32 -LANG:std -I/usr/include/CC.sgi
-LDFLAGS = $(OPT)
-SOFLAGS = -shared
-endif
-
-ifeq ($(ARCH),sgicc64)
-# SGI
-CXX = CC -64 -I/usr/include/CC.sgi
-CXXFLAGS = $(OPT)
-LD = CC -64 -LANG:std -I/usr/include/CC.sgi
-LDFLAGS = $(OPT)
-SOFLAGS = -shared
-endif
-
-ifeq ($(ARCH),sgigcc)
-# SGI 6.x with gcc
-CXX = g++
-CXXFLAGS = $(OPT) -Wall -fPIC
-LD = g++
-LDFLAGS = $(OPT) -Wl,-u,__builtin_new -Wl,-u,__builtin_delete -Wl,-u,__nw__FUiPv
-SOFLAGS = -shared
-endif
-
-ifeq ($(ARCH),sgin32gcc)
-# SGI 6.x with gcc for n32 ABI
-CXX = g++
-CXXFLAGS = $(OPT) -Wall -fPIC
-LD = g++
-LDFLAGS = $(OPT) -L/usr/lib32 -Wl,-woff,134 -lgen
-SOFLAGS = -shared
-endif
-
-ifeq ($(ARCH),alphagcc)
-# Alpha/OSF with gcc
-CXX = g++
-CXXFLAGS = $(OPT2) -Wall -fPIC
-LD = g++
-LDFLAGS = $(OPT2)
-SOFLAGS = -Wl,-expect_unresolved,* -shared
-endif
-
-ifeq ($(ARCH),alphacxx6)
-# Alpha/OSF with cxx6
-CXX = cxx
-CXXFLAGS = $(OPT)
-LD = cxx
-LDFLAGS = $(OPT)
-SOFLAGS = -shared -nocxxstd -Wl,-expect_unresolved,*,-msym
-endif
-
-ifeq ($(ARCH),linux)
-# Linux with egcs, gcc 2.9x, gcc 3.x
-CXX = g++
-CXXFLAGS = $(OPT2) -Wall -fPIC
-LD = g++
-LDFLAGS = $(OPT2)
-SOFLAGS = -shared
endif
ifeq ($(ARCH),linuxicc)
@@ -249,6 +92,15 @@ LDFLAGS = $(OPT2)
SOFLAGS = -shared
endif
+ifeq ($(ARCH),linuxia64sgi)
+# Itanium Linux with sgiCC
+CXX = sgiCC
+CXXFLAGS = $(OPT) -Wall -fPIC
+LD = gsgiCC
+LDFLAGS = $(OPT)
+SOFLAGS = -shared
+endif
+
ifeq ($(ARCH),linuxia64ecc)
# Itanium Linux with Intel icc (was ecc)
ICC_MAJOR := $(shell icc -v 2>&1 | awk '{ if (NR==1) print $$2 }' | \
@@ -325,48 +177,14 @@ LDFLAGS = $(OPT)
SOFLAGS = -shared
endif
-ifeq ($(ARCH),freebsd4)
-# FreeBSD with glibc
-CXX = g++
-CXXFLAGS = $(OPT) -W -Wall -fPIC
-LD = $(CXX)
-LDFLAGS = $(OPT)
-SOFLAGS = -shared -Wl,-x
-endif
-
-ifeq ($(ARCH),freebsd5)
-# FreeBSD with glibc
-CXX = g++
-CXXFLAGS = $(OPT) -W -Wall -fPIC
-LD = $(CXX)
-LDFLAGS = $(OPT)
-SOFLAGS = -shared -Wl,-x
-endif
-
-ifeq ($(ARCH),freebsd7)
-# FreeBSD with libc
-CXX = g++
-CXXFLAGS = $(OPT) -W -Wall -fPIC
-LD = $(CXX)
-LDFLAGS = $(OPT)
-SOFLAGS = -shared -Wl,-x
-endif
-
-ifeq ($(ARCH),openbsd)
-# OpenBSD with libc
-CXX = g++
-CXXFLAGS = $(OPT) -pipe -W -Wall -fPIC
-LD = g++
-LDFLAGS = $(OPT)
-SOFLAGS = -shared -Wl,-x
-endif
-
ifeq ($(ARCH),macosx)
# MacOS X with cc (GNU cc 2.95.2 and gcc 3.3)
+MACOSX_MINOR := $(shell sw_vers | sed -n 's/ProductVersion://p' | cut -d . -f 2)
+MACOSXTARGET := MACOSX_DEPLOYMENT_TARGET=10.$(MACOSX_MINOR)
CXX = g++
CXXFLAGS = $(OPT2) -pipe -Wall -W -Woverloaded-virtual
-LD = g++
-LDFLAGS = $(OPT2) -mmacosx-version-min=$(MACOSXTARGET) -bind_at_load
+LD = $(MACOSXTARGET) g++
+LDFLAGS = $(OPT2) -bind_at_load
# The SOFLAGS will be used to create the .dylib,
# the .so will be created separately
ifeq ($(subst $(MACOSX_MINOR),,1234),1234)
@@ -379,19 +197,20 @@ ifneq ($(subst $(MACOSX_MINOR),,12),12)
UNDEFOPT = suppress
LD = g++
endif
-#SOFLAGS = -dynamiclib -single_module -undefined $(UNDEFOPT) -install_name $(CURDIR)/
SOFLAGS = -dynamiclib -single_module -undefined $(UNDEFOPT)
endif
ifeq ($(ARCH),macosxicc)
# MacOS X with Intel icc compiler
+MACOSX_MINOR := $(shell sw_vers | sed -n 's/ProductVersion://p' | cut -d . -f 2)
+MACOSXTARGET := MACOSX_DEPLOYMENT_TARGET=10.$(MACOSX_MINOR)
ifeq ($(MACOSX_MINOR),5)
MACOSX_MINOR := 4
endif
CXX = icc
CXXFLAGS = $(OPT) -fPIC -wd1476
-LD = icpc
-LDFLAGS = $(OPT2) -mmacosx-version-min=$(MACOSXTARGET)
+LD = $(MACOSXTARGET) icpc
+LDFLAGS = $(OPT)
# The SOFLAGS will be used to create the .dylib,
# the .so will be created separately
ifeq ($(subst $(MACOSX_MINOR),,1234),1234)
@@ -399,17 +218,18 @@ DllSuf = so
else
DllSuf = dylib
endif
-#SOFLAGS = -dynamiclib -single_module -undefined dynamic_lookup -install_name $(CURDIR)/
-SOFLAGS = -dynamiclib -single_module -undefined $(UNDEFOPT)
+SOFLAGS = -dynamiclib -single_module -undefined dynamic_lookup
endif
ifeq ($(ARCH),macosx64)
# MacOS X >= 10.4 with gcc 64 bit mode (GNU gcc 4.*)
# Only specific option (-m64) comes from root-config
+MACOSX_MINOR := $(shell sw_vers | sed -n 's/ProductVersion://p' | cut -d . -f 2)
+MACOSXTARGET := MACOSX_DEPLOYMENT_TARGET=10.$(MACOSX_MINOR)
CXX = g++
CXXFLAGS = $(OPT2) -pipe -Wall -W -Woverloaded-virtual
-LD = g++
-LDFLAGS = $(OPT2) -mmacosx-version-min=$(MACOSXTARGET) -bind_at_load
+LD = $(MACOSXTARGET) g++ -m64
+LDFLAGS = $(OPT2) -bind_at_load -e start
# The SOFLAGS will be used to create the .dylib,
# the .so will be created separately
ifeq ($(subst $(MACOSX_MINOR),,1234),1234)
@@ -417,65 +237,26 @@ DllSuf = so
else
DllSuf = dylib
endif
-#SOFLAGS = -dynamiclib -single_module -undefined dynamic_lookup -install_name $(CURDIR)/
-SOFLAGS = -dynamiclib -single_module -undefined $(UNDEFOPT)
-endif
-
-ifeq ($(ARCH),win32)
-# Windows with the VC++ compiler
-VC_MAJOR := $(shell unset VS_UNICODE_OUTPUT; cl.exe 2>&1 | awk '{ if (NR==1) print $$(NF-2) }' | \
- cut -d'.' -f1)
-ObjSuf = obj
-SrcSuf = cxx
-ExeSuf = .exe
-DllSuf = dll
-OutPutOpt = -out:
-CXX = cl
-ifeq (debug,$(findstring debug,$(ROOTBUILD)))
-CXXOPT = -Z7
-LDOPT = -debug
-else
-ifneq ($(findstring debug, $(strip $(shell $(ROOTCONFIG) --config))),)
-CXXOPT = -Z7
-LDOPT = -debug
-else
-CXXOPT = -O2
-LDOPT = -opt:ref
-endif
-endif
-ROOTINCDIR := -I$(shell cygpath -m `$(ROOTCONFIG) --incdir`)
-CXXFLAGS = $(CXXOPT) -nologo $(ROOTINCDIR) -FIw32pragma.h
-LD = link
-LDFLAGS = $(LDOPT) -nologo
-SOFLAGS = -DLL
-
-EXPLLINKLIBS = $(ROOTLIBS) $(ROOTGLIBS)
-ifneq (,$(findstring $(VC_MAJOR),14 15 16))
-MT_EXE = mt -nologo -manifest $@.manifest -outputresource:$@\;1; rm -f $@.manifest
-MT_DLL = mt -nologo -manifest $@.manifest -outputresource:$@\;2; rm -f $@.manifest
-else
-MT_EXE =
-MT_DLL =
-endif
+SOFLAGS = -m64 -dynamiclib -single_module -undefined dynamic_lookup
endif
-ifeq ($(ARCH),win32gcc)
-# Windows with gcc
-DllSuf = dll
-ExeSuf = .exe
-CXX = g++
-CXXFLAGS = $(OPT) -pipe -Wall -Woverloaded-virtual -I/usr/X11R6/include
-LD = g++
-LDFLAGS = $(OPT) -Wl,--enable-auto-import \
- -Wl,--enable-runtime-pseudo-reloc \
- -L/usr/X11R6/lib
-SOFLAGS = -shared -Wl,--enable-auto-image-base \
- -Wl,--export-all-symbols
-EXPLLINKLIBS = $(ROOTLIBS) $(ROOTGLIBS)
+ifeq ($(ARCH),macosxxlc)
+# MacOS X with IBM xlC compiler
+MACOSX_MINOR := $(shell sw_vers | sed -n 's/ProductVersion://p' | cut -d . -f 2)
+MACOSXTARGET := MACOSX_DEPLOYMENT_TARGET=10.$(MACOSX_MINOR)
+CXX = xlC
+CXXFLAGS = $(OPT)
+LD = $(MACOSXTARGET) xlC
+LDFLAGS = $(OPT) -Wl,-bind_at_load
+# The SOFLAGS will be used to create the .dylib,
+# the .so will be created separately
+DllSuf = dylib
+UNDEFOPT = dynamic_lookup
+ifneq ($(subst $(MACOSX_MINOR),,12),12)
+UNDEFOPT = suppress
+LD = xlC
endif
-
-ifeq ($(CXX),)
-$(error $(ARCH) invalid architecture)
+SOFLAGS = -qmkshrobj -single_module -undefined $(UNDEFOPT)
endif
CXXFLAGS += $(ROOTCFLAGS)
@@ -485,20 +266,3 @@ GLIBS = $(ROOTGLIBS) $(SYSLIBS)
INCLUDE = -I$(CLHEP_BASE_DIR)/include -I$(NPTOOL)/NPLib/include
-ifneq ($(ALTCC),)
- CC = $(ALTCC)
-endif
-ifneq ($(ALTCXX),)
- CXX = $(ALTCXX)
-endif
-ifneq ($(ALTF77),)
- F77 = $(ALTF77)
-endif
-ifneq ($(ALTLD),)
- LD = $(ALTLD)
-endif
-
-ifneq ($(findstring g++, $(CXX)),)
-GCC_MAJOR := $(shell $(CXX) -dumpversion 2>&1 | cut -d'.' -f1)
-GCC_MINOR := $(shell $(CXX) -dumpversion 2>&1 | cut -d'.' -f2)
-endif
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 781822ee74debc7753028f1e8a4fd508735d42e1..4937e14f509c17f99916921a0b3171ed554e6564 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -80,7 +80,9 @@ namespace NPL
void SetExcitationHeavy(double exci) {fExcitationHeavy = exci; initializePrecomputeVariable();}
double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;}
- double GetExcitation() const {return fExcitation;}
+ double GetExcitationHeavy() const {return fExcitationHeavy;}
+ double GetExcitation() const {return fExcitationHeavy;}
+ double GetExcitationLight() const {return fExcitationLight;}
double GetQValue() const {return fQValue;}
Nucleus* GetNucleus1() const {return fNuclei1;}
Nucleus* GetNucleus2() const {return fNuclei2;}
diff --git a/NPLib/scripts/buildliblist.sh b/NPLib/scripts/buildliblist.sh
index b12480117a53c7d0c3272e10fb8a331f4d7f26f3..b19d77cf21c3578ade02e8747219a87a07583000 100755
--- a/NPLib/scripts/buildliblist.sh
+++ b/NPLib/scripts/buildliblist.sh
@@ -50,7 +50,7 @@ do
# add trailing \
name="$name \\"
# add tab at the beginning
- name=$(echo -e "\t $name")
- echo -e "\t $name" >> $outfile
+ name=$(echo "\t $name")
+ echo "\t $name" >> $outfile
fi ;
done
diff --git a/NPSimulation/include/EventGeneratorTransfert.hh b/NPSimulation/include/EventGeneratorTransfert.hh
index c9db4789b4dcacb3baa093095d00f7378979f5e7..813a57dbac778ed4943db99dc848a3588b1e211b 100644
--- a/NPSimulation/include/EventGeneratorTransfert.hh
+++ b/NPSimulation/include/EventGeneratorTransfert.hh
@@ -32,9 +32,6 @@
#include "VEventGenerator.hh"
#include "Target.hh"
-// NPLib
-#include "TInitialConditions.h"
-
// NPLib header
#include "NPReaction.h"
@@ -73,7 +70,7 @@ class EventGeneratorTransfert : public VEventGenerator
public: // Inherit from VEventGenerator class
void ReadConfiguration(string);
- void GenerateEvent(G4Event*, G4ParticleGun*);
+ void GenerateEvent(G4Event*);
void SetTarget(Target* Target) ;
@@ -100,11 +97,6 @@ class EventGeneratorTransfert : public VEventGenerator
double m_SigmaThetaX;
double m_SigmaPhiY;
-
- private: // TTree to store initial value of beam and reaction
- TInitialConditions* m_InitConditions;
-
-
// Other methods
void Print() const;
void InitializeRootOutput();
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index c7207bf29d40fee9e0d7df8433f3a05ae514ec0c..feb67d4bf4f8c6b7081886821ef9c24813df6437 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -33,7 +33,6 @@
// G4 headers
#include "G4ParticleTable.hh"
-#include "G4ParticleGun.hh"
#include "G4RotationMatrix.hh"
// G4 headers including CLHEP headers
@@ -43,7 +42,8 @@
// NPTool headers
#include "EventGeneratorTransfert.hh"
#include "RootOutput.h"
-
+#include "Particle.hh"
+#include "ParticleStack.hh"
using namespace CLHEP;
@@ -60,8 +60,7 @@ EventGeneratorTransfert::EventGeneratorTransfert()
m_SigmaX(0),
m_SigmaY(0),
m_SigmaThetaX(0),
- m_SigmaPhiY(0),
- m_InitConditions(new TInitialConditions)
+ m_SigmaPhiY(0)
{
//------------- Default Constructor -------------
}
@@ -80,7 +79,6 @@ void EventGeneratorTransfert::SetTarget(Target* Target)
EventGeneratorTransfert::~EventGeneratorTransfert()
{
//------------- Default Destructor ------------
- delete m_InitConditions;
delete m_Reaction;
}
@@ -126,9 +124,7 @@ EventGeneratorTransfert::EventGeneratorTransfert( string name1
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTransfert::InitializeRootOutput()
{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
+
}
@@ -357,7 +353,7 @@ void EventGeneratorTransfert::ReadConfiguration(string Path)
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
+void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent)
{
// If first time, write the DeDx table
if (anEvent->GetEventID() == 0) {
@@ -374,9 +370,6 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
}
}
- // Clear contents of Precedent event (now stored in ROOTOutput)
- m_InitConditions->Clear();
-
//////////////////////////////////////////////////
//////Define the kind of particle to shoot////////
//////////////////////////////////////////////////
@@ -385,14 +378,14 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
G4int LightA = m_Reaction->GetNucleus3()->GetA() ;
G4ParticleDefinition* LightName
- = G4ParticleTable::GetParticleTable()->GetIon(LightZ, LightA, 0.);
+ = G4ParticleTable::GetParticleTable()->GetIon(LightZ, LightA, m_Reaction->GetExcitationLight()*MeV);
// Recoil
G4int HeavyZ = m_Reaction->GetNucleus4()->GetZ() ;
G4int HeavyA = m_Reaction->GetNucleus4()->GetA() ;
G4ParticleDefinition* HeavyName
- = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitation()*MeV);
+ = G4ParticleTable::GetParticleTable()->GetIon(HeavyZ, HeavyA, m_Reaction->GetExcitationHeavy()*MeV);
// Beam
G4int BeamZ = m_Reaction->GetNucleus1()->GetZ();
@@ -404,141 +397,119 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
///// of interaction in target and Energy Loss of the beam within /////
///// the target. /////
///////////////////////////////////////////////////////////////////////
- G4ThreeVector InterCoord;
-
- G4double Beam_thetaX = 0, Beam_phiY = 0;
- G4double Beam_theta = 0, Beam_phi = 0;
- G4double FinalBeamEnergy = 0 ;
- G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
-
- m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
+ G4ThreeVector InterCoord;
+
+ G4double Beam_thetaX = 0, Beam_phiY = 0;
+ G4double Beam_theta = 0, Beam_phi = 0;
+ G4double FinalBeamEnergy = 0 ;
+ G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
+
+ m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
0, m_SigmaY, 0, m_SigmaPhiY,
InitialBeamEnergy,
BeamName,
InterCoord, Beam_thetaX, Beam_phiY,
Beam_theta, Beam_phi,
FinalBeamEnergy);
-
- m_Reaction->SetBeamEnergy(FinalBeamEnergy);
- m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
-
- // write vertex position to ROOT file
- G4double x0 = InterCoord.x();
- G4double y0 = InterCoord.y();
- G4double z0 = InterCoord.z();
- m_InitConditions->SetICPositionX(x0);//
- m_InitConditions->SetICPositionY(y0);//
- m_InitConditions->SetICPositionZ(z0);//
-
- // write emittance angles to ROOT file
- m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
- m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
-
- // write angles to ROOT file
- m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
- m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
-
- //////////////////////////////////////////////////////////
- ///// Build rotation matrix to go from the incident //////
- ///// beam frame to the "world" frame //////
- //////////////////////////////////////////////////////////
- G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
- cos(Beam_theta) * sin(Beam_phi),
- -sin(Beam_theta));
- G4ThreeVector col2(-sin(Beam_phi),
- cos(Beam_phi),
- 0);
- G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
- sin(Beam_theta) * sin(Beam_phi),
- cos(Beam_theta));
- G4RotationMatrix BeamToWorld(col1, col2, col3);
-
- /////////////////////////////////////////////////////////////////
- ///// Angles for emitted particles following Cross Section //////
- ///// Angles are in the beam frame //////
- /////////////////////////////////////////////////////////////////
-
- // Angles
- RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
- G4double ThetaCM = (m_Reaction->GetCrossSectionAngleMin() + CrossSectionShoot.shoot() * (m_Reaction->GetCrossSectionAngleMax() - m_Reaction->GetCrossSectionAngleMin())) * deg;
- G4double phi = RandFlat::shoot() * 2*pi;
- // write angles to ROOT file
- m_InitConditions->SetICEmittedAngleThetaCM(ThetaCM / deg);
- m_InitConditions->SetICEmittedAnglePhiIncidentFrame(phi / deg);
-
- //////////////////////////////////////////////////
- ///// Momentum and angles from kinematics /////
- ///// Angles are in the beam frame /////
- //////////////////////////////////////////////////
- // Variable where to store results
- G4double ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy;
- // Set the Theta angle of reaction
- m_Reaction->SetThetaCM(ThetaCM);
- // Compute Kinematic using previously defined ThetaCM
- m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
- // Momentum in beam frame for light particle
- G4ThreeVector momentum_kineLight_beam(sin(ThetaLight) * cos(phi),
- sin(ThetaLight) * sin(phi),
- cos(ThetaLight));
- // Momentum in beam frame for heavy particle
- G4ThreeVector momentum_kineHeavy_beam(sin(ThetaHeavy) * cos(phi+pi),
- sin(ThetaHeavy) * sin(phi+pi),
- cos(ThetaHeavy));
-
- // write angles/energy to ROOT file
- m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(ThetaLight / deg);
- m_InitConditions->SetICEmittedEnergy(EnergyLight/MeV);
-
- //////////////////////////////////////////////////
- ///////// Set up everything for shooting /////////
- //////////////////////////////////////////////////
- if (m_ShootLight) { // Case of light particle
- // Particle type
- particleGun->SetParticleDefinition(LightName);
- // Particle energy
- particleGun->SetParticleEnergy(EnergyLight);
- // Particle vertex position
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
- // Particle direction
- // Kinematical angles in the beam frame are transformed
- // to the "world" frame
- G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineLight_beam;
- // get theta and phi in the world frame
- G4double theta_world = momentum_kine_world.theta();
- G4double phi_world = momentum_kine_world.phi();
- if (phi_world < 1e-6) phi_world += 2*pi;
- // write angles in ROOT file
- m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
- m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
-
- //Set the gun to shoot
- particleGun->SetParticleMomentumDirection(momentum_kine_world);
- //Shoot the light particle
- particleGun->GeneratePrimaryVertex(anEvent);
- }
- if (m_ShootHeavy) { // Case of heavy particle
- // Particle type
- particleGun->SetParticleDefinition(HeavyName);
- // Particle energy
- particleGun->SetParticleEnergy(EnergyHeavy);
- // Particle vertex position
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
- // Particle direction
- // Kinematical angles in the beam frame are transformed
- // to the "world" frame
- G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineHeavy_beam;
- // get theta and phi in the world frame
- G4double theta_world = momentum_kine_world.theta();
- G4double phi_world = momentum_kine_world.phi();
- if (phi_world < 1e-6) phi_world += 2*pi;
- // write angles in ROOT file
- m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
- m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
- //Set the gun to shoot
- particleGun->SetParticleMomentumDirection(momentum_kine_world);
- //Shoot the light particle
- particleGun->GeneratePrimaryVertex(anEvent);
- }
+
+ m_Reaction->SetBeamEnergy(FinalBeamEnergy);
+
+
+ // write vertex position to ROOT file
+ G4double x0 = InterCoord.x();
+ G4double y0 = InterCoord.y();
+ G4double z0 = InterCoord.z();
+
+
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+ G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ cos(Beam_theta) * sin(Beam_phi),
+ -sin(Beam_theta));
+ G4ThreeVector col2(-sin(Beam_phi),
+ cos(Beam_phi),
+ 0);
+ G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ sin(Beam_theta) * sin(Beam_phi),
+ cos(Beam_theta));
+ G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+
+ // Angles
+ RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
+ G4double ThetaCM = (m_Reaction->GetCrossSectionAngleMin() + CrossSectionShoot.shoot() * (m_Reaction->GetCrossSectionAngleMax() - m_Reaction->GetCrossSectionAngleMin())) * deg;
+ G4double phi = RandFlat::shoot() * 2*pi;
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ G4double ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy;
+ // Set the Theta angle of reaction
+ m_Reaction->SetThetaCM(ThetaCM);
+ // Compute Kinematic using previously defined ThetaCM
+ m_Reaction->KineRelativistic(ThetaLight, EnergyLight, ThetaHeavy, EnergyHeavy);
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kineLight_beam(sin(ThetaLight) * cos(phi),
+ sin(ThetaLight) * sin(phi),
+ cos(ThetaLight));
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kineHeavy_beam(sin(ThetaHeavy) * cos(phi+pi),
+ sin(ThetaHeavy) * sin(phi+pi),
+ cos(ThetaHeavy));
+
+ //////////////////////////////////////////////////
+ ///////// Set up everything for shooting /////////
+ //////////////////////////////////////////////////
+ // Case of light particle
+ // Instantiate a new particle
+ Particle LightParticle;
+
+ // Particle type
+ LightParticle.SetParticleDefinition(LightName);
+ // Particle energy
+ LightParticle.SetParticleKineticEnergy(EnergyLight);
+ // Particle vertex position
+ LightParticle.SetParticlePosition(G4ThreeVector(x0, y0, z0));
+ // Particle direction
+ // Kinematical angles in the beam frame are transformed
+ // to the "world" frame
+ G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineLight_beam;
+ //Set the Momentum Direction
+ LightParticle.SetParticleMomentumDirection(momentum_kine_world);
+ // Set the shoot status
+ LightParticle.SetShootStatus(m_ShootLight) ;
+ //Add the particle to the particle stack
+ ParticleStack::getInstance()->AddParticleToStack(LightParticle);
+
+
+ // Case of heavy particle
+ // Instantiate a new particle
+ Particle HeavyParticle;
+
+ // Particle type
+ HeavyParticle.SetParticleDefinition(HeavyName);
+ // Particle energy
+ HeavyParticle.SetParticleKineticEnergy(EnergyHeavy);
+ // Particle vertex position
+ HeavyParticle.SetParticlePosition(G4ThreeVector(x0, y0, z0));
+ // Particle direction
+ // Kinematical angles in the beam frame are transformed
+ // to the "world" frame
+ momentum_kine_world = BeamToWorld * momentum_kineHeavy_beam;
+ //Set the Momentum Direction
+ HeavyParticle.SetParticleMomentumDirection(momentum_kine_world);
+ // Set the shoot status
+ HeavyParticle.SetShootStatus(m_ShootHeavy) ;
+ //Add the particle to the particle stack
+ ParticleStack::getInstance()->AddParticleToStack(HeavyParticle);
}
diff --git a/NPSimulation/src/PhysicsList.cc b/NPSimulation/src/PhysicsList.cc
index 16f2cbb2e39056cd16d69f60f2eb17a5941b93be..35ed164e77d9789fa7cca7b8719e28cfb59088c0 100644
--- a/NPSimulation/src/PhysicsList.cc
+++ b/NPSimulation/src/PhysicsList.cc
@@ -1,348 +1,349 @@
-/*****************************************************************************
- * Copyright (C) 2009-2010 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
- * *
- * Creation Date : January 2009 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: *
- * A quite standard, non-modulable Geant4 PPhysicis list. *
- * Well suited for low energy ions physics. *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * A good improvement should be a modular physicis list in order to deal *
- * accuratly with different physics cases. *
- *****************************************************************************/
- #include "PhysicsList.hh"
-
-// I/O
-#include "G4ios.hh"
-#include <iomanip>
-
-//Particle Definition
-#include "G4ParticleTypes.hh"
-#include "G4IonConstructor.hh"
-#include "G4ParticleDefinition.hh"
-#include "G4ParticleTable.hh"
-
-
-//Process
-#include "G4Transportation.hh"
-
-// Standard Geant4
-#include "G4ComptonScattering.hh"
-#include "G4GammaConversion.hh"
-#include "G4PhotoElectricEffect.hh"
-#include "G4eMultipleScattering.hh"
-#include "G4eIonisation.hh"
-#include "G4eBremsstrahlung.hh"
-#include "G4eplusAnnihilation.hh"
-
-// Penelope
-#include "G4PenelopeIonisation.hh"
-#include "G4PenelopeBremsstrahlung.hh"
-#include "G4PenelopeAnnihilation.hh"
-
-#include "G4PenelopeCompton.hh"
-#include "G4PenelopeGammaConversion.hh"
-#include "G4PenelopePhotoElectric.hh"
-#include "G4PenelopeRayleigh.hh"
-
-// Low energy ~ Penelope
-#include "G4LowEnergyIonisation.hh"
-#include "G4LowEnergyBremsstrahlung.hh"
-
-#include "G4LowEnergyCompton.hh"
-#include "G4LowEnergyGammaConversion.hh"
-#include "G4LowEnergyPhotoElectric.hh"
-#include "G4LowEnergyRayleigh.hh"
-
-
-#include "G4MuIonisation.hh"
-#include "G4MuMultipleScattering.hh"
-#include "G4MuBremsstrahlung.hh"
-#include "G4MuPairProduction.hh"
-
-#include "G4hIonisation.hh"
-#include "G4ionIonisation.hh"
-#include "G4hMultipleScattering.hh"
-#include "G4hLowEnergyIonisation.hh"
-
-#include "G4EmProcessOptions.hh"
-#include "G4ProcessManager.hh"
-#include "G4ProcessVector.hh"
-
-// Cut
-#include "G4ParticleWithCuts.hh"
-#include "G4UserSpecialCuts.hh"
-
-// Decay
-#include "G4Decay.hh"
-#include "G4DecayTable.hh"
-#include "G4VDecayChannel.hh"
-#include "G4NuclearDecayChannel.hh"
-#include "G4BetaMinusDecayChannel.hh"
-
-
-//#include "G4Gamma.hh"
-//#include "G4Electron.hh"
-//#include "G4Positron.hh"
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PhysicsList::PhysicsList()
-{
- // ie: no secondaries
- defaultCutValue = 1000 * pc;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PhysicsList::~PhysicsList()
-{
- ;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructParticle()
-{
- // In this method, static member functions should be called
- // for all particles which you want to use.
- // This ensures that objects of these particle types will be
- // created in the program.
-
- //Usefull to test geometry
- G4Geantino::GeantinoDefinition();
-
- //Usefull for Gamma
- ConstructBosons();
-
- //Usefull for betaDecay
- ConstructLeptons();
-
- //Needed by G4 (4.9.2) to run on mac os X ;-)
- ConstructMesons();
-
- //usefull for p and n
- ConstructBaryons();
-
- //Usefull of course :p
- ConstructIons();
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructBosons()
-{
- G4Geantino::GeantinoDefinition() ;
- G4ChargedGeantino::ChargedGeantinoDefinition() ;
- G4Gamma::GammaDefinition() ;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructLeptons()
-{
- G4Electron::ElectronDefinition() ;
- G4Positron::PositronDefinition() ;
- G4NeutrinoE::NeutrinoEDefinition() ;
- G4AntiNeutrinoE::AntiNeutrinoEDefinition() ;
- //G4NeutrinoMu::NeutrinoMuDefinition() ;
- //G4AntiNeutrinoMu::AntiNeutrinoMuDefinition() ;
- //G4MuonPlus::MuonPlusDefinition() ;
- //G4MuonMinus::MuonMinusDefinition() ;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructBaryons()
-{
- G4Proton::ProtonDefinition() ;
- G4Neutron::NeutronDefinition();
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructMesons()
-{
- // mesons
- G4PionPlus ::PionPlusDefinition();
- G4PionMinus ::PionMinusDefinition();
- G4PionZero ::PionZeroDefinition();
- G4Eta ::EtaDefinition();
- G4EtaPrime ::EtaPrimeDefinition();
- // G4RhoZero ::RhoZeroDefinition();
- G4KaonPlus ::KaonPlusDefinition();
- G4KaonMinus ::KaonMinusDefinition();
- G4KaonZero ::KaonZeroDefinition();
- G4AntiKaonZero ::AntiKaonZeroDefinition();
- G4KaonZeroLong ::KaonZeroLongDefinition();
- G4KaonZeroShort::KaonZeroShortDefinition();
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructIons()
-{
- G4He3::He3Definition() ;
- G4Deuteron::DeuteronDefinition() ;
- G4Triton::TritonDefinition() ;
- G4Alpha::AlphaDefinition() ;
-
- G4IonConstructor iConstructor ;
- iConstructor.ConstructParticle() ;
-
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructProcess()
-{
- AddTransportation() ;
- ConstructEM() ;
-
- SetCuts() ;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructEM()
-{
-
- theParticleIterator->reset();
-
- while ((*theParticleIterator)()) {
- G4ParticleDefinition* particle = theParticleIterator->value() ;
- G4ProcessManager* pmanager = particle->GetProcessManager() ;
- G4String particleName = particle->GetParticleName() ;
-
- if (particleName == "gamma") {
- // gamma
- //standard Geant4
- pmanager->AddDiscreteProcess(new G4PhotoElectricEffect) ;
- pmanager->AddDiscreteProcess(new G4ComptonScattering) ;
- pmanager->AddDiscreteProcess(new G4GammaConversion) ;
- //Low energy
- //pmanager->AddDiscreteProcess(new G4LowEnergyPhotoElectric) ;
- //pmanager->AddDiscreteProcess(new G4LowEnergyCompton) ;
- //pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion) ;
- //pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh) ;
- // Penelope
- //pmanager->AddDiscreteProcess(new G4PenelopePhotoElectric) ;
- //pmanager->AddDiscreteProcess(new G4PenelopeCompton) ;
- //pmanager->AddDiscreteProcess(new G4PenelopeGammaConversion) ;
- //pmanager->AddDiscreteProcess(new G4PenelopeRayleigh) ;
-
- } else if (particleName == "e-") {
- //electron
- pmanager->AddProcess(new G4eMultipleScattering , -1, 1, 1) ;
- //standard geant4:
- pmanager->AddProcess(new G4eIonisation , -1, 2, 2) ;
- pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3) ;
- // Low energy:
- //pmanager->AddProcess(new G4LowEnergyIonisation , -1, 2, 2) ;
- //pmanager->AddProcess(new G4LowEnergyBremsstrahlung , -1, -1, 3) ;
- // Penelope:
- // pmanager->AddProcess(new G4PenelopeIonisation , -1, 2, 2) ;
- // pmanager->AddProcess(new G4PenelopeBremsstrahlung , -1, -1, 3) ;
-
-
- } else if (particleName == "e+") {
- //positron
- pmanager->AddProcess(new G4eMultipleScattering , -1, 1, 1 );
- // standard Geant4 and Low energy
- pmanager->AddProcess(new G4eIonisation , -1, 2, 2 );
- pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3 );
- pmanager->AddProcess(new G4eplusAnnihilation , 0, -1, 4 );
- //Penelope:
- //pmanager->AddProcess(new G4PenelopeIonisation , -1, 2, 2 );
- //pmanager->AddProcess(new G4PenelopeBremsstrahlung , -1, -1, 3 );
- //pmanager->AddProcess(new G4PenelopeAnnihilation , 0, -1, 4 );
-
-
-
- } else if (particleName == "mu+" ||
- particleName == "mu-") {
-
- //muon
- /* pmanager->AddProcess(new G4muMultipleScattering , -1, 1, 1 ) ;
- pmanager->AddProcess(new G4MuIonisation , -1, 2, 2 ) ;
- pmanager->AddProcess(new G4MuBremsstrahlung , -1, -1, 3 ) ;
- pmanager->AddProcess(new G4MuPairProduction , -1, -1, 4 ) ;*/
-
- } else if (particleName == "GenericIon") {
- pmanager->AddProcess(new G4hMultipleScattering(), -1, 1, 1) ;
- G4ionIonisation* iI = new G4ionIonisation ;
- // mod by Nicolas [07/05/09]
- iI->ActivateNuclearStopping(true) ;
- iI->ActivateStoppingData(true) ;
- pmanager->AddProcess(iI , -1, 2, 2) ;
-
- //all others charged particles except geantino
- } else if ((!particle->IsShortLived()) &&
- (particle->GetPDGCharge() != 0.0) &&
- (particleName != "chargedgeantino")) {
-
- G4hIonisation* hI = new G4hIonisation ;
- // mod by Nicolas [07/05/09]
-// hI->ActivateNuclearStopping(true) ;
- pmanager->AddProcess(new G4hMultipleScattering , -1, 1, 1) ;
- pmanager->AddProcess(hI , -1, 2, 2) ;
-
-
- }//end else if
- }//end while particle
-
- G4EmProcessOptions opt ;
- opt.SetSubCutoff(true) ;
- opt.SetMinEnergy(0.001*eV) ;
- opt.SetMaxEnergy(1000.*MeV) ;
- opt.SetDEDXBinning(1000) ;
- opt.SetLambdaBinning(1000) ;
- // mod by Nicolas [07/05/09]
-// opt.SetLinearLossLimit(1.e-3) ;
-
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::SetCuts()
-{
- // uppress error messages even in case e/gamma/proton do not exist
- G4int temp = GetVerboseLevel();
- SetVerboseLevel(0);
- // " G4VUserPhysicsList::SetCutsWithDefault" method sets
- // the default cut value for all particle types
- SetCutsWithDefault();
-
- // for gamma-rays
-// SetCutValue(0.01*mm, "gamma");
-// SetCutValue(0.01*mm, "e-");
-// SetCutValue(0.01*mm, "e+");
-
- // Retrieve verbose level
- SetVerboseLevel(temp);
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PhysicsList::ConstructDecay()
-{
-
-// Add Decay Process
- G4Decay* theDecayProcess = new G4Decay() ;
- theParticleIterator->reset() ;
- while ((*theParticleIterator)()) {
- G4ParticleDefinition* particle = theParticleIterator->value() ;
- G4ProcessManager* pmanager = particle->GetProcessManager();
- if (theDecayProcess->IsApplicable(*particle)) {
- pmanager ->AddProcess(theDecayProcess);
- // set ordering for PostStepDoIt and AtRestDoIt
- pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
- pmanager->SetProcessOrdering(theDecayProcess, idxAtRest);
- }
- }
-//end Add Decay Process
-
-}
-
-void PhysicsList::MyOwnConstruction()
-{
- ConstructDecay();
-}
-
-
-
+/*****************************************************************************
+ * Copyright (C) 2009-2010 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
+ * *
+ * Creation Date : January 2009 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * A quite standard, non-modulable Geant4 PPhysicis list. *
+ * Well suited for low energy ions physics. *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * A good improvement should be a modular physicis list in order to deal *
+ * accuratly with different physics cases. *
+ *****************************************************************************/
+ #include "PhysicsList.hh"
+
+// I/O
+#include "G4ios.hh"
+#include <iomanip>
+
+//Particle Definition
+#include "G4ParticleTypes.hh"
+#include "G4IonConstructor.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+
+
+//Process
+#include "G4Transportation.hh"
+
+// Standard Geant4
+#include "G4ComptonScattering.hh"
+#include "G4GammaConversion.hh"
+#include "G4PhotoElectricEffect.hh"
+#include "G4eMultipleScattering.hh"
+#include "G4eIonisation.hh"
+#include "G4eBremsstrahlung.hh"
+#include "G4eplusAnnihilation.hh"
+
+// Penelope
+//#include "G4PenelopeIonisation.hh"
+//#include "G4PenelopeBremsstrahlung.hh"
+//#include "G4PenelopeAnnihilation.hh"
+
+//#include "G4PenelopeCompton.hh"
+//#include "G4PenelopeGammaConversion.hh"
+//#include "G4PenelopePhotoElectric.hh"
+//#include "G4PenelopeRayleigh.hh"
+
+// Low energy ~ Penelope
+//#include "G4LowEnergyIonisation.hh"
+//#include "G4LowEnergyBremsstrahlung.hh"
+
+//#include "G4LowEnergyCompton.hh"
+//#include "G4LowEnergyGammaConversion.hh"
+//#include "G4LowEnergyPhotoElectric.hh"
+//#include "G4LowEnergyRayleigh.hh"
+
+
+#include "G4MuIonisation.hh"
+#include "G4MuMultipleScattering.hh"
+#include "G4MuBremsstrahlung.hh"
+#include "G4MuPairProduction.hh"
+
+#include "G4hIonisation.hh"
+#include "G4ionIonisation.hh"
+#include "G4hMultipleScattering.hh"
+//#include "G4hLowEnergyIonisation.hh"
+
+#include "G4EmProcessOptions.hh"
+#include "G4ProcessManager.hh"
+#include "G4ProcessVector.hh"
+
+// Cut
+#include "G4ParticleWithCuts.hh"
+#include "G4UserSpecialCuts.hh"
+
+// Decay
+#include "G4Decay.hh"
+#include "G4DecayTable.hh"
+#include "G4VDecayChannel.hh"
+#include "G4NuclearDecayChannel.hh"
+#include "G4BetaMinusDecayChannel.hh"
+
+
+//#include "G4Gamma.hh"
+//#include "G4Electron.hh"
+//#include "G4Positron.hh"
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PhysicsList::PhysicsList()
+{
+ // ie: no secondaries
+ defaultCutValue = 1000 * pc;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PhysicsList::~PhysicsList()
+{
+ ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructParticle()
+{
+ // In this method, static member functions should be called
+ // for all particles which you want to use.
+ // This ensures that objects of these particle types will be
+ // created in the program.
+
+ //Usefull to test geometry
+ G4Geantino::GeantinoDefinition();
+
+ //Usefull for Gamma
+ ConstructBosons();
+
+ //Usefull for betaDecay
+ ConstructLeptons();
+
+ //Needed by G4 (4.9.2) to run on mac os X ;-)
+ ConstructMesons();
+
+ //usefull for p and n
+ ConstructBaryons();
+
+ //Usefull of course :p
+ ConstructIons();
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructBosons()
+{
+ G4Geantino::GeantinoDefinition() ;
+ G4ChargedGeantino::ChargedGeantinoDefinition() ;
+ G4Gamma::GammaDefinition() ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructLeptons()
+{
+ G4Electron::ElectronDefinition() ;
+ G4Positron::PositronDefinition() ;
+ G4NeutrinoE::NeutrinoEDefinition() ;
+ G4AntiNeutrinoE::AntiNeutrinoEDefinition() ;
+ //G4NeutrinoMu::NeutrinoMuDefinition() ;
+ //G4AntiNeutrinoMu::AntiNeutrinoMuDefinition() ;
+ //G4MuonPlus::MuonPlusDefinition() ;
+ //G4MuonMinus::MuonMinusDefinition() ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructBaryons()
+{
+ G4Proton::ProtonDefinition() ;
+ G4AntiProton::AntiProtonDefinition() ;
+ G4Neutron::NeutronDefinition();
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructMesons()
+{
+ // mesons
+ G4PionPlus ::PionPlusDefinition();
+ G4PionMinus ::PionMinusDefinition();
+ G4PionZero ::PionZeroDefinition();
+ G4Eta ::EtaDefinition();
+ G4EtaPrime ::EtaPrimeDefinition();
+ // G4RhoZero ::RhoZeroDefinition();
+ G4KaonPlus ::KaonPlusDefinition();
+ G4KaonMinus ::KaonMinusDefinition();
+ G4KaonZero ::KaonZeroDefinition();
+ G4AntiKaonZero ::AntiKaonZeroDefinition();
+ G4KaonZeroLong ::KaonZeroLongDefinition();
+ G4KaonZeroShort::KaonZeroShortDefinition();
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructIons()
+{
+ G4He3::He3Definition() ;
+ G4Deuteron::DeuteronDefinition() ;
+ G4Triton::TritonDefinition() ;
+ G4Alpha::AlphaDefinition() ;
+
+ G4IonConstructor iConstructor ;
+ iConstructor.ConstructParticle() ;
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructProcess()
+{
+ AddTransportation() ;
+ ConstructEM() ;
+
+ SetCuts() ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructEM()
+{
+
+ theParticleIterator->reset();
+
+ while ((*theParticleIterator)()) {
+ G4ParticleDefinition* particle = theParticleIterator->value() ;
+ G4ProcessManager* pmanager = particle->GetProcessManager() ;
+ G4String particleName = particle->GetParticleName() ;
+
+ if (particleName == "gamma") {
+ // gamma
+ //standard Geant4
+ pmanager->AddDiscreteProcess(new G4PhotoElectricEffect) ;
+ pmanager->AddDiscreteProcess(new G4ComptonScattering) ;
+ pmanager->AddDiscreteProcess(new G4GammaConversion) ;
+ //Low energy
+ //pmanager->AddDiscreteProcess(new G4LowEnergyPhotoElectric) ;
+ //pmanager->AddDiscreteProcess(new G4LowEnergyCompton) ;
+ //pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion) ;
+ //pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh) ;
+ // Penelope
+ //pmanager->AddDiscreteProcess(new G4PenelopePhotoElectric) ;
+ //pmanager->AddDiscreteProcess(new G4PenelopeCompton) ;
+ //pmanager->AddDiscreteProcess(new G4PenelopeGammaConversion) ;
+ //pmanager->AddDiscreteProcess(new G4PenelopeRayleigh) ;
+
+ } else if (particleName == "e-") {
+ //electron
+ pmanager->AddProcess(new G4eMultipleScattering , -1, 1, 1) ;
+ //standard geant4:
+ pmanager->AddProcess(new G4eIonisation , -1, 2, 2) ;
+ pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3) ;
+ // Low energy:
+ //pmanager->AddProcess(new G4LowEnergyIonisation , -1, 2, 2) ;
+ //pmanager->AddProcess(new G4LowEnergyBremsstrahlung , -1, -1, 3) ;
+ // Penelope:
+ // pmanager->AddProcess(new G4PenelopeIonisation , -1, 2, 2) ;
+ // pmanager->AddProcess(new G4PenelopeBremsstrahlung , -1, -1, 3) ;
+
+
+ } else if (particleName == "e+") {
+ //positron
+ pmanager->AddProcess(new G4eMultipleScattering , -1, 1, 1 );
+ // standard Geant4 and Low energy
+ pmanager->AddProcess(new G4eIonisation , -1, 2, 2 );
+ pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3 );
+ pmanager->AddProcess(new G4eplusAnnihilation , 0, -1, 4 );
+ //Penelope:
+ //pmanager->AddProcess(new G4PenelopeIonisation , -1, 2, 2 );
+ //pmanager->AddProcess(new G4PenelopeBremsstrahlung , -1, -1, 3 );
+ //pmanager->AddProcess(new G4PenelopeAnnihilation , 0, -1, 4 );
+
+
+
+ } else if (particleName == "mu+" ||
+ particleName == "mu-") {
+
+ //muon
+ /* pmanager->AddProcess(new G4muMultipleScattering , -1, 1, 1 ) ;
+ pmanager->AddProcess(new G4MuIonisation , -1, 2, 2 ) ;
+ pmanager->AddProcess(new G4MuBremsstrahlung , -1, -1, 3 ) ;
+ pmanager->AddProcess(new G4MuPairProduction , -1, -1, 4 ) ;*/
+
+ } else if (particleName == "GenericIon") {
+ pmanager->AddProcess(new G4hMultipleScattering(), -1, 1, 1) ;
+ G4ionIonisation* iI = new G4ionIonisation ;
+ // mod by Nicolas [07/05/09]
+// iI->ActivateNuclearStopping(true) ;
+ iI->ActivateStoppingData(true) ;
+ pmanager->AddProcess(iI , -1, 2, 2) ;
+
+ //all others charged particles except geantino
+ } else if ((!particle->IsShortLived()) &&
+ (particle->GetPDGCharge() != 0.0) &&
+ (particleName != "chargedgeantino")) {
+
+ G4hIonisation* hI = new G4hIonisation ;
+ // mod by Nicolas [07/05/09]
+// hI->ActivateNuclearStopping(true) ;
+ pmanager->AddProcess(new G4hMultipleScattering , -1, 1, 1) ;
+ pmanager->AddProcess(hI , -1, 2, 2) ;
+
+
+ }//end else if
+ }//end while particle
+
+ G4EmProcessOptions opt ;
+ opt.SetSubCutoff(true) ;
+ opt.SetMinEnergy(0.001*eV) ;
+ opt.SetMaxEnergy(1000.*MeV) ;
+ opt.SetDEDXBinning(1000) ;
+ opt.SetLambdaBinning(1000) ;
+ // mod by Nicolas [07/05/09]
+// opt.SetLinearLossLimit(1.e-3) ;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::SetCuts()
+{
+ // uppress error messages even in case e/gamma/proton do not exist
+ G4int temp = GetVerboseLevel();
+ SetVerboseLevel(0);
+ // " G4VUserPhysicsList::SetCutsWithDefault" method sets
+ // the default cut value for all particle types
+ SetCutsWithDefault();
+
+ // for gamma-rays
+// SetCutValue(0.01*mm, "gamma");
+// SetCutValue(0.01*mm, "e-");
+// SetCutValue(0.01*mm, "e+");
+
+ // Retrieve verbose level
+ SetVerboseLevel(temp);
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PhysicsList::ConstructDecay()
+{
+
+// Add Decay Process
+ G4Decay* theDecayProcess = new G4Decay() ;
+ theParticleIterator->reset() ;
+ while ((*theParticleIterator)()) {
+ G4ParticleDefinition* particle = theParticleIterator->value() ;
+ G4ProcessManager* pmanager = particle->GetProcessManager();
+ if (theDecayProcess->IsApplicable(*particle)) {
+ pmanager ->AddProcess(theDecayProcess);
+ // set ordering for PostStepDoIt and AtRestDoIt
+ pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
+ pmanager->SetProcessOrdering(theDecayProcess, idxAtRest);
+ }
+ }
+//end Add Decay Process
+
+}
+
+void PhysicsList::MyOwnConstruction()
+{
+ ConstructDecay();
+}
+
+
+