Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4NuMuNucleusCcModel.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: G4NuMuNucleusCcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
 27 //
 28 // Geant4 Header : G4NuMuNucleusCcModel
 29 //
 30 // Author : V.Grichine 12.2.19
 31 //  
 32 
 33 #include <iostream>
 34 #include <fstream>
 35 #include <sstream>
 36 
 37 #include "G4NuMuNucleusCcModel.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 "G4MuonMinus.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 G4NuMuNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
 91 #endif     
 92 
 93 
 94 G4NuMuNucleusCcModel::G4NuMuNucleusCcModel(const G4String& name) 
 95   : G4NeutrinoNucleusModel(name)
 96 {
 97   fData = fMaster = false;
 98   InitialiseModel();  
 99 }
100 
101 
102 G4NuMuNucleusCcModel::~G4NuMuNucleusCcModel()
103 {}
104 
105 
106 void G4NuMuNucleusCcModel::ModelDescription(std::ostream& outFile) const
107 {
108 
109     outFile << "G4NuMuNucleusCcModel is a neutrino-nucleus (charge current)  scattering\n"
110             << "model which uses the standard model \n"
111             << "transfer parameterization.  The model is fully relativistic\n";
112 
113 }
114 
115 /////////////////////////////////////////////////////////
116 //
117 // Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
118 
119 void G4NuMuNucleusCcModel::InitialiseModel()
120 {
121   G4String pName  = "nu_mu";
122   // G4String pName  = "anti_nu_mu";
123   
124   G4int nSize(0), i(0), j(0), k(0);
125 
126   if(!fData)
127   { 
128 #ifdef G4MULTITHREADED
129     G4MUTEXLOCK(&numuNucleusModel);
130     if(!fData)
131     { 
132 #endif     
133       fMaster = true;
134 #ifdef G4MULTITHREADED
135     }
136     G4MUTEXUNLOCK(&numuNucleusModel);
137 #endif
138   }
139   
140   if(fMaster)
141   {  
142     const char* path = G4FindDataDir("G4PARTICLEXSDATA");
143     std::ostringstream ost1, ost2, ost3, ost4;
144     ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraycckr";
145 
146     std::ifstream filein1( ost1.str().c_str() );
147 
148     // filein.open("$PARTICLEXSDATA/");
149 
150     filein1>>nSize;
151 
152     for( k = 0; k < fNbin; ++k )
153     {
154       for( i = 0; i <= fNbin; ++i )
155       {
156         filein1 >> fNuMuXarrayKR[k][i];
157         // G4cout<< fNuMuXarrayKR[k][i] << "  ";
158       }
159     }
160     // G4cout<<G4endl<<G4endl;
161 
162     ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrcckr";
163     std::ifstream  filein2( ost2.str().c_str() );
164 
165     filein2>>nSize;
166 
167     for( k = 0; k < fNbin; ++k )
168     {
169       for( i = 0; i < fNbin; ++i )
170       {
171         filein2 >> fNuMuXdistrKR[k][i];
172         // G4cout<< fNuMuXdistrKR[k][i] << "  ";
173       }
174     }
175     // G4cout<<G4endl<<G4endl;
176 
177     ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraycckr";
178     std::ifstream  filein3( ost3.str().c_str() );
179 
180     filein3>>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           filein3 >> fNuMuQarrayKR[k][i][j];
189           // G4cout<< fNuMuQarrayKR[k][i][j] << "  ";
190         }
191       }
192     }
193     // G4cout<<G4endl<<G4endl;
194 
195     ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrcckr";
196     std::ifstream  filein4( ost4.str().c_str() );
197 
198     filein4>>nSize;
199 
200     for( k = 0; k < fNbin; ++k )
201     {
202       for( i = 0; i <= fNbin; ++i )
203       {
204         for( j = 0; j < fNbin; ++j )
205         {
206           filein4 >> fNuMuQdistrKR[k][i][j];
207           // G4cout<< fNuMuQdistrKR[k][i][j] << "  ";
208         }
209       }
210     }
211     fData = true;
212   }
213 }
214 
215 /////////////////////////////////////////////////////////
216 
217 G4bool G4NuMuNucleusCcModel::IsApplicable(const G4HadProjectile & aPart, 
218                   G4Nucleus & )
219 {
220   G4bool result  = false;
221   G4String pName = aPart.GetDefinition()->GetParticleName();
222   G4double energy = aPart.GetTotalEnergy();
223   
224   if(  pName == "nu_mu"  // || pName == "anti_nu_mu" )  
225         &&
226         energy > fMinNuEnergy                                )
227   {
228     result = true;
229   }
230 
231   return result;
232 }
233 
234 /////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
235 //
236 //
237 
238 G4HadFinalState* G4NuMuNucleusCcModel::ApplyYourself(
239      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
240 {
241   theParticleChange.Clear();
242   fProton = f2p2h = fBreak = false;
243   fCascade = fString  = false;
244   fLVh = fLVl = fLVt = fLVcpi = G4LorentzVector(0.,0.,0.,0.);
245 
246   const G4HadProjectile* aParticle = &aTrack;
247   G4double energy = aParticle->GetTotalEnergy();
248 
249   G4String pName  = aParticle->GetDefinition()->GetParticleName();
250 
251   if( energy < fMinNuEnergy ) 
252   {
253     theParticleChange.SetEnergyChange(energy);
254     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
255     return &theParticleChange;
256   }
257 
258   SampleLVkr( aTrack, targetNucleus);
259 
260   if( fBreak == true || fEmu < fMu ) // ~5*10^-6
261   {
262     // G4cout<<"ni, ";
263     theParticleChange.SetEnergyChange(energy);
264     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
265     return &theParticleChange;
266   }
267 
268   // LVs of initial state
269 
270   G4LorentzVector lvp1 = aParticle->Get4Momentum();
271   G4LorentzVector lvt1( 0., 0., 0., fM1 );
272   G4double mPip = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
273 
274   // 1-pi by fQtransfer && nu-energy
275   G4LorentzVector lvpip1( 0., 0., 0., mPip );
276   G4LorentzVector lvsum, lv2, lvX;
277   G4ThreeVector eP;
278   G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.), massX(0.), massR(0.), eCut(0.);
279   G4DynamicParticle* aLept = nullptr; // lepton lv
280 
281   G4int Z = targetNucleus.GetZ_asInt();
282   G4int A = targetNucleus.GetA_asInt();
283   G4double  mTarg = targetNucleus.AtomicMass(A,Z);
284   G4int pdgP(0), qB(0);
285   // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip;
286 
287   G4int iPi     = GetOnePionIndex(energy);
288   G4double p1pi = GetNuMuOnePionProb( iPi, energy);
289 
290   if( p1pi > G4UniformRand()  && fCosTheta > 0.9  ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus
291   {
292     // lvsum = lvp1 + lvpip1;
293     lvsum = lvp1 + lvt1;
294     // cost = fCosThetaPi;
295     cost = fCosTheta;
296     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
297     phi  = G4UniformRand()*CLHEP::twopi;
298     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
299 
300     // muMom = sqrt(fEmuPi*fEmuPi-fMu*fMu);
301     muMom = sqrt(fEmu*fEmu-fMu*fMu);
302 
303     eP *= muMom;
304 
305     // lv2 = G4LorentzVector( eP, fEmuPi );
306     // lv2 = G4LorentzVector( eP, fEmu );
307     lv2 = fLVl;
308 
309     // lvX = lvsum - lv2;
310     lvX = fLVh;
311     massX2 = lvX.m2();
312     massX = lvX.m();
313     massR = fLVt.m();
314     
315     if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
316     {
317       fCascade = true;
318       theParticleChange.SetEnergyChange(energy);
319       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
320       return &theParticleChange;
321     }
322     fW2 = massX2;
323 
324     if(  pName == "nu_mu" )         aLept = new G4DynamicParticle( theMuonMinus, lv2 );  
325     // else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus,  lv2 );
326     else
327     {
328       theParticleChange.SetEnergyChange(energy);
329       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
330       return &theParticleChange;
331     }
332     if( pName == "nu_mu" ) pdgP =  211;
333     // else                   pdgP = -211;
334     // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mTarg; // massX -> fMpi
335 
336     if( A > 1 )
337     {
338       eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR);
339       eCut /= 2.*massR;
340       eCut += massX;
341     }
342     else  eCut = fM1 + fMpi;
343 
344     if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) // 
345     {
346       CoherentPion( lvX, pdgP, targetNucleus);
347     }
348     else
349     {
350       fCascade = true;
351       theParticleChange.SetEnergyChange(energy);
352       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
353       return &theParticleChange;
354     } 
355     theParticleChange.AddSecondary( aLept, fSecID );
356 
357     return &theParticleChange;
358   }
359   else // lepton part in lab
360   { 
361     lvsum = lvp1 + lvt1;
362     cost = fCosTheta;
363     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
364     phi  = G4UniformRand()*CLHEP::twopi;
365     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
366 
367     muMom = sqrt(fEmu*fEmu-fMu*fMu);
368 
369     eP *= muMom;
370 
371     lv2 = G4LorentzVector( eP, fEmu );
372     lv2 = fLVl;
373     lvX = lvsum - lv2;
374     lvX = fLVh;
375     massX2 = lvX.m2();
376 
377     if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
378     {
379       fCascade = true;
380       theParticleChange.SetEnergyChange(energy);
381       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
382       return &theParticleChange;
383     }
384     fW2 = massX2;
385 
386     if(  pName == "nu_mu" )         aLept = new G4DynamicParticle( theMuonMinus, lv2 );  
387     // else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus,  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_mu" ) 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_mu" ) // (++) 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_mu" ) // (+) 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 // if( pName == "anti_nu_mu" ) 
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_mu" )      qB =  2;
515     // else if(  fProton && pName == "anti_nu_mu" ) qB =  0;
516     else if( !fProton && pName == "nu_mu" )      qB =  1;
517     // else if( !fProton && pName == "anti_nu_mu" ) qB = -1;
518 
519 
520       ClusterDecay( lvX, qB );
521   }
522   return &theParticleChange;
523 }
524 
525 
526 /////////////////////////////////////////////////////////////////////
527 ////////////////////////////////////////////////////////////////////
528 ///////////////////////////////////////////////////////////////////
529 
530 /////////////////////////////////////////////////
531 //
532 // sample x, then Q2
533 
534 void G4NuMuNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
535 {
536   fBreak = false;
537   G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100); 
538   G4int Z = targetNucleus.GetZ_asInt(); 
539   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
540   G4double Ex(0.), ei(0.), nm2(0.);
541   G4double cost(1.), sint(0.), phi(0.), muMom(0.); 
542   G4ThreeVector eP, bst;
543   const G4HadProjectile* aParticle = &aTrack;
544   G4LorentzVector lvp1 = aParticle->Get4Momentum();
545 
546   if( A == 1 ) // hydrogen, no Fermi motion ???
547   {
548     fNuEnergy = aParticle->GetTotalEnergy();
549     iTer = 0;
550 
551     do
552     {
553       fXsample = SampleXkr(fNuEnergy);
554       fQtransfer = SampleQkr(fNuEnergy, fXsample);
555       fQ2 = fQtransfer*fQtransfer;
556 
557      if( fXsample > 0. )
558       {
559         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
560         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
561       }
562       else
563       {
564         fW2 = fM1*fM1;
565         fEmu = fNuEnergy;
566       }
567       e3 = fNuEnergy + fM1 - fEmu;
568 
569       if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
570     
571       pMu2 = fEmu*fEmu - fMu*fMu;
572 
573       if(pMu2 < 0.) { fBreak = true; return; }
574 
575       pX2  = e3*e3 - fW2;
576 
577       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
578       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
579       iTer++;
580     }
581     while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax );
582 
583     if( iTer >= iTerMax ) { fBreak = true; return; }
584 
585     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
586     { 
587       G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
588       // fCosTheta = -1. + 2.*G4UniformRand(); 
589       if(fCosTheta < -1.) fCosTheta = -1.;
590       if(fCosTheta >  1.) fCosTheta =  1.;
591     }
592     // LVs
593 
594     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
595     G4LorentzVector lvsum = lvp1 + lvt1;
596 
597     cost = fCosTheta;
598     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
599     phi  = G4UniformRand()*CLHEP::twopi;
600     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
601     muMom = sqrt(fEmu*fEmu-fMu*fMu);
602     eP *= muMom;
603     fLVl = G4LorentzVector( eP, fEmu );
604 
605     fLVh = lvsum - fLVl;
606     fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
607   }
608   else // Fermi motion, Q2 in nucleon rest frame
609   {
610     G4Nucleus recoil1( A-1, Z );
611     rM = recoil1.AtomicMass(A-1,Z);   
612     do
613     {
614       // nMom = NucleonMomentumBR( targetNucleus ); // BR
615       nMom = GgSampleNM( targetNucleus ); // Gg
616       Ex = GetEx(A-1, fProton);
617       ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom );
618       //   ei = 0.5*( tM - s2M - 2*eX );
619     
620       nm2 = ei*ei - nMom*nMom;
621       iTer++;
622     }
623     while( nm2 < 0. && iTer < iTerMax ); 
624 
625     if( iTer >= iTerMax ) { fBreak = true; return; }
626     
627     G4ThreeVector nMomDir = nMom*G4RandomDirection();
628 
629     if( !f2p2h || A < 3 ) // 1p1h
630     {
631       // hM = tM - rM;
632 
633       fLVt = G4LorentzVector( -nMomDir, sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ) ); // rM ); //
634       fLVh = G4LorentzVector(  nMomDir, ei ); // hM); //
635     }
636     else // 2p2h
637     {
638       G4Nucleus recoil(A-2,Z-1);
639       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
640       hM = tM - rM;
641 
642       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
643       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom )  ); 
644     }
645     // G4cout<<hM<<", ";
646     // bst = fLVh.boostVector();
647 
648     // lvp1.boost(-bst); // -> nucleon rest system, where Q2 transfer is ???
649 
650     fNuEnergy  = lvp1.e();
651     // G4double mN = fLVh.m(); // better mN = fM1 !? vmg
652     iTer = 0;
653 
654     do // no FM!?, 5.4.20 vmg
655     {
656       fXsample = SampleXkr(fNuEnergy);
657       fQtransfer = SampleQkr(fNuEnergy, fXsample);
658       fQ2 = fQtransfer*fQtransfer;
659 
660       // G4double mR = mN + fM1*(A-1.)*std::exp(-2.0*fQtransfer/mN); // recoil mass in+el
661 
662       if( fXsample > 0. )
663       {
664         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
665 
666         // fW2 = mN*mN - fQ2 + fQ2/fXsample; // sample excited hadron mass
667         // fEmu = fNuEnergy - fQ2/2./mR/fXsample; // fM1->mN
668 
669         fEmu = fNuEnergy - fQ2/2./fM1/fXsample; // fM1->mN
670       }
671       else
672       {
673         // fW2 = mN*mN;
674 
675         fW2 = fM1*fM1; 
676         fEmu = fNuEnergy;
677       }
678       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
679       // e3 = fNuEnergy + mR - fEmu;
680 
681       e3 = fNuEnergy + fM1 - fEmu;
682 
683       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
684     
685       pMu2 = fEmu*fEmu - fMu*fMu;
686       pX2  = e3*e3 - fW2;
687 
688       if(pMu2 < 0.) { fBreak = true; return; }
689 
690       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
691       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
692       iTer++;
693     }
694     while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax );
695 
696     if( iTer >= iTerMax ) { fBreak = true; return; }
697 
698     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
699     { 
700       G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
701       // fCosTheta = -1. + 2.*G4UniformRand(); 
702       if( fCosTheta < -1.) fCosTheta = -1.;
703       if( fCosTheta >  1.) fCosTheta =  1.;
704     }
705     // LVs
706     // G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., mN ); // fM1 );
707 
708     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 ); // fM1 );
709     G4LorentzVector lvsum = lvp1 + lvt1;
710 
711     cost = fCosTheta;
712     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
713     phi  = G4UniformRand()*CLHEP::twopi;
714     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
715     muMom = sqrt(fEmu*fEmu-fMu*fMu);
716     eP *= muMom;
717     fLVl = G4LorentzVector( eP, fEmu );
718     fLVh = lvsum - fLVl;
719 
720     // if( fLVh.e() < mN || fLVh.m2() < 0.) { fBreak = true; return; }
721 
722     if( fLVh.e() < fM1 || fLVh.m2() < 0.) { fBreak = true; return; }
723 
724     // back to lab system
725 
726     // fLVl.boost(bst);
727     // fLVh.boost(bst);
728   }
729   //G4cout<<iTer<<", "<<fBreak<<"; ";
730 }
731 
732 //
733 //
734 ///////////////////////////
735