diff --git a/NPLib/Physics/NPPhaseSpace.cxx b/NPLib/Physics/NPPhaseSpace.cxx
index 02d1cd99fdfb3f9c9eb50a936158b9c41680872e..42ff53bc8f5193d027b5ea7659224192dc621727 100644
--- a/NPLib/Physics/NPPhaseSpace.cxx
+++ b/NPLib/Physics/NPPhaseSpace.cxx
@@ -65,13 +65,12 @@ void NPL::PhaseSpace::ReadConfigurationFile(NPL::InputParser parser){
fTargetLV=TLorentzVector(0,0,0,fTarget.Mass());
vector<string> vDaughters= blocks[i]->GetVectorString("Daughters");
fExcitations= blocks[i]->GetVectorDouble("ExcitationEnergies","MeV");
+ fFermi=blocks[i]->GetInt("Fermi");
for(auto d:vDaughters){
fDaughters.push_back(GetParticle(d,parser));
- fMasses.push_back(fDaughters.rend()->Mass());
+ fMasses.push_back(fDaughters.rbegin()->Mass());
}
-
- SetBeamLV(0,0,1,fBeam.GetGamma()*fBeam.Mass());
-
+ SetBeamLV(0,0,1,fBeamEnergy);
}
else{
@@ -121,17 +120,18 @@ bool NPL::PhaseSpace::SetBeamLV(const TLorentzVector& LV){
option="";
if(fFermi)
option="Fermi";
- return fPhaseSpace.SetDecay(fTotalLV,fMasses.size(),&fMasses[0],option.c_str());
+
+ return fPhaseSpace.SetDecay(fTotalLV,fMasses.size(),fMasses.data(),option.c_str());
}
////////////////////////////////////////////////////////////////////////////////
bool NPL::PhaseSpace::SetBeamLV(double dirX,double dirY,double dirZ,double Kinetic){
+ fBeamEnergy=Kinetic;
fBeam.SetKineticEnergy(Kinetic);
double E = fBeam.GetGamma()*fBeam.Mass();
double pc=sqrt(E*E-fBeam.Mass()*fBeam.Mass());
TVector3 p(dirX,dirY,dirZ);
- p.Unit();
- p*=pc;
+ p=pc*p.Unit();
TLorentzVector beamLV(p,E);
- return SetBeamLV(fBeamLV);
+ return SetBeamLV(beamLV);
}
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index a09f5fcccb4310b95c7a5c2633f4aeed02c2ca6e..67f2e57262b78fba39c689cc013d48dd239c2e13 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -47,7 +47,7 @@ NPS::BeamReaction::BeamReaction(G4String modelName, G4Region* envelope)
m_shoot=false;
m_rand=0;
m_Z=0;
-
+ m_event_weight=1;
ABLA = new G4AblaInterface();
}
@@ -108,9 +108,8 @@ void NPS::BeamReaction::ReadConfiguration() {
m_active = true;
m_ReactionConditions = new TReactionConditions();
AttachReactionConditions();
- if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
- RootOutput::getInstance()->GetTree()->Branch(
- "ReactionConditions", "TReactionConditions", &m_ReactionConditions);
+ if (!RootOutput::getInstance()->GetTree()->FindBranch("EventWeight"))
+ RootOutput::getInstance()->GetTree()->Branch("EventWeight", &m_event_weight);
}
// Fusion
@@ -671,10 +670,15 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
//////////////////////////
else if(m_ReactionType==PhaseSpace){
// Prepare beam LV
+ static int d_Z,d_A;
+ static double d_Ex;
if(m_PhaseSpace.SetBeamLV(pdirection.x(),pdirection.y(),pdirection.z(),energy)){
- for(unsigned int i = 0,size=m_PhaseSpace.GetDecaySize() ; i < size ; i++){
- int d_Z = m_QFS.GetParticle1()->GetZ();
- int d_A = m_QFS.GetParticle1()->GetA();
+ m_event_weight= m_PhaseSpace.Generate();
+ unsigned int size = m_PhaseSpace.GetDecaySize();
+ for(unsigned int i = 0; i < size ; i++){
+ d_Z = m_PhaseSpace.GetParticle(i)->GetZ();
+ d_A = m_PhaseSpace.GetParticle(i)->GetA();
+ d_Ex =m_PhaseSpace.GetExcitation(i);
static G4IonTable* IonTable
= G4ParticleTable::GetParticleTable()->GetIonTable();
@@ -684,7 +688,18 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
if (d_Z == 0 && d_A == 1) // neutron is special case
dName = G4Neutron::Definition();
else
- dName = IonTable->GetIon(d_Z, d_A);
+ dName = IonTable->GetIon(d_Z, d_A,d_Ex);
+
+ auto dLV = m_PhaseSpace.GetDecayLV(i);
+
+ G4ThreeVector dir(dLV->Px(),dLV->Py(),dLV->Pz());
+ dir=dir.unit();
+ m_PhaseSpace.GetParticle(i)->SetBeta(dLV->Beta());
+
+ // Emmit daughter
+ G4DynamicParticle particle(dName, dir,m_PhaseSpace.GetParticle(i)->GetEnergy());
+ fastStep.CreateSecondaryTrack(particle, localPosition, time);
+
}
}
}
@@ -693,14 +708,6 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
// Fusion Case //
////////////////////
else if(m_ReactionType==Fusion){
- // Set the end of the step conditions
- fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
- fastStep.SetPrimaryTrackFinalPosition(worldPosition);
- fastStep.SetTotalEnergyDeposited(0);
- fastStep.SetPrimaryTrackFinalTime(time);// FIXME
- fastStep.KillPrimaryTrack();
- fastStep.SetPrimaryTrackPathLength(0.0);
-
static G4IonTable* IonTable
= G4ParticleTable::GetParticleTable()->GetIonTable();
//////Define the kind of particle to shoot////////
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index eec08e3790d87f9cd5298a93faba03dabb4dc79a..de59daf728e9277ea2de1bc0fabbc71fc7d0a3ad 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -70,6 +70,7 @@ namespace NPS{
double m_rand;
double m_length;
int m_Parent_ID;
+ double m_event_weight;
double SlowDownBeam(const G4ParticleDefinition* Beam, double IncidentEnergy, double Thickness,G4Material* Material);
private:// specific for the simple case of fusion
diff --git a/Projects/e793s/Reaction/47Kdp_PhaseSpace.reaction b/Projects/e793s/Reaction/47Kdp_PhaseSpace.reaction
new file mode 100755
index 0000000000000000000000000000000000000000..118fe932c5c68ec081c252dc1e7c758ae9c39648
--- /dev/null
+++ b/Projects/e793s/Reaction/47Kdp_PhaseSpace.reaction
@@ -0,0 +1,26 @@
+%%%%%%%%%%%%%%%%%%%%%% E793S %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 47K
+ ExcitationEnergy= 0 MeV
+ Energy= 362 MeV
+ SigmaEnergy= 0. MeV
+ SigmaThetaX= 0.0 deg
+ SigmaPhiY= 0.0 deg
+ SigmaX= 2 mm
+ SigmaY= 2 mm
+ MeanThetaX= 0 deg
+ MeanPhiY= 0 deg
+ MeanX= 0 mm
+ MeanY= 0 mm
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+PhaseSpace
+ Beam= 47K
+ Target= 2H
+ Daughters= 1H 47K n
+ ExcitationEnergies= 0 0 0 MeV
+ Fermi= 1