* Inelastic break up now working

NPSimulation/Process/BeamReaction.cc

@@ -68,6 +68,8 @@ void NPS::BeamReaction::ReadConfiguration() {
 
   if (input.GetAllBlocksWithToken("TwoBodyReaction").size() > 0)
     m_ReactionType = TwoBody;
+  if (input.GetAllBlocksWithToken("InelasticBreakup").size() > 0)
+    m_ReactionType = InelasticBreakup;
   else if (input.GetAllBlocksWithToken("QFSReaction").size() > 0)
     m_ReactionType = QFS;
   else if (input.GetAllBlocksWithToken("PhaseSpace").size() > 0)
@@ -89,6 +91,20 @@ void NPS::BeamReaction::ReadConfiguration() {
     }
   }
 
+  // Inelastic Breakup
+  if (m_ReactionType == InelasticBreakup) {
+    m_InelasticBreakup.ReadConfigurationFile(input);
+    m_BeamName = NPL::ChangeNameToG4Standard(m_InelasticBreakup.GetParticleBeam()->GetName());
+    if (m_InelasticBreakup.GetParticleLight()->GetName() != "") {
+      m_active = true;
+      m_ReactionConditions = new TReactionConditions();
+      AttachReactionConditions();
+      if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+        RootOutput::getInstance()->GetTree()->Branch("ReactionConditions", "TReactionConditions",
+                                                     &m_ReactionConditions);
+    }
+  }
+
   // QFS
   else if (m_ReactionType == QFS) {
     m_QFS.ReadConfigurationFile(input);
@@ -209,6 +225,9 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
         std::cout << "Two body reaction not allowed" << std::endl;
       }
     }
+    else if (m_ReactionType == InelasticBreakup) {
+      return true;
+    }
     else if (m_ReactionType == Fusion) {
       return true;
     }
@@ -415,6 +434,169 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
     m_ReactionConditions->SetMomentumDirectionZ(momentum_kine3_world.z());
     m_ReactionConditions->SetMomentumDirectionZ(momentum_kine4_world.z());
 
