Geant4 Cross Reference

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