Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4NuMuNucleusNcModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // $Id: G4NuMuNucleusNcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
 27 //
 28 // Geant4 Header : G4NuMuNucleusNcModel
 29 //
 30 // Author : V.Grichine 12.2.19
 31 //  
 32 
 33 #include "G4NuMuNucleusNcModel.hh"
 34 #include "G4NeutrinoNucleusModel.hh" 
 35 
 36 // #include "G4NuMuResQX.hh" 
 37 
 38 #include "G4SystemOfUnits.hh"
 39 #include "G4ParticleTable.hh"
 40 #include "G4ParticleDefinition.hh"
 41 #include "G4IonTable.hh"
 42 #include "Randomize.hh"
 43 #include "G4RandomDirection.hh"
 44 
 45 // #include "G4Integrator.hh"
 46 #include "G4DataVector.hh"
 47 #include "G4PhysicsTable.hh"
 48 #include "G4KineticTrack.hh"
 49 #include "G4DecayKineticTracks.hh"
 50 #include "G4KineticTrackVector.hh"
 51 #include "G4Fragment.hh"
 52 #include "G4ReactionProductVector.hh"
 53 
 54 
 55 #include "G4NeutrinoMu.hh"
 56 #include "G4AntiNeutrinoMu.hh"
 57 #include "G4Nucleus.hh"
 58 #include "G4LorentzVector.hh"
 59 
 60 using namespace std;
 61 using namespace CLHEP;
 62 
 63 #ifdef G4MULTITHREADED
 64     G4Mutex G4NuMuNucleusNcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
 65 #endif     
 66 
 67 
 68 G4NuMuNucleusNcModel::G4NuMuNucleusNcModel(const G4String& name) 
 69   : G4NeutrinoNucleusModel(name)
 70 {
 71   SetMinEnergy( 0.0*GeV );
 72   SetMaxEnergy( 100.*TeV );
 73   SetMinEnergy(1.e-6*eV);
 74 
 75   theNuMu =  G4NeutrinoMu::NeutrinoMu();
 76   theANuMu =  G4AntiNeutrinoMu::AntiNeutrinoMu();
 77 
 78   fMnumu = 0.; 
 79   fData = fMaster = false;
 80   InitialiseModel();  
 81      
 82 }
 83 
 84 
 85 G4NuMuNucleusNcModel::~G4NuMuNucleusNcModel()
 86 {}
 87 
 88 
 89 void G4NuMuNucleusNcModel::ModelDescription(std::ostream& outFile) const
 90 {
 91 
 92     outFile << "G4NuMuNucleusNcModel is a neutrino-nucleus (neutral current) scattering\n"
 93             << "model which uses the standard model \n"
 94             << "transfer parameterization.  The model is fully relativistic\n";
 95 
 96 }
 97 
 98 /////////////////////////////////////////////////////////
 99 //
