* Progress on Phase space

NPLib/Physics/NPPhaseSpace.cxx

@@ -53,30 +53,34 @@ void NPL::PhaseSpace::ReadConfigurationFile(NPL::InputParser parser){
 
   vector<string> token = {"Beam","Target","Daughters","ExcitationEnergies","Fermi"};
   fDaughters.clear(); 
-
+  fMasses.clear();
   for(unsigned int i = 0 ; i < blocks.size() ; i++){
     if(blocks[i]->HasTokenList(token)){
       int v = NPOptionManager::getInstance()->GetVerboseLevel();
       NPOptionManager::getInstance()->SetVerboseLevel(0);
       fBeam.ReadConfigurationFile(parser);
       NPOptionManager::getInstance()->SetVerboseLevel(v);
-
       fBeamEnergy= fBeam.GetEnergy();
       fTarget= GetParticle(blocks[i]->GetString("Target"),parser);
+      fTargetLV=TLorentzVector(0,0,0,fTarget.Mass());
       vector<string> vDaughters= blocks[i]->GetVectorString("Daughters");
       fExcitations= blocks[i]->GetVectorDouble("ExcitationEnergies","MeV");
-
-      for(auto d : vDaughters){
+      for(auto d:vDaughters){
         fDaughters.push_back(GetParticle(d,parser));
+        fMasses.push_back(fDaughters.rend()->Mass());
       }
+
+     SetBeamLV(0,0,1,fBeam.GetGamma()*fBeam.Mass());
+
     }
 
     else{
       cout << "ERROR: check your input file formatting \033[0m" << endl;
       exit(1);
     }
+
+  }
   }
-}
 ////////////////////////////////////////////////////////////////////////////////
 Particle NPL::PhaseSpace::GetParticle(string name, NPL::InputParser parser){
   vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("DefineParticle",name);
@@ -109,3 +113,25 @@ Particle NPL::PhaseSpace::GetParticle(string name, NPL::InputParser parser){
 
   return (NPL::Particle());
 }
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::PhaseSpace::SetBeamLV(const TLorentzVector& LV){
+    fBeamLV=LV;
+    fTotalLV = fBeamLV+fTargetLV;
+    static std::string option;
+    option="";
+    if(fFermi)
+      option="Fermi";
+    return fPhaseSpace.SetDecay(fTotalLV,fMasses.size(),&fMasses[0],option.c_str());
+
+  }
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::PhaseSpace::SetBeamLV(double dirX,double dirY,double dirZ,double 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;
+     TLorentzVector beamLV(p,E);
+     return SetBeamLV(fBeamLV);
+  }

NPLib/Physics/NPPhaseSpace.h

@@ -63,13 +63,22 @@ namespace NPL{
 
     public:
       // Getters and Setters
-      void     SetBeamEnergy(const double& eBeam)      {fBeamEnergy = eBeam;     initializePrecomputeVariable();}
+      bool SetBeamLV(const TLorentzVector& LV);
+      bool SetBeamLV(double dirX,double dirY,double dirZ,double Kinetic);
       Particle* GetBeam(){return &fBeam;};
-
+      double Generate(){return fPhaseSpace.Generate();}
+      TLorentzVector* GetDecayLV(unsigned int& n){return fPhaseSpace.GetDecay(n);}
+      size_t GetDecaySize(){return fDaughters.size();}
+      Particle* GetParticle(unsigned int& n){return &fDaughters[n];}
+      double  GetExcitation(unsigned int& n){return fExcitations[n];}
+      
     private: // intern precompute variable
-      void initializePrecomputeVariable();
-      TLorentzVector fInitialEnergyImpulsion;
-      std::vector<double> fmasses;
+      void Init();
+      TLorentzVector fBeamLV;
+      TLorentzVector fTargetLV;
+      TLorentzVector fTotalLV;
+
+      std::vector<double> fMasses;
 
       ClassDef(PhaseSpace,0)
 

NPSimulation/Process/BeamReaction.cc

@@ -197,7 +197,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
   if (m_shoot && m_length < m_StepSize) {
     if(m_ReactionType==QFS){
       //if ( m_QFS.IsAllowed() ) {
-        return true;
+      return true;
       //}
       //else{
       //  m_shoot=false;
@@ -505,16 +505,16 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
     int j=0;
     m_QFS.SetIsAllowed(false);
     while(!m_QFS.IsAllowed()){
-        // Shoot internal momentum for the removed cluster
-        // if a momentum Sigma is given then shoot in 3 indep. Gaussians
-        // if input files are given for distributions use them instead
-        m_QFS.ShootInternalMomentum();
+      // Shoot internal momentum for the removed cluster
+      // if a momentum Sigma is given then shoot in 3 indep. Gaussians
+      // if input files are given for distributions use them instead
+      m_QFS.ShootInternalMomentum();
 
-        // Go from CM to Lab
-        m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
+      // Go from CM to Lab
+      m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
 
-        j++;
-        if(j>100) cout<< "ERROR: too many iteration and QFS kinematical conditions not allowed"<<endl;    
+      j++;
+      if(j>100) cout<< "ERROR: too many iteration and QFS kinematical conditions not allowed"<<endl;    
     }
 
     //---------------------------------------------------------
@@ -666,6 +666,29 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
 
   }// end QFS
 
+  //////////////////////////
+  //  Phase space  Case   //
+  //////////////////////////
+  else if(m_ReactionType==PhaseSpace){
+    // Prepare beam LV
+    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();
+
+        static G4IonTable* IonTable
+          = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+        G4ParticleDefinition* dName;
+
+        if (d_Z == 0 && d_A == 1) // neutron is special case
+          dName = G4Neutron::Definition();
+        else 
+          dName = IonTable->GetIon(d_Z, d_A);
+      }
+    }
+  }
+
   ////////////////////
   //  Fusion Case   //
   ////////////////////