Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4NuTauNucleusCcModel.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 ]

Diff markup

Differences between /processes/hadronic/models/lepto_nuclear/src/G4NuTauNucleusCcModel.cc (Version 11.3.0) and /processes/hadronic/models/lepto_nuclear/src/G4NuTauNucleusCcModel.cc (Version 10.4)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // $Id: G4NuTauNucleusCcModel.cc 91806 2015-08    
 27 //                                                
 28 // Geant4 Header : G4NuTauNucleusCcModel          
 29 //                                                
 30 // Author : V.Grichine 12.2.19                    
 31 //                                                
 32                                                   
 33 #include <iostream>                               
 34 #include <fstream>                                
 35 #include <sstream>                                
 36                                                   
 37 #include "G4NuTauNucleusCcModel.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 "G4TauMinus.hh"                          
 82 #include "G4TauPlus.hh"                           
 83 #include "G4Nucleus.hh"                           
 84 #include "G4LorentzVector.hh"                     
 85                                                   
 86 using namespace std;                              
 87 using namespace CLHEP;                            
 88                                                   
 89 #ifdef G4MULTITHREADED                            
 90     G4Mutex G4NuTauNucleusCcModel::numuNucleus    
 91 #endif                                            
 92                                                   
 93                                                   
 94 G4NuTauNucleusCcModel::G4NuTauNucleusCcModel(c    
 95   : G4NeutrinoNucleusModel(name)                  
 96 {                                                 
 97   fData = fMaster = false;                        
 98   fMtau = 1.77686*GeV;                            
 99   theTauMinus = G4TauMinus::TauMinus();           
100   theTauPlus = G4TauPlus::TauPlus();              
101   InitialiseModel();                              
102 }                                                 
103                                                   
104                                                   
105 G4NuTauNucleusCcModel::~G4NuTauNucleusCcModel(    
106 {}                                                
107                                                   
108                                                   
109 void G4NuTauNucleusCcModel::ModelDescription(s    
110 {                                                 
111                                                   
112     outFile << "G4NuTauNucleusCcModel is a tau    
113             << "model which uses the standard     
114             << "transfer parameterization.  Th    
115                                                   
116 }                                                 
117                                                   
118 //////////////////////////////////////////////    
119 //                                                
120 // Read data from G4PARTICLEXSDATA (locally PA    
121                                                   
122 void G4NuTauNucleusCcModel::InitialiseModel()     
123 {                                                 
124   G4String pName  = "nu_mu";                      
125   // G4String pName  = "anti_nu_mu";              
126                                                   
127   G4int nSize(0), i(0), j(0), k(0);               
128                                                   
129   if(!fData)                                      
130   {                                               
131 #ifdef G4MULTITHREADED                            
132     G4MUTEXLOCK(&numuNucleusModel);               
133     if(!fData)                                    
134     {                                             
135 #endif                                            
136       fMaster = true;                             
137 #ifdef G4MULTITHREADED                            
138     }                                             
139     G4MUTEXUNLOCK(&numuNucleusModel);             
140 #endif                                            
141   }                                               
142                                                   
143   if(fMaster)                                     
144   {                                               
145     const char* path = G4FindDataDir("G4PARTIC    
146     std::ostringstream ost1, ost2, ost3, ost4;    
147     ost1 << path << "/" << "neutrino" << "/" <    
148                                                   
149     std::ifstream filein1( ost1.str().c_str()     
150                                                   
151     // filein.open("$PARTICLEXSDATA/");           
152                                                   
153     filein1>>nSize;                               
154                                                   
155     for( k = 0; k < fNbin; ++k )                  
156     {                                             
157       for( i = 0; i <= fNbin; ++i )               
158       {                                           
159         filein1 >> fNuMuXarrayKR[k][i];           
160         // G4cout<< fNuMuXarrayKR[k][i] << "      
161       }                                           
162     }                                             
163     // G4cout<<G4endl<<G4endl;                    
164                                                   
165     ost2 << path << "/" << "neutrino" << "/" <    
166     std::ifstream  filein2( ost2.str().c_str()    
167                                                   
168     filein2>>nSize;                               
169                                                   
170     for( k = 0; k < fNbin; ++k )                  
171     {                                             
172       for( i = 0; i < fNbin; ++i )                
173       {                                           
174         filein2 >> fNuMuXdistrKR[k][i];           
175         // G4cout<< fNuMuXdistrKR[k][i] << "      
176       }                                           
177     }                                             
178     // G4cout<<G4endl<<G4endl;                    
179                                                   
180     ost3 << path << "/" << "neutrino" << "/" <    
181     std::ifstream  filein3( ost3.str().c_str()    
182                                                   
183     filein3>>nSize;                               
184                                                   
185     for( k = 0; k < fNbin; ++k )                  
186     {                                             
187       for( i = 0; i <= fNbin; ++i )               
188       {                                           
189         for( j = 0; j <= fNbin; ++j )             
190         {                                         
191           filein3 >> fNuMuQarrayKR[k][i][j];      
192           // G4cout<< fNuMuQarrayKR[k][i][j] <    
193         }                                         
194       }                                           
195     }                                             
196     // G4cout<<G4endl<<G4endl;                    
197                                                   
198     ost4 << path << "/" << "neutrino" << "/" <    
199     std::ifstream  filein4( ost4.str().c_str()    
200                                                   
201     filein4>>nSize;                               
202                                                   
203     for( k = 0; k < fNbin; ++k )                  
204     {                                             
205       for( i = 0; i <= fNbin; ++i )               
206       {                                           
207         for( j = 0; j < fNbin; ++j )              
208         {                                         
209           filein4 >> fNuMuQdistrKR[k][i][j];      
210           // G4cout<< fNuMuQdistrKR[k][i][j] <    
211         }                                         
212       }                                           
213     }                                             
214     fData = true;                                 
215   }                                               
216 }                                                 
217                                                   
218 //////////////////////////////////////////////    
219                                                   
220 G4bool G4NuTauNucleusCcModel::IsApplicable(con    
221                   G4Nucleus & )                   
222 {                                                 
223   G4bool result  = false;                         
224   G4String pName = aPart.GetDefinition()->GetP    
225   G4double energy = aPart.GetTotalEnergy();       
226                                                   
227   if(  pName == "nu_tau"  // || pName == "anti    
228         &&                                        
229         energy > fMinNuEnergy                     
230   {                                               
231     result = true;                                
232   }                                               
233                                                   
234   return result;                                  
235 }                                                 
236                                                   
237 /////////////////////////////////////////// Cl    
238 //                                                
239 //                                                
240                                                   
241 G4HadFinalState* G4NuTauNucleusCcModel::ApplyY    
242      const G4HadProjectile& aTrack, G4Nucleus&    
243 {                                                 
244   theParticleChange.Clear();                      
245   fProton = f2p2h = fBreak = false;               
246   fCascade = fString  = false;                    
247   fLVh = fLVl = fLVt = fLVcpi = G4LorentzVecto    
248                                                   
249   const G4HadProjectile* aParticle = &aTrack;     
250   G4double energy = aParticle->GetTotalEnergy(    
251                                                   
252   G4String pName  = aParticle->GetDefinition()    
253                                                   
254   if( energy < fMinNuEnergy )                     
255   {                                               
256     theParticleChange.SetEnergyChange(energy);    
257     theParticleChange.SetMomentumChange(aTrack    
258     return &theParticleChange;                    
259   }                                               
260                                                   
261   SampleLVkr( aTrack, targetNucleus);             
262                                                   
263   if( fBreak == true || fEmu < fMtau ) // ~5*1    
264   {                                               
265     // G4cout<<"ni, ";                            
266     theParticleChange.SetEnergyChange(energy);    
267     theParticleChange.SetMomentumChange(aTrack    
268     return &theParticleChange;                    
269   }                                               
270                                                   
271   // LVs of initial state                         
272                                                   
273   G4LorentzVector lvp1 = aParticle->Get4Moment    
274   G4LorentzVector lvt1( 0., 0., 0., fM1 );        
275   G4double mPip = G4ParticleTable::GetParticle    
276                                                   
277   // 1-pi by fQtransfer && nu-energy              
278   G4LorentzVector lvpip1( 0., 0., 0., mPip );     
279   G4LorentzVector lvsum, lv2, lvX;                
280   G4ThreeVector eP;                               
281   G4double cost(1.), sint(0.), phi(0.), muMom(    
282   G4DynamicParticle* aLept = nullptr; // lepto    
283                                                   
284   G4int Z = targetNucleus.GetZ_asInt();           
285   G4int A = targetNucleus.GetA_asInt();           
286   G4double  mTarg = targetNucleus.AtomicMass(A    
287   G4int pdgP(0), qB(0);                           
288   // G4double mSum = G4ParticleTable::GetParti    
289                                                   
290   G4int iPi     = GetOnePionIndex(energy);        
291   G4double p1pi = GetNuMuOnePionProb( iPi, ene    
292                                                   
293   if( p1pi > G4UniformRand()  && fCosTheta > 0    
294   {                                               
295     // lvsum = lvp1 + lvpip1;                     
296     lvsum = lvp1 + lvt1;                          
297     // cost = fCosThetaPi;                        
298     cost = fCosTheta;                             
299     sint = std::sqrt( (1.0 - cost)*(1.0 + cost    
300     phi  = G4UniformRand()*CLHEP::twopi;          
301     eP   = G4ThreeVector( sint*std::cos(phi),     
302                                                   
303     // muMom = sqrt(fEmuPi*fEmuPi-fMtau*fMtau)    
304     muMom = sqrt(fEmu*fEmu-fMtau*fMtau);          
305                                                   
306     eP *= muMom;                                  
307                                                   
308     // lv2 = G4LorentzVector( eP, fEmuPi );       
309     // lv2 = G4LorentzVector( eP, fEmu );         
310     lv2 = fLVl;                                   
311                                                   
312     // lvX = lvsum - lv2;                         
313     lvX = fLVh;                                   
314     massX2 = lvX.m2();                            
315     massX = lvX.m();                              
316     massR = fLVt.m();                             
317                                                   
318     if ( massX2 <= 0. ) // vmg: very rarely ~     
319     {                                             
320       fCascade = true;                            
321       theParticleChange.SetEnergyChange(energy    
322       theParticleChange.SetMomentumChange(aTra    
323       return &theParticleChange;                  
324     }                                             
325     fW2 = massX2;                                 
326                                                   
327     if(  pName == "nu_tau" )         aLept = n    
328     // else if( pName == "anti_nu_tau") aLept     
329     else                                          
330     {                                             
331       theParticleChange.SetEnergyChange(energy    
332       theParticleChange.SetMomentumChange(aTra    
333       return &theParticleChange;                  
334     }                                             
335     if( pName == "nu_tau" ) pdgP =  211;          
336     // else                   pdgP = -211;        
337     // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mT    
338                                                   
339     if( A > 1 )                                   
340     {                                             
341       eCut = (fMpi + mTarg)*(fMpi + mTarg) - (    
342       eCut /= 2.*massR;                           
343       eCut += massX;                              
344     }                                             
345     else  eCut = fM1 + fMpi;                      
346                                                   
347     if ( lvX.e() > eCut ) // && sqrt( GetW2()     
348     {                                             
349       CoherentPion( lvX, pdgP, targetNucleus);    
350     }                                             
351     else                                          
352     {                                             
353       fCascade = true;                            
354       theParticleChange.SetEnergyChange(energy    
355       theParticleChange.SetMomentumChange(aTra    
356       return &theParticleChange;                  
357     }                                             
358     theParticleChange.AddSecondary( aLept, fSe    
359                                                   
360     return &theParticleChange;                    
361   }                                               
362   else // lepton part in lab                      
363   {                                               
364     lvsum = lvp1 + lvt1;                          
365     cost = fCosTheta;                             
366     sint = std::sqrt( (1.0 - cost)*(1.0 + cost    
367     phi  = G4UniformRand()*CLHEP::twopi;          
368     eP   = G4ThreeVector( sint*std::cos(phi),     
369                                                   
370     muMom = sqrt(fEmu*fEmu-fMtau*fMtau);          
371                                                   
372     eP *= muMom;                                  
373                                                   
374     lv2 = G4LorentzVector( eP, fEmu );            
375     lv2 = fLVl;                                   
376     lvX = lvsum - lv2;                            
377     lvX = fLVh;                                   
378     massX2 = lvX.m2();                            
379                                                   
380     if ( massX2 <= 0. ) // vmg: very rarely ~     
381     {                                             
382       fCascade = true;                            
383       theParticleChange.SetEnergyChange(energy    
384       theParticleChange.SetMomentumChange(aTra    
385       return &theParticleChange;                  
386     }                                             
387     fW2 = massX2;                                 
388                                                   
389     if(  pName == "nu_tau" )         aLept = n    
390     // else if( pName == "anti_nu_tau") aLept     
391     else                                          
392     {                                             
393       theParticleChange.SetEnergyChange(energy    
394       theParticleChange.SetMomentumChange(aTra    
395       return &theParticleChange;                  
396     }                                             
397     theParticleChange.AddSecondary( aLept, fSe    
398   }                                               
399                                                   
400   // hadron part                                  
401                                                   
402   fRecoil  = nullptr;                             
403                                                   
404   if( A == 1 )                                    
405   {                                               
406     if( pName == "nu_tau" ) qB = 2;               
407     // else                   qB = 0;             
408                                                   
409     // if( G4UniformRand() > 0.1 ) //  > 0.999    
410     {                                             
411       ClusterDecay( lvX, qB );                    
412     }                                             
413     return &theParticleChange;                    
414   }                                               
415     /*                                            
416     // else                                       
417     {                                             
418       if( pName == "nu_tau" ) pdgP =  211;        
419       else                   pdgP = -211;         
420                                                   
421                                                   
422       if ( fQtransfer < 0.95*GeV ) // < 0.35*G    
423       {                                           
424   if( lvX.m() > mSum ) CoherentPion( lvX, pdgP    
425       }                                           
426     }                                             
427     return &theParticleChange;                    
428   }                                               
429   */                                              
430   G4Nucleus recoil;                               
431   G4double ratio = G4double(Z)/G4double(A);       
432                                                   
433   if( ratio > G4UniformRand() ) // proton is e    
434   {                                               
435     fProton = true;                               
436     recoil = G4Nucleus(A-1,Z-1);                  
437     fRecoil = &recoil;                            
438     if( pName == "nu_tau" ) // (++) state -> p    
439     {                                             
440       fMt = G4ParticleTable::GetParticleTable(    
441           + G4ParticleTable::GetParticleTable(    
442     }                                             
443     else // (0) state -> p + pi-, n + pi0         
444     {                                             
445       // fMt = G4ParticleTable::GetParticleTab    
446       //     + G4ParticleTable::GetParticleTab    
447     }                                             
448   }                                               
449   else // excited neutron                         
450   {                                               
451     fProton = false;                              
452     recoil = G4Nucleus(A-1,Z);                    
453     fRecoil = &recoil;                            
454     if( pName == "nu_tau" ) // (+) state -> n     
455     {                                             
456       fMt = G4ParticleTable::GetParticleTable(    
457           + G4ParticleTable::GetParticleTable(    
458     }                                             
459     else // (-) state -> n + pi-, // n + pi0      
460     {                                             
461       // fMt = G4ParticleTable::GetParticleTab    
462       //     + G4ParticleTable::GetParticleTab    
463     }                                             
464   }                                               
465   // G4int       index = GetEnergyIndex(energy    
466   G4int nepdg = aParticle->GetDefinition()->Ge    
467                                                   
468   G4double qeTotRat; // = GetNuMuQeTotRat(inde    
469   qeTotRat = CalculateQEratioA( Z, A, energy,     
470                                                   
471   G4ThreeVector dX = (lvX.vect()).unit();         
472   G4double eX   = lvX.e();  // excited nucleon    
473   G4double mX   = sqrt(massX2);                   
474   // G4double pX   = sqrt( eX*eX - mX*mX );       
475   // G4double sumE = eX + rM;                     
476                                                   
477   if( qeTotRat > G4UniformRand() || mX <= fMt     
478   {                                               
479     fString = false;                              
480                                                   
481     G4double rM;                                  
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_tau" )          
491     {                                             
492       fPDGencoding = 2112;                        
493       fMr =   G4ParticleTable::GetParticleTabl    
494   FindParticle(fPDGencoding)->GetPDGMass(); //    
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)    
501                                                   
502     if( eX <= eTh ) // vmg, very rarely out of    
503     {                                             
504       fString = true;                             
505       theParticleChange.SetEnergyChange(energy    
506       theParticleChange.SetMomentumChange(aTra    
507       return &theParticleChange;                  
508     }                                             
509     // FinalBarion( fLVh, 0, fPDGencoding ); /    
510     FinalBarion( lvX, 0, fPDGencoding ); // p(    
511   }                                               
512   else // if ( eX < 9500000.*GeV ) // <  25.*G    
513   {                                               
514     if     (  fProton && pName == "nu_tau" )      
515     // else if(  fProton && pName == "anti_nu_    
516     else if( !fProton && pName == "nu_tau" )      
517     // else if( !fProton && pName == "anti_nu_    
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 G4NuTauNucleusCcModel::SampleLVkr(const G    
535 {                                                 
536   fBreak = false;                                 
537   G4int A = targetNucleus.GetA_asInt(), iTer(0    
538   G4int Z = targetNucleus.GetZ_asInt();           
539   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.)    
540   G4double Ex(0.), ei(0.), nm2(0.);               
541   G4double cost(1.), sint(0.), phi(0.), muMom(    
542   G4ThreeVector eP, bst;                          
543   const G4HadProjectile* aParticle = &aTrack;     
544   G4LorentzVector lvp1 = aParticle->Get4Moment    
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, fXsamp    
555       fQ2 = fQtransfer*fQtransfer;                
556                                                   
557      if( fXsample > 0. )                          
558       {                                           
559         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; //    
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 =    
570                                                   
571       pMu2 = fEmu*fEmu - fMtau*fMtau;             
572                                                   
573       if(pMu2 < 0.) { fBreak = true; return; }    
574                                                   
575       pX2  = e3*e3 - fW2;                         
576                                                   
577       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2    
578       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);       
579       iTer++;                                     
580     }                                             
581     while( ( abs(fCosTheta) > 1. || fEmu < fMt    
582                                                   
583     if( iTer >= iTerMax ) { fBreak = true; ret    
584                                                   
585     if( abs(fCosTheta) > 1.) // vmg: due to bi    
586     {                                             
587       G4cout<<"H2: fCosTheta = "<<fCosTheta<<"    
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    
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),     
601     muMom = sqrt(fEmu*fEmu-fMtau*fMtau);          
602     eP *= muMom;                                  
603     fLVl = G4LorentzVector( eP, fEmu );           
604                                                   
605     fLVh = lvsum - fLVl;                          
606     fLVt = G4LorentzVector( 0., 0., 0., 0. );     
607   }                                               
608   else // Fermi motion, Q2 in nucleon rest fra    
609   {                                               
610     G4Nucleus recoil1( A-1, Z );                  
611     rM = recoil1.AtomicMass(A-1,Z);               
612     do                                            
613     {                                             
614       // nMom = NucleonMomentumBR( targetNucle    
615       nMom = GgSampleNM( targetNucleus ); // G    
616       Ex = GetEx(A-1, fProton);                   
617       ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nM    
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; ret    
626                                                   
627     G4ThreeVector nMomDir = nMom*G4RandomDirec    
628                                                   
629     if( !f2p2h || A < 3 ) // 1p1h                 
630     {                                             
631       // hM = tM - rM;                            
632                                                   
633       fLVt = G4LorentzVector( -nMomDir, sqrt(     
634       fLVh = G4LorentzVector(  nMomDir, ei );     
635     }                                             
636     else // 2p2h                                  
637     {                                             
638       G4Nucleus recoil(A-2,Z-1);                  
639       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMo    
640       hM = tM - rM;                               
641                                                   
642       fLVt = G4LorentzVector( nMomDir, sqrt( r    
643       fLVh = G4LorentzVector(-nMomDir, sqrt( h    
644     }                                             
645     // G4cout<<hM<<", ";                          
646     // bst = fLVh.boostVector();                  
647                                                   
648     // lvp1.boost(-bst); // -> nucleon rest sy    
649                                                   
650     fNuEnergy  = lvp1.e();                        
651     // G4double mN = fLVh.m(); // better mN =     
652     iTer = 0;                                     
653                                                   
654     do // no FM!?, 5.4.20 vmg                     
655     {                                             
656       fXsample = SampleXkr(fNuEnergy);            
657       fQtransfer = SampleQkr(fNuEnergy, fXsamp    
658       fQ2 = fQtransfer*fQtransfer;                
659                                                   
660       // G4double mR = mN + fM1*(A-1.)*std::ex    
661                                                   
662       if( fXsample > 0. )                         
663       {                                           
664         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; //    
665                                                   
666         // fW2 = mN*mN - fQ2 + fQ2/fXsample; /    
667         // fEmu = fNuEnergy - fQ2/2./mR/fXsamp    
668                                                   
669         fEmu = fNuEnergy - fQ2/2./fM1/fXsample    
670       }                                           
671       else                                        
672       {                                           
673         // fW2 = mN*mN;                           
674                                                   
675         fW2 = fM1*fM1;                            
676         fEmu = fNuEnergy;                         
677       }                                           
678       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu    
679       // e3 = fNuEnergy + mR - fEmu;              
680                                                   
681       e3 = fNuEnergy + fM1 - fEmu;                
682                                                   
683       // if( e3 < sqrt(fW2) )  G4cout<<"energy    
684                                                   
685       pMu2 = fEmu*fEmu - fMtau*fMtau;             
686       pX2  = e3*e3 - fW2;                         
687                                                   
688       if(pMu2 < 0.) { fBreak = true; return; }    
689                                                   
690       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2    
691       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);       
692       iTer++;                                     
693     }                                             
694     while( ( abs(fCosTheta) > 1. || fEmu < fMt    
695                                                   
696     if( iTer >= iTerMax ) { fBreak = true; ret    
697                                                   
698     if( abs(fCosTheta) > 1.) // vmg: due to bi    
699     {                                             
700       G4cout<<"FM: fCosTheta = "<<fCosTheta<<"    
701       // fCosTheta = -1. + 2.*G4UniformRand();    
702       if( fCosTheta < -1.) fCosTheta = -1.;       
703       if( fCosTheta >  1.) fCosTheta =  1.;       
704     }                                             
705     // LVs                                        
706     // G4LorentzVector lvt1  = G4LorentzVector    
707                                                   
708     G4LorentzVector lvt1  = G4LorentzVector( 0    
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),     
715     muMom = sqrt(fEmu*fEmu-fMtau*fMtau);          
716     eP *= muMom;                                  
717     fLVl = G4LorentzVector( eP, fEmu );           
718     fLVh = lvsum - fLVl;                          
719                                                   
720     // if( fLVh.e() < mN || fLVh.m2() < 0.) {     
721                                                   
722     if( fLVh.e() < fM1 || fLVh.m2() < 0.) { fB    
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