100 // Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
101 
102 void G4NuMuNucleusNcModel::InitialiseModel()
103 {
104   G4String pName  = "nu_mu";
105   
106   G4int nSize(0), i(0), j(0), k(0);
107 
108   if(!fData)
109   { 
110 #ifdef G4MULTITHREADED
111     G4MUTEXLOCK(&numuNucleusModel);
112     if(!fData)
113     { 
114 #endif     
115       fMaster = true;
116 #ifdef G4MULTITHREADED
117     }
118     G4MUTEXUNLOCK(&numuNucleusModel);
119 #endif
120   }
121 
122   if(fMaster)
123   {  
124     const char* path = G4FindDataDir("G4PARTICLEXSDATA");
125     std::ostringstream ost1, ost2, ost3, ost4;
126     ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraynckr";
127 
128     std::ifstream filein1( ost1.str().c_str() );
129 
130     // filein.open("$PARTICLEXSDATA/");
131 
132     filein1>>nSize;
133 
134     for( k = 0; k < fNbin; ++k )
135     {
136       for( i = 0; i <= fNbin; ++i )
137       {
138         filein1 >> fNuMuXarrayKR[k][i];
139         // G4cout<< fNuMuXarrayKR[k][i] << "  ";
140       }
141     }
142     // G4cout<<G4endl<<G4endl;
143 
144     ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrnckr";
145     std::ifstream  filein2( ost2.str().c_str() );
146 
147     filein2>>nSize;
148 
149     for( k = 0; k < fNbin; ++k )
150     {
151       for( i = 0; i < fNbin; ++i )
152       {
153         filein2 >> fNuMuXdistrKR[k][i];
154         // G4cout<< fNuMuXdistrKR[k][i] << "  ";
155       }
156     }
157     // G4cout<<G4endl<<G4endl;
158 
159     ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraynckr";
160     std::ifstream  filein3( ost3.str().c_str() );
161 
162     filein3>>nSize;
163 
164     for( k = 0; k < fNbin; ++k )
165     {
166       for( i = 0; i <= fNbin; ++i )
167       {
168         for( j = 0; j <= fNbin; ++j )
169         {
170           filein3 >> fNuMuQarrayKR[k][i][j];
171           // G4cout<< fNuMuQarrayKR[k][i][j] << "  ";
172         }
173       }
174     }
175     // G4cout<<G4endl<<G4endl;
176 
177     ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrnckr";
178     std::ifstream  filein4( ost4.str().c_str() );
179 
180     filein4>>nSize;
181 
182     for( k = 0; k < fNbin; ++k )
183     {
184       for( i = 0; i <= fNbin; ++i )
185       {
186         for( j = 0; j < fNbin; ++j )
187         {
188           filein4 >> fNuMuQdistrKR[k][i][j];
189           // G4cout<< fNuMuQdistrKR[k][i][j] << "  ";
190         }
191       }
192     }
193     fData = true;
194   }
195 }
196 
197 /////////////////////////////////////////////////////////
198 
199 G4bool G4NuMuNucleusNcModel::IsApplicable(const G4HadProjectile & aPart, 
200                   G4Nucleus & )
201 {
202   G4bool result  = false;
203   G4String pName = aPart.GetDefinition()->GetParticleName();
204   G4double energy = aPart.GetTotalEnergy();
205   
206   if(  pName == "nu_mu" // || pName == "anti_nu_mu"   ) 
207         &&
208         energy > fMinNuEnergy                                )
209   {
210     result = true;
211   }
212 
213   return result;
214 }
215 
216 /////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
217 //
218 //
219 
220 G4HadFinalState* G4NuMuNucleusNcModel::ApplyYourself(
221      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
222 {
223   theParticleChange.Clear();
224   fProton = f2p2h = fBreak = false;
225   const G4HadProjectile* aParticle = &aTrack;
226   G4double energy = aParticle->GetTotalEnergy();
227 
228   G4String pName  = aParticle->GetDefinition()->GetParticleName();
229 
230   if( energy < fMinNuEnergy ) 
231   {
232     theParticleChange.SetEnergyChange(energy);
233     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
234     return &theParticleChange;
235   }
236   SampleLVkr( aTrack, targetNucleus);
237 
238   if( fBreak == true || fEmu < fMnumu ) // ~5*10^-6
239   {
240     // G4cout<<"ni, ";
241     theParticleChange.SetEnergyChange(energy);
242     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
243     return &theParticleChange;
244   }
245 
246   // LVs of initial state
247 
248   G4LorentzVector lvp1 = aParticle->Get4Momentum();
249   G4LorentzVector lvt1( 0., 0., 0., fM1 );
250   G4double mPip = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
251 
252   // 1-pi by fQtransfer && nu-energy
253   G4LorentzVector lvpip1( 0., 0., 0., mPip );
254   G4LorentzVector lvsum, lv2, lvX;
255   G4ThreeVector eP;
256   G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.);
257   G4DynamicParticle* aLept = nullptr; // lepton lv
258 
259   G4int Z = targetNucleus.GetZ_asInt();
260   G4int A = targetNucleus.GetA_asInt();
261   G4double  mTarg = targetNucleus.AtomicMass(A,Z);
262   G4int pdgP(0), qB(0);
263   // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip;
264 
265   G4int iPi     = GetOnePionIndex(energy);
266   G4double p1pi = GetNuMuOnePionProb( iPi, energy);
267 
268   if( p1pi > G4UniformRand() && fCosTheta > 0.9  ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus
269   {
270     // lvsum = lvp1 + lvpip1;
271     lvsum = lvp1 + lvt1;
272     // cost = fCosThetaPi;
273     cost = fCosTheta;
274     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
275     phi  = G4UniformRand()*CLHEP::twopi;
276     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
277 
278     // muMom = sqrt(fEmuPi*fEmuPi-fMnumu*fMnumu);
279     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
280 
281     eP *= muMom;
282 
283     // lv2 = G4LorentzVector( eP, fEmuPi );
284     lv2 = G4LorentzVector( eP, fEmu );
285     lv2 = fLVl;
286 
287     lvX = lvsum - lv2;
288     lvX = fLVh;
289     massX2 = lvX.m2();
290     G4double massX = lvX.m();
291     G4double massR = fLVt.m();
292 
293     // if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
294     if ( massX2 <= fM1*fM1 ) // 9-3-20 vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
295       if ( lvX.e() <= fM1 ) // 9-3-20 vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
296     {
297       theParticleChange.SetEnergyChange(energy);
298       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
299       return &theParticleChange;
300     }
301     fW2 = massX2;
302 
303     if(  pName == "nu_mu" )         aLept = new G4DynamicParticle( theNuMu, lv2 );  
304     else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theANuMu,  lv2 );
305     else
306     {
307       theParticleChange.SetEnergyChange(energy);
308       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
309       return &theParticleChange;
310     } 
311  
312     pdgP = 111;
313 
314     G4double eCut; // = fMpi + 0.5*(fMpi*fMpi - massX2)/mTarg; // massX -> fMpi
315 
316     if( A > 1 )
317     {
318       eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR);
319       eCut /= 2.*massR;
320       eCut += massX;
321     }
322     else  eCut = fM1 + fMpi;
323 
324     if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) // 
325     {
326       CoherentPion( lvX, pdgP, targetNucleus);
327     }
328     else
329     {
330       theParticleChange.SetEnergyChange(energy);
331       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
332       return &theParticleChange;
333     } 
334     theParticleChange.AddSecondary( aLept, fSecID );
335 
336     return &theParticleChange;
337   }
338   else // lepton part in lab
339   { 
340     lvsum = lvp1 + lvt1;
341     cost = fCosTheta;
342     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
343     phi  = G4UniformRand()*CLHEP::twopi;
344     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
345 
346     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
347 
348     eP *= muMom;
349 
350     lv2 = G4LorentzVector( eP, fEmu );
351 
352     lvX = lvsum - lv2;
353 
354     massX2 = lvX.m2();
355 
356     if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
357     {
358       theParticleChange.SetEnergyChange(energy);
359       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
360       return &theParticleChange;
361     }
362     fW2 = massX2;
363 
364     aLept = new G4DynamicParticle( theNuMu, lv2 );  
365   
366     theParticleChange.AddSecondary( aLept, fSecID );
367   }
368 
369   // hadron part
370 
371   fRecoil  = nullptr;
372   fCascade = false;
373   fString  = false;
374   
375   if( A == 1 )
376   {
377     qB = 1;
378 
379     // if( G4UniformRand() > 0.1 ) //  > 0.9999 ) // > 0.0001 ) //
380     {
381       ClusterDecay( lvX, qB );
382     }
383     return &theParticleChange;
384   }
385   G4Nucleus recoil;
386   G4double rM(0.), ratio = G4double(Z)/G4double(A);
387 
388   if( ratio > G4UniformRand() ) // proton is excited
389   {
390     fProton = true;
391     recoil = G4Nucleus(A-1,Z-1);
392     fRecoil = &recoil;
393     rM = recoil.AtomicMass(A-1,Z-1);
394 
395     fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
396           + G4ParticleTable::GetParticleTable()->FindParticle(111)->GetPDGMass();
397   }
398   else // excited neutron
399   {
400     fProton = false;
401     recoil = G4Nucleus(A-1,Z);
402     fRecoil = &recoil;
403     rM = recoil.AtomicMass(A-1,Z);
404 
405     fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
406           + G4ParticleTable::GetParticleTable()->FindParticle(111)->GetPDGMass(); 
407   }
408   // G4int       index = GetEnergyIndex(energy);
409   G4int nepdg = aParticle->GetDefinition()->GetPDGEncoding();
410 
411   G4double qeTotRat; //  = GetNuMuQeTotRat(index, energy);
412   qeTotRat = CalculateQEratioA( Z, A, energy, nepdg);
413 
414   G4ThreeVector dX = (lvX.vect()).unit();
415   G4double eX   = lvX.e();  // excited nucleon
416   G4double mX   = sqrt(massX2);
417 
418   if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE
419   {  
420     fString = false;
421 
422     if( fProton ) 
423     {  
424       fPDGencoding = 2212;
425       fMr =  proton_mass_c2;
426       recoil = G4Nucleus(A-1,Z-1);
427       fRecoil = &recoil;
428       rM = recoil.AtomicMass(A-1,Z-1);
429     } 
430     else
431     {  
432       fPDGencoding = 2112;
433       fMr =   G4ParticleTable::GetParticleTable()->
434   FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV;
435       recoil = G4Nucleus(A-1,Z);
436       fRecoil = &recoil;
437       rM = recoil.AtomicMass(A-1,Z);
438     } 
439     G4double eTh = fMr+0.5*(fMr*fMr-mX*mX)/rM;
440 
441     if(eX <= eTh) // vmg, very rarely out of kinematics
442     {
443       theParticleChange.SetEnergyChange(energy);
444       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
445       return &theParticleChange;
446     } 
447     FinalBarion( lvX, 0, fPDGencoding ); // p(n)+deexcited recoil
448   }
449   else // if ( eX < 9500000.*GeV ) // < 25.*GeV) //  < 95.*GeV ) // < 2.5*GeV ) //cluster decay
450   {  
451     if     (  fProton && pName == "nu_mu" )      qB =  1;
452     else if( !fProton && pName == "nu_mu" )      qB =  0;
453 
454     ClusterDecay( lvX, qB );
455   }
456   return &theParticleChange;
457 }
458 
459 
460 /////////////////////////////////////////////////////////////////////
461 ////////////////////////////////////////////////////////////////////
462 ///////////////////////////////////////////////////////////////////
463 
464 /////////////////////////////////////////////////
465 //
466 // sample x, then Q2
467 
468 void G4NuMuNucleusNcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
469 {
470   fBreak = false;
471   G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100); 
472   G4int Z = targetNucleus.GetZ_asInt(); 
473   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
474   G4double cost(1.), sint(0.), phi(0.), muMom(0.); 
475   G4ThreeVector eP, bst;
476   const G4HadProjectile* aParticle = &aTrack;
477   G4LorentzVector lvp1 = aParticle->Get4Momentum();
478   nMom = NucleonMomentum( targetNucleus );
479 
480   if( A == 1 || nMom == 0. ) // hydrogen, no Fermi motion ???
481   {
482     fNuEnergy = aParticle->GetTotalEnergy();
483     iTer = 0;
484 
485     do
486     {
487       fXsample = SampleXkr(fNuEnergy);
488       fQtransfer = SampleQkr(fNuEnergy, fXsample);
489       fQ2 = fQtransfer*fQtransfer;
490 
491      if( fXsample > 0. )
492       {
493         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
494         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
495       }
496       else
497       {
498         fW2 = fM1*fM1;
499         fEmu = fNuEnergy;
500       }
501       e3 = fNuEnergy + fM1 - fEmu;
502 
503       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl; // vmg ~10^-5 for NC
504     
505       pMu2 = fEmu*fEmu - fMnumu*fMnumu;
506       pX2  = e3*e3 - fW2;
507 
508       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
509       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
510       iTer++;
511     }
512     while( ( abs(fCosTheta) > 1. || fEmu < fMnumu ) && iTer < iTerMax );
513 
514     if( iTer >= iTerMax ) { fBreak = true; return; }
515 
516     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
517     { 
518       G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
519       // fCosTheta = -1. + 2.*G4UniformRand(); 
520       if(fCosTheta < -1.) fCosTheta = -1.;
521       if(fCosTheta >  1.) fCosTheta =  1.;
522     }
523     // LVs
524 
525     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
526     G4LorentzVector lvsum = lvp1 + lvt1;
527 
528     cost = fCosTheta;
529     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
530     phi  = G4UniformRand()*CLHEP::twopi;
531     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
532     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
533     eP *= muMom;
534     fLVl = G4LorentzVector( eP, fEmu );
535 
536     fLVh = lvsum - fLVl;
537     fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
538   }
539   else // Fermi motion, Q2 in nucleon rest frame
540   {
541     G4ThreeVector nMomDir = nMom*G4RandomDirection();
542 
543     if( !f2p2h ) // 1p1h
544     {
545       G4Nucleus recoil(A-1,Z);
546       rM = sqrt( recoil.AtomicMass(A-1,Z)*recoil.AtomicMass(A-1,Z) + nMom*nMom );
547       hM = tM - rM;
548 
549       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
550       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) ); 
551     }
552     else // 2p2h
553     {
554       G4Nucleus recoil(A-2,Z-1);
555       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
556       hM = tM - rM;
557 
558       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
559       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) ); 
560     }
561     // G4cout<<hM<<", ";
562     // bst = fLVh.boostVector(); // 9-3-20
563 
564     // lvp1.boost(-bst); // 9-3-20 -> nucleon rest system, where Q2 transfer is ???
565 
566     fNuEnergy  = lvp1.e();
567     iTer = 0;
568 
569     do
570     {
571       fXsample = SampleXkr(fNuEnergy);
572       fQtransfer = SampleQkr(fNuEnergy, fXsample);
573       fQ2 = fQtransfer*fQtransfer;
574 
575       if( fXsample > 0. )
576       {
577         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
578         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
579       }
580       else
581       {
582         fW2 = fM1*fM1;
583         fEmu = fNuEnergy;
584       }
585 
586       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
587 
588       e3 = fNuEnergy + fM1 - fEmu;
589 
590       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
591     
592       pMu2 = fEmu*fEmu - fMnumu*fMnumu;
593       pX2  = e3*e3 - fW2;
594 
595       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
596       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
597       iTer++;
598     }
599     while( ( abs(fCosTheta) > 1. || fEmu < fMnumu ) && iTer < iTerMax );
600 
601     if( iTer >= iTerMax ) { fBreak = true; return; }
602 
603     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
604     { 
605       G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
606       // fCosTheta = -1. + 2.*G4UniformRand(); 
607       if(fCosTheta < -1.) fCosTheta = -1.;
608       if(fCosTheta >  1.) fCosTheta =  1.;
609     }
610     // LVs
611     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
612     G4LorentzVector lvsum = lvp1 + lvt1;
613 
614     cost = fCosTheta;
615     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
616     phi  = G4UniformRand()*CLHEP::twopi;
617     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
618     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
619     eP *= muMom;
620     fLVl = G4LorentzVector( eP, fEmu );
621     fLVh = lvsum - fLVl;
622     // back to lab system
623     // fLVl.boost(bst); // 9-3-20
624     // fLVh.boost(bst); // 9-3-20
625   }
626   //G4cout<<iTer<<", "<<fBreak<<"; ";
627 }
628 
629 //
630 //
631 ///////////////////////////
632