diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 108d319c9c4bb3bb41100e8a1baf71934662dd38..c5602d3b886d7abf2862b61c003c5d83d619041f 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -83,6 +83,9 @@ QFS::QFS(){
     fTheta2VsTheta1 = 0;
     fPhi2VsPhi1 = 0;
 
+    fPerpMomentumHist = NULL;
+    fParMomentumHist = NULL;
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -153,6 +156,16 @@ void QFS::ReadConfigurationFile(NPL::InputParser parser){
           fshoot1 = blocks[i]->GetInt("ShootLight");
           fshoot2 = blocks[i]->GetInt("ShootLight");
       }
+      if(blocks[i]->HasToken("PerpMomentumPath")){
+          vector<string> file_perp = blocks[i]->GetVectorString("PerpMomentumPath");
+          TH1F* Perptemp = Read1DProfile(file_perp[0], file_perp[1]);
+          SetPerpMomentumHist(Perptemp);
+      }
+      if(blocks[i]->HasToken("ParMomentumPath")){
+          vector<string> file_par = blocks[i]->GetVectorString("ParMomentumPath");
+          TH1F* Partemp = Read1DProfile(file_par[0], file_par[1]);
+          SetParMomentumHist(Partemp);
+      }
   }
 
   cout << "\033[0m" ;
@@ -215,14 +228,10 @@ void QFS::CalculateVariables(){
     fEnergyImpulsionLab_A = TLorentzVector(0.,0.,PA,EA);
     
     //Internal momentum of removed cluster/nucleon
-    //gRandom->SetSeed(0);
-    //Pa.SetX(gRandom->Gaus(0.,fMomentumSigma));
-    //Pa.SetY(gRandom->Gaus(0.,fMomentumSigma));
-    //Pa.SetZ(gRandom->Gaus(0.,fMomentumSigma));
     Pa.SetX(fInternalMomentum.X());
     Pa.SetY(fInternalMomentum.Y());
     Pa.SetZ(fInternalMomentum.Z());
-
+    
     //Internal momentum of heavy recoil after removal
     PB.SetXYZ( (-Pa.X()) , (-Pa.Y()) , (-Pa.Z()) );
 
@@ -400,10 +409,39 @@ void QFS::Dump(){
 }
 ////////////////////////////////////////////////////////////////////////////////////////////
 
-double QFS::ShootRandomThetaCM(){
-  // TO DO //
-  double theta =0;
-  return theta;
+TVector3 QFS::ShootInternalMomentum(){
+
+  //Shoot a momentum (vector) for the internal cluster in the beam-at-rest frame
+  // (1) if only a width is provided: shoot in 3 independant Gaussian
+  // (2) if input histos are provided: use them instead of option (1)  
+  // Remark : if both width and input histos are provided only histos are considered 
+
+  TVector3 momentum = {0,0,0};
+  double  PerpMomentum =0;
+  double  ParMomentum =0;
+  double  angle_tmp =0;
+
+  momentum.SetX(gRandom->Gaus(0.,fMomentumSigma));
+  momentum.SetY(gRandom->Gaus(0.,fMomentumSigma));
+  momentum.SetZ(gRandom->Gaus(0.,fMomentumSigma));
+
+  if(fPerpMomentumHist){
+      PerpMomentum=fPerpMomentumHist->GetRandom();
+      angle_tmp = gRandom->Rndm()*2*M_PI;
+      momentum.SetX(PerpMomentum * TMath::Cos(angle_tmp));
+      momentum.SetY(PerpMomentum * TMath::Sin(angle_tmp));
+  }
+  if(fParMomentumHist){
+      ParMomentum=fParMomentumHist->GetRandom();
+      momentum.SetZ(PerpMomentum);
+  }
+
+  //cout << " Shooting Random Momentum: "  << endl;
+  //cout<<"Px:"<<momentum.X() << endl;
+  //cout<<"Py:"<<momentum.Y() << endl;
+  //cout<<"Pz:"<<momentum.Z() << endl;
+  SetInternalMomentum(momentum);
+  return momentum;
 }
 
 ////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index fde81d853e41ddecd681fb4d88b03697627ceebc..0016793dde93bf66fcd5f72c8b332f444be2b78d 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -81,6 +81,8 @@ namespace NPL{
     double   fExcitationB;             // Excitation energy in MeV of beam-like ejectile
     double   fMomentumSigma;           // Width of the momentum distribution  (sigma)             
     TVector3 fInternalMomentum;        // Internal momentum of the removed cluster            
+    TH1F*    fPerpMomentumHist;        // Perpendicular momentum distribution in beam at rest frame
+    TH1F*    fParMomentumHist;         // Parallel momentum distribution in beam at rest frame
     double   fisotropic;               
 
     TGraph* fTheta2VsTheta1;
@@ -98,7 +100,7 @@ namespace NPL{
     void CalculateVariables();
     void TestR3B();
     void KineRelativistic(double &ThetaLab1, double &PhiLab1, double &KineticEnergyLab1, double &ThetaLab2, double &PhiLab2, double &KineticEnergyLab2);
-    double ShootRandomThetaCM();
+    TVector3 ShootInternalMomentum();
     bool IsAllowed();
     void Dump();
 
@@ -165,6 +167,10 @@ namespace NPL{
     void SetPhiCM(const double& angle) {fPhiCM = angle;}
     void SetInternalMomentum(const TVector3& mom) {fInternalMomentum = mom;}
     void SetMomentumSigma(const double& sigma) {fMomentumSigma = sigma;}
+    void SetPerpMomentumHist  (TH1F*  PerpMomentumHist)
+        {delete fPerpMomentumHist; fPerpMomentumHist   = PerpMomentumHist;}
+    void SetParMomentumHist  (TH1F*  ParMomentumHist)
+        {delete fParMomentumHist; fParMomentumHist   = ParMomentumHist;}
 
     //GETTERS
     Nucleus*  GetNucleusA()               {return &fNucleiA;}