Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4NuElNucleusCcModel.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: G4NuElNucleusCcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
 27 //
 28 // Geant4 Header : G4NuElNucleusCcModel
 29 //
 30 // Author : V.Grichine 12.2.19
 31 //  
 32 
 33 #include <iostream>
 34 #include <fstream>
 35 #include <sstream>
 36 
 37 #include "G4NuElNucleusCcModel.hh"
 38 // #include "G4NuMuNuclCcDistrKR.hh" 
 39 
 40 // #include "G4NuMuResQX.hh" 
 41 
 42 #include "G4SystemOfUnits.hh"
 43 #include "G4ParticleTable.hh"
 44 #include "G4ParticleDefinition.hh"
 45 #include "G4IonTable.hh"
 46 #include "Randomize.hh"
 47 #include "G4RandomDirection.hh"
 48 // #include "G4Threading.hh"
 49 
 50 // #include "G4Integrator.hh"
 51 #include "G4DataVector.hh"
 52 #include "G4PhysicsTable.hh"
 53 /*
 54 #include "G4CascadeInterface.hh"
 55 // #include "G4BinaryCascade.hh"
 56 #include "G4TheoFSGenerator.hh"
 57 #include "G4LundStringFragmentation.hh"
 58 #include "G4ExcitedStringDecay.hh"
 59 #include "G4FTFModel.hh"
 60 // #include "G4BinaryCascade.hh"
 61 #include "G4HadFinalState.hh"
 62 #include "G4HadSecondary.hh"
 63 #include "G4HadronicInteractionRegistry.hh"
 64 // #include "G4INCLXXInterface.hh"
 65 #include "G4QGSModel.hh"
 66 #include "G4QGSMFragmentation.hh"
 67 #include "G4QGSParticipants.hh"
 68 */
 69 #include "G4KineticTrack.hh"
 70 #include "G4DecayKineticTracks.hh"
 71 #include "G4KineticTrackVector.hh"
 72 #include "G4Fragment.hh"
 73 #include "G4NucleiProperties.hh"
 74 #include "G4ReactionProductVector.hh"
 75 
 76 #include "G4GeneratorPrecompoundInterface.hh"
 77 #include "G4PreCompoundModel.hh"
 78 #include "G4ExcitationHandler.hh"
 79 
 80 
 81 #include "G4Electron.hh"
 82 // #include "G4MuonPlus.hh"
 83 #include "G4Nucleus.hh"
 84 #include "G4LorentzVector.hh"
 85 
 86 using namespace std;
 87 using namespace CLHEP;
 88 
 89 #ifdef G4MULTITHREADED
 90     G4Mutex G4NuElNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
 91 #endif     
 92 
 93 
 94 G4NuElNucleusCcModel::G4NuElNucleusCcModel(const G4String& name) 
 95   : G4NeutrinoNucleusModel(name)
 96 {
 97   theElectron = G4Electron::Electron();
 98   fData = fMaster = false;
 99   fMel = electron_mass_c2;
100   InitialiseModel();  
101 }
102 
103 
104 G4NuElNucleusCcModel::~G4NuElNucleusCcModel()
105 {}
106 
107 
108 void G4NuElNucleusCcModel::ModelDescription(std::ostream& outFile) const
109 {
110 
111     outFile << "G4NuElNucleusCcModel is a neutrino-nucleus (charge current)  scattering\n"
112             << "model which uses the standard model \n"
113             << "transfer parameterization.  The model is fully relativistic\n";
114 
115 }
116 
117 /////////////////////////////////////////////////////////
118 //
119 // Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
120 
121 void G4NuElNucleusCcModel::InitialiseModel()
122 {
123   G4String pName  = "nu_e";
124   
125   G4int nSize(0), i(0), j(0), k(0);
126 
127   if(!fData)
128   { 
129 #ifdef G4MULTITHREADED
130     G4MUTEXLOCK(&numuNucleusModel);
131     if(!fData)
132     { 
133 #endif     
134       fMaster = true;
135 #ifdef G4MULTITHREADED
136     }
137     G4MUTEXUNLOCK(&numuNucleusModel);
138 #endif
139   }
140   
141   if(fMaster)
142   {  
143     const char* path = G4FindDataDir("G4PARTICLEXSDATA");
144     std::ostringstream ost1, ost2, ost3, ost4;
145     ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraycckr";
146 
147     std::ifstream filein1( ost1.str().c_str() );
148 
149     // filein.open("$PARTICLEXSDATA/");
150 
151     filein1>>nSize;
152 
153     for( k = 0; k < fNbin; ++k )
154     {
155       for( i = 0; i <= fNbin; ++i )
156       {
157         filein1 >> fNuMuXarrayKR[k][i];
158         // G4cout<< fNuMuXarrayKR[k][i] << "  ";
159       }
160     }
161     // G4cout<<G4endl<<G4endl;
162 
163     ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrcckr";
164     std::ifstream  filein2( ost2.str().c_str() );
165 
166     filein2>>nSize;
167 
168     for( k = 0; k < fNbin; ++k )
169     {
170       for( i = 0; i < fNbin; ++i )
171       {
172         filein2 >> fNuMuXdistrKR[k][i];
173         // G4cout<< fNuMuXdistrKR[k][i] << "  ";
174       }
175     }
176     // G4cout<<G4endl<<G4endl;
177 
178     ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraycckr";
179     std::ifstream  filein3( ost3.str().c_str() );
180 
181     filein3>>nSize;
182 
183     for( k = 0; k < fNbin; ++k )
184     {
185       for( i = 0; i <= fNbin; ++i )
186       {
187         for( j = 0; j <= fNbin; ++j )
188         {
189           filein3 >> fNuMuQarrayKR[k][i][j];
190           // G4cout<< fNuMuQarrayKR[k][i][j] << "  ";
191         }
192       }
193     }
194     // G4cout<<G4endl<<G4endl;
195 
196     ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrcckr";
197     std::ifstream  filein4( ost4.str().c_str() );
198 
199     filein4>>nSize;
200 
201     for( k = 0; k < fNbin; ++k )
202     {
203       for( i = 0; i <= fNbin; ++i )
204       {
205         for( j = 0; j < fNbin; ++j )
206         {
207           filein4 >> fNuMuQdistrKR[k][i][j];
208           // G4cout<< fNuMuQdistrKR[k][i][j] << "  ";
209         }
210       }
211     }
212     fData = true;
213   }
214 }
215 
216 /////////////////////////////////////////////////////////
217 
218 G4bool G4NuElNucleusCcModel::IsApplicable(const G4HadProjectile & aPart, 
219                   G4Nucleus & )
220 {
221   G4bool result  = false;
222   G4String pName = aPart.GetDefinition()->GetParticleName();
223   G4double energy = aPart.GetTotalEnergy();
224   fMinNuEnergy = GetMinNuElEnergy();
225   
226   if(  pName == "nu_e"    
227         &&
228         energy > fMinNuEnergy                                )
229   {
230     result = true;
231   }
232 
233   return result;
234 }
235 
236 /////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
237 //
238 //
239 
240 G4HadFinalState* G4NuElNucleusCcModel::ApplyYourself(
241      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
242 {
243   theParticleChange.Clear();
244   fProton = f2p2h = fBreak = false;
245   fCascade = fString  = false;
246   fLVh = fLVl = fLVt = fLVcpi = G4LorentzVector(0.,0.,0.,0.);
247 
248   const G4HadProjectile* aParticle = &aTrack;
249   G4double energy = aParticle->GetTotalEnergy();
250 
251   G4String pName  = aParticle->GetDefinition()->GetParticleName();
252 
253   if( energy < fMinNuEnergy ) 
254   {
255     theParticleChange.SetEnergyChange(energy);
256     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
257     return &theParticleChange;
258   }
259 
260   SampleLVkr( aTrack, targetNucleus);
261 
262   if( fBreak == true || fEmu < fMel ) // ~5*10^-6
263   {
264     // G4cout<<"ni, ";
265     theParticleChange.SetEnergyChange(energy);
266     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
267     return &theParticleChange;
268   }
269 
270   // LVs of initial state
271 
272   G4LorentzVector lvp1 = aParticle->Get4Momentum();
273   G4LorentzVector lvt1( 0., 0., 0., fM1 );
274   G4double mPip = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
275 
276   // 1-pi by fQtransfer && nu-energy
277   G4LorentzVector lvpip1( 0., 0., 0., mPip );
278   G4LorentzVector lvsum, lv2, lvX;
279   G4ThreeVector eP;
280   G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.), massX(0.), massR(0.), eCut(0.);
281   G4DynamicParticle* aLept = nullptr; // lepton lv
282 
283   G4int Z = targetNucleus.GetZ_asInt();
284   G4int A = targetNucleus.GetA_asInt();
285   G4double  mTarg = targetNucleus.AtomicMass(A,Z);
286   G4int pdgP(0), qB(0);
287   // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip;
288 
289   G4int iPi     = GetOnePionIndex(energy);
290   G4double p1pi = GetNuMuOnePionProb( iPi, energy);
291 
292   if( p1pi > G4UniformRand()  && fCosTheta > 0.9  ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus
293   {
294     // lvsum = lvp1 + lvpip1;
295     lvsum = lvp1 + lvt1;
296     // cost = fCosThetaPi;
297     cost = fCosTheta;
298     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
299     phi  = G4UniformRand()*CLHEP::twopi;
300     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
301 
302     // muMom = sqrt(fEmuPi*fEmuPi-fMel*fMel);
303     muMom = sqrt(fEmu*fEmu-fMel*fMel);
304 
305     eP *= muMom;
306 
307     // lv2 = G4LorentzVector( eP, fEmuPi );
308     // lv2 = G4LorentzVector( eP, fEmu );
309     lv2 = fLVl;
310 
311     // lvX = lvsum - lv2;
312     lvX = fLVh;
313     massX2 = lvX.m2();
314     massX = lvX.m();
315     massR = fLVt.m();
316     
317     if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
318     {
319       fCascade = true;
320       theParticleChange.SetEnergyChange(energy);
321       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
322       return &theParticleChange;
323     }
324     fW2 = massX2;
325 
326     if(  pName == "nu_e" )         aLept = new G4DynamicParticle( theElectron, lv2 );  
327     else
328     {
329       theParticleChange.SetEnergyChange(energy);
330       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
331       return &theParticleChange;
332     }
333     if( pName == "nu_e" ) pdgP =  211;
334     // else                   pdgP = -211;
335     // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mTarg; // massX -> fMpi
336 
337     if( A > 1 )
338     {
339       eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR);
340       eCut /= 2.*massR;
341       eCut += massX;
342     }
343     else  eCut = fM1 + fMpi;
344 
345     if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) // 
346     {
347       CoherentPion( lvX, pdgP, targetNucleus);
348     }
349     else
350     {
351       fCascade = true;
352       theParticleChange.SetEnergyChange(energy);
353       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
354       return &theParticleChange;
355     } 
356     theParticleChange.AddSecondary( aLept, fSecID );
357 
358     return &theParticleChange;
359   }
360   else // lepton part in lab
361   { 
362     lvsum = lvp1 + lvt1;
363     cost = fCosTheta;
364     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
365     phi  = G4UniformRand()*CLHEP::twopi;
366     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
367 
368     muMom = sqrt(fEmu*fEmu-fMel*fMel);
369 
370     eP *= muMom;
371 
372     lv2 = G4LorentzVector( eP, fEmu );
373     lv2 = fLVl;
374     lvX = lvsum - lv2;
375     lvX = fLVh;
376     massX2 = lvX.m2();
377 
378     if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
379     {
380       fCascade = true;
381       theParticleChange.SetEnergyChange(energy);
382       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
383       return &theParticleChange;
384     }
385     fW2 = massX2;
386 
387     if(  pName == "nu_e" )         aLept = new G4DynamicParticle( theElectron, lv2 );  
388     else
389     {
390       theParticleChange.SetEnergyChange(energy);
391       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
392       return &theParticleChange;
393     }
394     theParticleChange.AddSecondary( aLept, fSecID );
395   }
396 
397   // hadron part
398 
399   fRecoil  = nullptr;
400   
401   if( A == 1 )
402   {
403     if( pName == "nu_e" ) qB = 2;
404     // else                   qB = 0;
405 
406     // if( G4UniformRand() > 0.1 ) //  > 0.9999 ) // > 0.0001 ) //
407     {
408       ClusterDecay( lvX, qB );
409     }
410     return &theParticleChange;
411   }
412     /*
413     // else
414     {
415       if( pName == "nu_mu" ) pdgP =  211;
416       else                   pdgP = -211;
417 
418 
419       if ( fQtransfer < 0.95*GeV ) // < 0.35*GeV ) //
420       {
421   if( lvX.m() > mSum ) CoherentPion( lvX, pdgP, targetNucleus);
422       }
423     }
424     return &theParticleChange;
425   }
426   */
427   G4Nucleus recoil;
428   G4double rM(0.), ratio = G4double(Z)/G4double(A);
429 
430   if( ratio > G4UniformRand() ) // proton is excited
431   {
432     fProton = true;
433     recoil = G4Nucleus(A-1,Z-1);
434     fRecoil = &recoil;
435     rM = recoil.AtomicMass(A-1,Z-1);
436 
437     if( pName == "nu_e" ) // (++) state -> p + pi+
438     { 
439       fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
440           + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
441     }
442     else // (0) state -> p + pi-, n + pi0
443     {
444       // fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
445       //     + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass();
446     } 
447   }
448   else // excited neutron
449   {
450     fProton = false;
451     recoil = G4Nucleus(A-1,Z);
452     fRecoil = &recoil;
453     rM = recoil.AtomicMass(A-1,Z);
454 
455     if( pName == "nu_e" ) // (+) state -> n + pi+
456     {      
457       fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
458           + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
459     }
460     else // (-) state -> n + pi-, // n + pi0
461     {
462       // fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
463       //     + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass();
464     } 
465   }
466   // G4int       index = GetEnergyIndex(energy);
467   G4int nepdg = aParticle->GetDefinition()->GetPDGEncoding();
468 
469   G4double qeTotRat; //  = GetNuMuQeTotRat(index, energy);
470   qeTotRat = CalculateQEratioA( Z, A, energy, nepdg);
471 
472   G4ThreeVector dX = (lvX.vect()).unit();
473   G4double eX   = lvX.e();  // excited nucleon
474   G4double mX   = sqrt(massX2);
475   // G4double pX   = sqrt( eX*eX - mX*mX );
476   // G4double sumE = eX + rM;
477 
478   if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE
479   {  
480     fString = false;
481 
482     if( fProton ) 
483     {  
484       fPDGencoding = 2212;
485       fMr =  proton_mass_c2;
486       recoil = G4Nucleus(A-1,Z-1);
487       fRecoil = &recoil;
488       rM = recoil.AtomicMass(A-1,Z-1);
489     } 
490     else 
491     {  
492       fPDGencoding = 2112;
493       fMr =   G4ParticleTable::GetParticleTable()->
494   FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV;
495       recoil = G4Nucleus(A-1,Z);
496       fRecoil = &recoil;
497       rM = recoil.AtomicMass(A-1,Z);
498     } 
499     // sumE = eX + rM;   
500     G4double eTh = fMr + 0.5*(fMr*fMr - mX*mX)/rM;
501 
502     if( eX <= eTh ) // vmg, very rarely out of kinematics
503     {
504       fString = true;
505       theParticleChange.SetEnergyChange(energy);
506       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
507       return &theParticleChange;
508     }
509     // FinalBarion( fLVh, 0, fPDGencoding ); // p(n)+deexcited recoil
510     FinalBarion( lvX, 0, fPDGencoding ); // p(n)+deexcited recoil
511   }
512   else // if ( eX < 9500000.*GeV ) // <  25.*GeV) // < 95.*GeV ) // < 2.5*GeV ) //cluster decay
513   {  
514     if     (  fProton && pName == "nu_e" )      qB =  2;
515     else if( !fProton && pName == "nu_e" )      qB =  1;
516 
517     ClusterDecay( lvX, qB );
518   }
519   return &theParticleChange;
520 }
521 
522 
523 /////////////////////////////////////////////////////////////////////
524 ////////////////////////////////////////////////////////////////////
525 ///////////////////////////////////////////////////////////////////
526 
527 /////////////////////////////////////////////////
528 //
529 // sample x, then Q2
530 
531 void G4NuElNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
532 {
533   fBreak = false;
534   G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100); 
535   G4int Z = targetNucleus.GetZ_asInt(); 
536   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
537   G4double Ex(0.), ei(0.), nm2(0.);
538   G4double cost(1.), sint(0.), phi(0.), muMom(0.); 
539   G4ThreeVector eP, bst;
540   const G4HadProjectile* aParticle = &aTrack;
541   G4LorentzVector lvp1 = aParticle->Get4Momentum();
542 
543   if( A == 1 ) // hydrogen, no Fermi motion ???
544   {
545     fNuEnergy = aParticle->GetTotalEnergy();
546     iTer = 0;
547 
548     do
549     {
550       fXsample = SampleXkr(fNuEnergy);
551       fQtransfer = SampleQkr(fNuEnergy, fXsample);
552       fQ2 = fQtransfer*fQtransfer;
553 
554      if( fXsample > 0. )
555       {
556         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
557         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
558       }
559       else
560       {
561         fW2 = fM1*fM1;
562         fEmu = fNuEnergy;
563       }
564       e3 = fNuEnergy + fM1 - fEmu;
565 
566       if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
567     
568       pMu2 = fEmu*fEmu - fMel*fMel;
569 
570       if(pMu2 < 0.) { fBreak = true; return; }
571 
572       pX2  = e3*e3 - fW2;
573 
574       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
575       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
576       iTer++;
577     }
578     while( ( abs(fCosTheta) > 1. || fEmu < fMel ) && iTer < iTerMax );
579 
580     if( iTer >= iTerMax ) { fBreak = true; return; }
581 
582     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
583     { 
584       G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
585       // fCosTheta = -1. + 2.*G4UniformRand(); 
586       if(fCosTheta < -1.) fCosTheta = -1.;
587       if(fCosTheta >  1.) fCosTheta =  1.;
588     }
589     // LVs
590 
591     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
592     G4LorentzVector lvsum = lvp1 + lvt1;
593 
594     cost = fCosTheta;
595     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
596     phi  = G4UniformRand()*CLHEP::twopi;
597     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
598     muMom = sqrt(fEmu*fEmu-fMel*fMel);
599     eP *= muMom;
600     fLVl = G4LorentzVector( eP, fEmu );
601 
602     fLVh = lvsum - fLVl;
603     fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
604   }
605   else // Fermi motion, Q2 in nucleon rest frame
606   {
607     G4Nucleus recoil1( A-1, Z );
608     rM = recoil1.AtomicMass(A-1,Z);   
609     do
610     {
611       // nMom = NucleonMomentumBR( targetNucleus ); // BR
612       nMom = GgSampleNM( targetNucleus ); // Gg
613       Ex = GetEx(A-1, fProton);
614       ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom );
615       //   ei = 0.5*( tM - s2M - 2*eX );
616     
617       nm2 = ei*ei - nMom*nMom;
618       iTer++;
619     }
620     while( nm2 < 0. && iTer < iTerMax ); 
621 
622     if( iTer >= iTerMax ) { fBreak = true; return; }
623     
624     G4ThreeVector nMomDir = nMom*G4RandomDirection();
625 
626     if( !f2p2h || A < 3 ) // 1p1h
627     {
628       // hM = tM - rM;
629 
630       fLVt = G4LorentzVector( -nMomDir, sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ) ); // rM ); //
631       fLVh = G4LorentzVector(  nMomDir, ei ); // hM); //
632     }
633     else // 2p2h
634     {
635       G4Nucleus recoil(A-2,Z-1);
636       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
637       hM = tM - rM;
638 
639       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
640       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom )  ); 
641     }
642     // G4cout<<hM<<", ";
643     // bst = fLVh.boostVector();
644 
645     // lvp1.boost(-bst); // -> nucleon rest system, where Q2 transfer is ???
646 
647     fNuEnergy  = lvp1.e();
648     // G4double mN = fLVh.m(); // better mN = fM1 !? vmg
649     iTer = 0;
650 
651     do // no FM!?, 5.4.20 vmg
652     {
653       fXsample = SampleXkr(fNuEnergy);
654       fQtransfer = SampleQkr(fNuEnergy, fXsample);
655       fQ2 = fQtransfer*fQtransfer;
656 
657       // G4double mR = mN + fM1*(A-1.)*std::exp(-2.0*fQtransfer/mN); // recoil mass in+el
658 
659       if( fXsample > 0. )
660       {
661         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
662 
663         // fW2 = mN*mN - fQ2 + fQ2/fXsample; // sample excited hadron mass
664         // fEmu = fNuEnergy - fQ2/2./mR/fXsample; // fM1->mN
665 
666         fEmu = fNuEnergy - fQ2/2./fM1/fXsample; // fM1->mN
667       }
668       else
669       {
670         // fW2 = mN*mN;
671 
672         fW2 = fM1*fM1; 
673         fEmu = fNuEnergy;
674       }
675       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
676       // e3 = fNuEnergy + mR - fEmu;
677 
678       e3 = fNuEnergy + fM1 - fEmu;
679 
680       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
681     
682       pMu2 = fEmu*fEmu - fMel*fMel;
683       pX2  = e3*e3 - fW2;
684 
685       if(pMu2 < 0.) { fBreak = true; return; }
686 
687       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
688       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
689       iTer++;
690     }
691     while( ( abs(fCosTheta) > 1. || fEmu < fMel ) && iTer < iTerMax );
692 
693     if( iTer >= iTerMax ) { fBreak = true; return; }
694 
695     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
696     { 
697       G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
698       // fCosTheta = -1. + 2.*G4UniformRand(); 
699       if( fCosTheta < -1.) fCosTheta = -1.;
700       if( fCosTheta >  1.) fCosTheta =  1.;
701     }
702     // LVs
703     // G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., mN ); // fM1 );
704 
705     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 ); // fM1 );
706     G4LorentzVector lvsum = lvp1 + lvt1;
707 
708     cost = fCosTheta;
709     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
710     phi  = G4UniformRand()*CLHEP::twopi;
711     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
712     muMom = sqrt(fEmu*fEmu-fMel*fMel);
713     eP *= muMom;
714     fLVl = G4LorentzVector( eP, fEmu );
715     fLVh = lvsum - fLVl;
716 
717     // if( fLVh.e() < mN || fLVh.m2() < 0.) { fBreak = true; return; }
718 
719     if( fLVh.e() < fM1 || fLVh.m2() < 0.) { fBreak = true; return; }
720 
721     // back to lab system
722 
723     // fLVl.boost(bst);
724     // fLVh.boost(bst);
725   }
726   //G4cout<<iTer<<", "<<fBreak<<"; ";
727 }
728 
729 //
730 //
731 ///////////////////////////
732