Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4ANuElNucleusCcModel.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: G4ANuElNucleusCcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
 27 //
 28 // Geant4 Header : G4ANuElNucleusCcModel
 29 //
 30 // Author : V.Grichine 12.2.19
 31 //  
 32 
 33 #include <iostream>
 34 #include <fstream>
 35 #include <sstream>
 36 
 37 #include "G4ANuElNucleusCcModel.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 "G4Positron.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 G4ANuElNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
 91 #endif     
 92 
 93 
 94 G4ANuElNucleusCcModel::G4ANuElNucleusCcModel(const G4String& name) 
 95   : G4NeutrinoNucleusModel(name)
 96 {
 97   thePositron = G4Positron::Positron();
 98   fData = fMaster = false;
 99   fMel = electron_mass_c2;
100   InitialiseModel();  
101 }
102 
103 
104 G4ANuElNucleusCcModel::~G4ANuElNucleusCcModel()
105 {}
106 
107 
108 void G4ANuElNucleusCcModel::ModelDescription(std::ostream& outFile) const
109 {
110 
111     outFile << "G4ANuElNucleusCcModel 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 G4ANuElNucleusCcModel::InitialiseModel()
122 {
123   G4String pName  = "anti_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 G4ANuElNucleusCcModel::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 == "anti_nu_e"    
227         &&
228         energy > fMinNuEnergy                                )
229   {
230     result = true;
231   }
232 
233   return result;
234 }
235 
236 /////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
237 //
238 //
239 
240 G4HadFinalState* G4ANuElNucleusCcModel::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 == "anti_nu_e" )         aLept = new G4DynamicParticle( thePositron, lv2 );  
327     else
328     {
329       theParticleChange.SetEnergyChange(energy);
330       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
331       return &theParticleChange;
332     }
333     if( pName == "anti_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 == "anti_nu_e" )         aLept = new G4DynamicParticle( thePositron, 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 == "anti_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 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     if( pName == "anti_nu_e" ) // (++) state -> p + pi+
436     { 
437       fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
438           + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
439     }
440     else // (0) state -> p + pi-, n + pi0
441     {
442       // fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
443       //     + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass();
444     } 
445   }
446   else // excited neutron
447   {
448     fProton = false;
449     recoil = G4Nucleus(A-1,Z);
450     fRecoil = &recoil;
451     if( pName == "anti_nu_e" ) // (+) state -> n + pi+
452     {      
453       fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
454           + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
455     }
456     else // (-) state -> n + pi-, // n + pi0
457     {
458       // fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
459       //     + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass();
460     } 
461   }
462   // G4int       index = GetEnergyIndex(energy);
463   G4int nepdg = aParticle->GetDefinition()->GetPDGEncoding();
464 
465   G4double qeTotRat; //  = GetNuMuQeTotRat(index, energy);
466   qeTotRat = CalculateQEratioA( Z, A, energy, nepdg);
467 
468   G4ThreeVector dX = (lvX.vect()).unit();
469   G4double eX   = lvX.e();  // excited nucleon
470   G4double mX   = sqrt(massX2);
471   // G4double pX   = sqrt( eX*eX - mX*mX );
472   // G4double sumE = eX + rM;
473 
474   if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE
475   {  
476     fString = false;
477 
478     G4double rM;
479     if( fProton ) 
480     {  
481       fPDGencoding = 2212;
482       fMr =  proton_mass_c2;
483       recoil = G4Nucleus(A-1,Z-1);
484       fRecoil = &recoil;
485       rM = recoil.AtomicMass(A-1,Z-1);
486     } 
487     else 
488     {  
489       fPDGencoding = 2112;
490       fMr =   G4ParticleTable::GetParticleTable()->
491   FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV;
492       recoil = G4Nucleus(A-1,Z);
493       fRecoil = &recoil;
494       rM = recoil.AtomicMass(A-1,Z);
495     } 
496     // sumE = eX + rM;   
497     G4double eTh = fMr + 0.5*(fMr*fMr - mX*mX)/rM;
498 
499     if( eX <= eTh ) // vmg, very rarely out of kinematics
500     {
501       fString = true;
502       theParticleChange.SetEnergyChange(energy);
503       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
504       return &theParticleChange;
505     }
506     // FinalBarion( fLVh, 0, fPDGencoding ); // p(n)+deexcited recoil
507     FinalBarion( lvX, 0, fPDGencoding ); // p(n)+deexcited recoil
508   }
509   else // if ( eX < 9500000.*GeV ) // <  25.*GeV) // < 95.*GeV ) // < 2.5*GeV ) //cluster decay
510   {  
511     if     (  fProton && pName == "anti_nu_e" )      qB =  2;
512     else if( !fProton && pName == "anti_nu_e" )      qB =  1;
513 
514     ClusterDecay( lvX, qB );
515   }
516   return &theParticleChange;
517 }
518 
519 
520 /////////////////////////////////////////////////////////////////////
521 ////////////////////////////////////////////////////////////////////
522 ///////////////////////////////////////////////////////////////////
523 
524 /////////////////////////////////////////////////
525 //
526 // sample x, then Q2
527 
528 void G4ANuElNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
529 {
530   fBreak = false;
531   G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100); 
532   G4int Z = targetNucleus.GetZ_asInt(); 
533   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
534   G4double Ex(0.), ei(0.), nm2(0.);
535   G4double cost(1.), sint(0.), phi(0.), muMom(0.); 
536   G4ThreeVector eP, bst;
537   const G4HadProjectile* aParticle = &aTrack;
538   G4LorentzVector lvp1 = aParticle->Get4Momentum();
539 
540   if( A == 1 ) // hydrogen, no Fermi motion ???
541   {
542     fNuEnergy = aParticle->GetTotalEnergy();
543     iTer = 0;
544 
545     do
546     {
547       fXsample = SampleXkr(fNuEnergy);
548       fQtransfer = SampleQkr(fNuEnergy, fXsample);
549       fQ2 = fQtransfer*fQtransfer;
550 
551      if( fXsample > 0. )
552       {
553         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
554         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
555       }
556       else
557       {
558         fW2 = fM1*fM1;
559         fEmu = fNuEnergy;
560       }
561       e3 = fNuEnergy + fM1 - fEmu;
562 
563       if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
564     
565       pMu2 = fEmu*fEmu - fMel*fMel;
566 
567       if(pMu2 < 0.) { fBreak = true; return; }
568 
569       pX2  = e3*e3 - fW2;
570 
571       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
572       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
573       iTer++;
574     }
575     while( ( abs(fCosTheta) > 1. || fEmu < fMel ) && iTer < iTerMax );
576 
577     if( iTer >= iTerMax ) { fBreak = true; return; }
578 
579     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
580     { 
581       G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
582       // fCosTheta = -1. + 2.*G4UniformRand(); 
583       if(fCosTheta < -1.) fCosTheta = -1.;
584       if(fCosTheta >  1.) fCosTheta =  1.;
585     }
586     // LVs
587 
588     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
589     G4LorentzVector lvsum = lvp1 + lvt1;
590 
591     cost = fCosTheta;
592     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
593     phi  = G4UniformRand()*CLHEP::twopi;
594     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
595     muMom = sqrt(fEmu*fEmu-fMel*fMel);
596     eP *= muMom;
597     fLVl = G4LorentzVector( eP, fEmu );
598 
599     fLVh = lvsum - fLVl;
600     fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
601   }
602   else // Fermi motion, Q2 in nucleon rest frame
603   {
604     G4Nucleus recoil1( A-1, Z );
605     rM = recoil1.AtomicMass(A-1,Z);   
606     do
607     {
608       // nMom = NucleonMomentumBR( targetNucleus ); // BR
609       nMom = GgSampleNM( targetNucleus ); // Gg
610       Ex = GetEx(A-1, fProton);
611       ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom );
612       //   ei = 0.5*( tM - s2M - 2*eX );
613     
614       nm2 = ei*ei - nMom*nMom;
615       iTer++;
616     }
617     while( nm2 < 0. && iTer < iTerMax ); 
618 
619     if( iTer >= iTerMax ) { fBreak = true; return; }
620     
621     G4ThreeVector nMomDir = nMom*G4RandomDirection();
622 
623     if( !f2p2h || A < 3 ) // 1p1h
624     {
625       // hM = tM - rM;
626 
627       fLVt = G4LorentzVector( -nMomDir, sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ) ); // rM ); //
628       fLVh = G4LorentzVector(  nMomDir, ei ); // hM); //
629     }
630     else // 2p2h
631     {
632       G4Nucleus recoil(A-2,Z-1);
633       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
634       hM = tM - rM;
635 
636       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
637       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom )  ); 
638     }
639     // G4cout<<hM<<", ";
640     // bst = fLVh.boostVector();
641 
642     // lvp1.boost(-bst); // -> nucleon rest system, where Q2 transfer is ???
643 
644     fNuEnergy  = lvp1.e();
645     // G4double mN = fLVh.m(); // better mN = fM1 !? vmg
646     iTer = 0;
647 
648     do // no FM!?, 5.4.20 vmg
649     {
650       fXsample = SampleXkr(fNuEnergy);
651       fQtransfer = SampleQkr(fNuEnergy, fXsample);
652       fQ2 = fQtransfer*fQtransfer;
653 
654       // G4double mR = mN + fM1*(A-1.)*std::exp(-2.0*fQtransfer/mN); // recoil mass in+el
655 
656       if( fXsample > 0. )
657       {
658         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
659 
660         // fW2 = mN*mN - fQ2 + fQ2/fXsample; // sample excited hadron mass
661         // fEmu = fNuEnergy - fQ2/2./mR/fXsample; // fM1->mN
662 
663         fEmu = fNuEnergy - fQ2/2./fM1/fXsample; // fM1->mN
664       }
665       else
666       {
667         // fW2 = mN*mN;
668 
669         fW2 = fM1*fM1; 
670         fEmu = fNuEnergy;
671       }
672       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
673       // e3 = fNuEnergy + mR - fEmu;
674 
675       e3 = fNuEnergy + fM1 - fEmu;
676 
677       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
678     
679       pMu2 = fEmu*fEmu - fMel*fMel;
680       pX2  = e3*e3 - fW2;
681 
682       if(pMu2 < 0.) { fBreak = true; return; }
683 
684       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
685       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
686       iTer++;
687     }
688     while( ( abs(fCosTheta) > 1. || fEmu < fMel ) && iTer < iTerMax );
689 
690     if( iTer >= iTerMax ) { fBreak = true; return; }
691 
692     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
693     { 
694       G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
695       // fCosTheta = -1. + 2.*G4UniformRand(); 
696       if( fCosTheta < -1.) fCosTheta = -1.;
697       if( fCosTheta >  1.) fCosTheta =  1.;
698     }
699     // LVs
700     // G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., mN ); // fM1 );
701 
702     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 ); // fM1 );
703     G4LorentzVector lvsum = lvp1 + lvt1;
704 
705     cost = fCosTheta;
706     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
707     phi  = G4UniformRand()*CLHEP::twopi;
708     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
709     muMom = sqrt(fEmu*fEmu-fMel*fMel);
710     eP *= muMom;
711     fLVl = G4LorentzVector( eP, fEmu );
712     fLVh = lvsum - fLVl;
713 
714     // if( fLVh.e() < mN || fLVh.m2() < 0.) { fBreak = true; return; }
715 
716     if( fLVh.e() < fM1 || fLVh.m2() < 0.) { fBreak = true; return; }
717 
718     // back to lab system
719 
720     // fLVl.boost(bst);
721     // fLVh.boost(bst);
722   }
723   //G4cout<<iTer<<", "<<fBreak<<"; ";
724 }
725 
726 //
727 //
728 ///////////////////////////
729