+  } // end if TwoBodyReaction
+  ////////////////////////////////////////
+  // Inelastic Breakup Reaction Case /////
+  ////////////////////////////////////////
+  if (m_ReactionType == InelasticBreakup) {
+
+    static G4IonTable* IonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
+    //////Define the kind of particle to shoot////////
+    // Particle Light
+    G4ParticleDefinition* LightName;
+    int LightZ = m_InelasticBreakup.GetParticleLight()->GetZ();
+    int LightA = m_InelasticBreakup.GetParticleLight()->GetA();
+
+    if (LightZ == 0 && LightA == 1) {
+      LightName = G4Neutron::Definition();
+    }
+    else {
+      LightName = IonTable->GetIon(LightZ, LightA);
+    }
+
+    // Particle Heavy
+    G4int HeavyZ = m_InelasticBreakup.GetParticleHeavy()->GetZ();
+    G4int HeavyA = m_InelasticBreakup.GetParticleHeavy()->GetA();
+    G4ParticleDefinition* HeavyName = IonTable->GetIon(HeavyZ, HeavyA);
+
+    // Particle TargetLike
+    G4int TargetZ = m_InelasticBreakup.GetParticleTarget()->GetZ();
+    G4int TargetA = m_InelasticBreakup.GetParticleTarget()->GetA();
+    G4ParticleDefinition* TargetName = IonTable->GetIon(TargetZ, TargetA);
+
+    // Set the Energy of the reaction
+    m_InelasticBreakup.SetBeamEnergy(reac_energy);
+
+    double Beam_theta = pdirection.theta();
+    double Beam_phi = pdirection.phi();
+
+    ///////////////////////////
+    ///// Beam Parameters /////
+    ///////////////////////////
+    m_ReactionConditions->SetBeamParticleName(PrimaryTrack->GetParticleDefinition()->GetParticleName());
+
+    m_ReactionConditions->SetBeamReactionEnergy(reac_energy);
+    m_ReactionConditions->SetVertexPositionX(worldPosition.x());
+    m_ReactionConditions->SetVertexPositionY(worldPosition.y());
+    m_ReactionConditions->SetVertexPositionZ(worldPosition.z());
+
+    G4ThreeVector U(1, 0, 0);
+    G4ThreeVector V(0, 1, 0);
+    G4ThreeVector ZZ(0, 0, 1);
+    m_ReactionConditions->SetBeamEmittanceTheta(PrimaryTrack->GetMomentumDirection().theta() / deg);
+    m_ReactionConditions->SetBeamEmittancePhi(PrimaryTrack->GetMomentumDirection().phi() / deg);
+    m_ReactionConditions->SetBeamEmittanceThetaX(PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+    m_ReactionConditions->SetBeamEmittancePhiY(PrimaryTrack->GetMomentumDirection().angle(V) / 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                         //////
+    /////////////////////////////////////////////////////////////////
+    double phi = G4RandFlat::shoot() * 2. * pi;
+
+    //////////////////////////////////////////////////
+    /////  Momentum and angles from  kinematics  /////
+    /////  Angles are in the beam frame          /////
+    //////////////////////////////////////////////////
+    // Variable where to store results
+    double ThetaLight, PhiLight, EnergyLight, ThetaTarget, EnergyTarget, ThetaHeavy, EnergyHeavy;
+
+    // Compute Kinematic using previously defined ThetaCM
+    m_InelasticBreakup.GenerateEvent(ThetaLight, PhiLight, EnergyLight, ThetaTarget, EnergyTarget, 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 World frame //to go from the incident beam frame to the "world"
+    // frame
+    G4ThreeVector momentum_kineLight_world = momentum_kineLight_beam;
+    momentum_kineLight_world.rotate(Beam_theta,
+                                    V);            // rotation of Beam_theta on Y axis
+    momentum_kineLight_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+    // Momentum in beam frame for light particle
+    G4ThreeVector momentum_kineTarget_beam(sin(ThetaTarget) * cos(phi), sin(ThetaTarget) * sin(phi), cos(ThetaTarget));
+    // Momentum in World frame //to go from the incident beam frame to the "world"
+    // frame
+    G4ThreeVector momentum_kineTarget_world = momentum_kineTarget_beam;
+    momentum_kineTarget_world.rotate(Beam_theta,
+                                     V);            // rotation of Beam_theta on Y axis
+    momentum_kineTarget_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+    // Momentum in beam frame for heavy particle
+    G4ThreeVector momentum_kineHeavy_beam(sin(ThetaHeavy) * cos(phi + pi), sin(ThetaHeavy) * sin(phi + pi),
+                                          cos(ThetaHeavy));
+    // Momentum in World frame
+    G4ThreeVector momentum_kineHeavy_world = momentum_kineHeavy_beam;
+    momentum_kineHeavy_world.rotate(Beam_theta,
+                                    V);            // rotation of Beam_theta on Y axis
+    momentum_kineHeavy_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+    // Emitt secondary
+    G4DynamicParticle particleLight(LightName, momentum_kineLight_world, EnergyLight);
+    fastStep.CreateSecondaryTrack(particleLight, localPosition, time);
+
+    G4DynamicParticle particleTarget(LightName, momentum_kineTarget_world, EnergyTarget);
+    fastStep.CreateSecondaryTrack(particleTarget, localPosition, time);
+
+    G4DynamicParticle particleHeavy(HeavyName, momentum_kineHeavy_world, EnergyHeavy);
+    fastStep.CreateSecondaryTrack(particleHeavy, localPosition, time);
+
+    ///////////////////////////////////////
+    ///// Emitted particle Parameters /////
+    ///////////////////////////////////////
+    // Names 3 and 4//
+    m_ReactionConditions->SetParticleName(LightName->GetParticleName());
+    m_ReactionConditions->SetParticleName(TargetName->GetParticleName());
+    m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
+    // Angle 3 and 4 //
+    m_ReactionConditions->SetTheta(ThetaLight / deg);
+    m_ReactionConditions->SetTheta(ThetaTarget / deg);
+    m_ReactionConditions->SetTheta(ThetaHeavy / deg);
+
+    m_ReactionConditions->SetPhi(phi / deg);
+    if ((phi + pi) / deg > 360)
+      m_ReactionConditions->SetPhi((phi - pi) / deg);
+    else
+      m_ReactionConditions->SetPhi((phi + pi) / deg);
+
+    // Energy 3 and 4 //
+    m_ReactionConditions->SetKineticEnergy(EnergyLight);
+    m_ReactionConditions->SetKineticEnergy(EnergyTarget);
+    m_ReactionConditions->SetKineticEnergy(EnergyHeavy);
+    // ThetaCM and Ex//
+    m_ReactionConditions->SetThetaCM(m_InelasticBreakup.GetThetaCM() / deg);
+    m_ReactionConditions->SetExcitationEnergy3(0);
+    m_ReactionConditions->SetExcitationEnergy3(0);
+    m_ReactionConditions->SetExcitationEnergy4(0);
+    // Momuntum X 3 and 4 //
+    m_ReactionConditions->SetMomentumDirectionX(momentum_kineLight_world.x());
+    m_ReactionConditions->SetMomentumDirectionX(momentum_kineTarget_world.x());
+    m_ReactionConditions->SetMomentumDirectionX(momentum_kineHeavy_world.x());
+    // Momuntum Y 3 and 4 //
+    m_ReactionConditions->SetMomentumDirectionY(momentum_kineLight_world.y());
+    m_ReactionConditions->SetMomentumDirectionY(momentum_kineTarget_world.y());
+    m_ReactionConditions->SetMomentumDirectionY(momentum_kineHeavy_world.y());
+    // Momuntum Z 3 and 4 //
+    m_ReactionConditions->SetMomentumDirectionZ(momentum_kineLight_world.z());
+    m_ReactionConditions->SetMomentumDirectionZ(momentum_kineTarget_world.z());
+    m_ReactionConditions->SetMomentumDirectionZ(momentum_kineHeavy_world.z());
+
   } // end if TwoBodyReaction
 
   // QFS

NPSimulation/Process/BeamReaction.hh

@@ -28,6 +28,7 @@
 #include "G4AblaInterface.hh"
 #include "G4Fragment.hh"
 #include "G4VFastSimulationModel.hh"
+#include "NPInelasticBreakup.h"
 #include "NPPhaseSpace.h"
 #include "NPQFS.h"
 #include "NPReaction.h"
@@ -35,7 +36,7 @@
 #include "TReactionConditions.h"
 class G4VPhysicalVolume;
 namespace NPS {
-  enum ReactionType { TwoBody, QFS, PhaseSpace, Fusion };
+  enum ReactionType { TwoBody, InelasticBreakup, QFS, PhaseSpace, Fusion };
 
   class BeamReaction : public G4VFastSimulationModel {
    public:
@@ -51,6 +52,7 @@ namespace NPS {
 
    private:
     NPL::Reaction m_Reaction;
+    NPL::InelasticBreakup m_InelasticBreakup;
     NPL::QFS m_QFS;
     NPL::PhaseSpace m_PhaseSpace;
     string m_BeamName;