Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lepto_nuclear/src/G4NuMuNucleusNcModel.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/G4NuMuNucleusNcModel.cc (Version 11.3.0) and /processes/hadronic/models/lepto_nuclear/src/G4NuMuNucleusNcModel.cc (Version 11.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // $Id: G4NuMuNucleusNcModel.cc 91806 2015-08-     26 // $Id: G4NuMuNucleusNcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $
 27 //                                                 27 //
 28 // Geant4 Header : G4NuMuNucleusNcModel            28 // Geant4 Header : G4NuMuNucleusNcModel
 29 //                                                 29 //
 30 // Author : V.Grichine 12.2.19                     30 // Author : V.Grichine 12.2.19
 31 //                                                 31 //  
 32                                                    32 
 33 #include "G4NuMuNucleusNcModel.hh"                 33 #include "G4NuMuNucleusNcModel.hh"
 34 #include "G4NeutrinoNucleusModel.hh"               34 #include "G4NeutrinoNucleusModel.hh" 
 35                                                    35 
 36 // #include "G4NuMuResQX.hh"                       36 // #include "G4NuMuResQX.hh" 
 37                                                    37 
 38 #include "G4SystemOfUnits.hh"                      38 #include "G4SystemOfUnits.hh"
 39 #include "G4ParticleTable.hh"                      39 #include "G4ParticleTable.hh"
 40 #include "G4ParticleDefinition.hh"                 40 #include "G4ParticleDefinition.hh"
 41 #include "G4IonTable.hh"                           41 #include "G4IonTable.hh"
 42 #include "Randomize.hh"                            42 #include "Randomize.hh"
 43 #include "G4RandomDirection.hh"                    43 #include "G4RandomDirection.hh"
 44                                                    44 
 45 // #include "G4Integrator.hh"                      45 // #include "G4Integrator.hh"
 46 #include "G4DataVector.hh"                         46 #include "G4DataVector.hh"
 47 #include "G4PhysicsTable.hh"                       47 #include "G4PhysicsTable.hh"
 48 #include "G4KineticTrack.hh"                       48 #include "G4KineticTrack.hh"
 49 #include "G4DecayKineticTracks.hh"                 49 #include "G4DecayKineticTracks.hh"
 50 #include "G4KineticTrackVector.hh"                 50 #include "G4KineticTrackVector.hh"
 51 #include "G4Fragment.hh"                           51 #include "G4Fragment.hh"
 52 #include "G4ReactionProductVector.hh"              52 #include "G4ReactionProductVector.hh"
 53                                                    53 
 54                                                    54 
 55 #include "G4NeutrinoMu.hh"                         55 #include "G4NeutrinoMu.hh"
 56 #include "G4AntiNeutrinoMu.hh"                     56 #include "G4AntiNeutrinoMu.hh"
 57 #include "G4Nucleus.hh"                            57 #include "G4Nucleus.hh"
 58 #include "G4LorentzVector.hh"                      58 #include "G4LorentzVector.hh"
 59                                                    59 
 60 using namespace std;                               60 using namespace std;
 61 using namespace CLHEP;                             61 using namespace CLHEP;
 62                                                    62 
 63 #ifdef G4MULTITHREADED                             63 #ifdef G4MULTITHREADED
 64     G4Mutex G4NuMuNucleusNcModel::numuNucleusM     64     G4Mutex G4NuMuNucleusNcModel::numuNucleusModel = G4MUTEX_INITIALIZER;
 65 #endif                                             65 #endif     
 66                                                    66 
 67                                                    67 
 68 G4NuMuNucleusNcModel::G4NuMuNucleusNcModel(con     68 G4NuMuNucleusNcModel::G4NuMuNucleusNcModel(const G4String& name) 
 69   : G4NeutrinoNucleusModel(name)                   69   : G4NeutrinoNucleusModel(name)
 70 {                                                  70 {
 71   SetMinEnergy( 0.0*GeV );                         71   SetMinEnergy( 0.0*GeV );
 72   SetMaxEnergy( 100.*TeV );                        72   SetMaxEnergy( 100.*TeV );
 73   SetMinEnergy(1.e-6*eV);                          73   SetMinEnergy(1.e-6*eV);
 74                                                    74 
 75   theNuMu =  G4NeutrinoMu::NeutrinoMu();           75   theNuMu =  G4NeutrinoMu::NeutrinoMu();
 76   theANuMu =  G4AntiNeutrinoMu::AntiNeutrinoMu     76   theANuMu =  G4AntiNeutrinoMu::AntiNeutrinoMu();
 77                                                    77 
 78   fMnumu = 0.;                                     78   fMnumu = 0.; 
 79   fData = fMaster = false;                         79   fData = fMaster = false;
 80   InitialiseModel();                               80   InitialiseModel();  
 81                                                    81      
 82 }                                                  82 }
 83                                                    83 
 84                                                    84 
 85 G4NuMuNucleusNcModel::~G4NuMuNucleusNcModel()      85 G4NuMuNucleusNcModel::~G4NuMuNucleusNcModel()
 86 {}                                                 86 {}
 87                                                    87 
 88                                                    88 
 89 void G4NuMuNucleusNcModel::ModelDescription(st     89 void G4NuMuNucleusNcModel::ModelDescription(std::ostream& outFile) const
 90 {                                                  90 {
 91                                                    91 
 92     outFile << "G4NuMuNucleusNcModel is a neut     92     outFile << "G4NuMuNucleusNcModel is a neutrino-nucleus (neutral current) scattering\n"
 93             << "model which uses the standard      93             << "model which uses the standard model \n"
 94             << "transfer parameterization.  Th     94             << "transfer parameterization.  The model is fully relativistic\n";
 95                                                    95 
 96 }                                                  96 }
 97                                                    97 
 98 //////////////////////////////////////////////     98 /////////////////////////////////////////////////////////
 99 //                                                 99 //
100 // Read data from G4PARTICLEXSDATA (locally PA    100 // Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA)
101                                                   101 
102 void G4NuMuNucleusNcModel::InitialiseModel()      102 void G4NuMuNucleusNcModel::InitialiseModel()
103 {                                                 103 {
104   G4String pName  = "nu_mu";                      104   G4String pName  = "nu_mu";
105                                                   105   
106   G4int nSize(0), i(0), j(0), k(0);               106   G4int nSize(0), i(0), j(0), k(0);
107                                                   107 
108   if(!fData)                                      108   if(!fData)
109   {                                               109   { 
110 #ifdef G4MULTITHREADED                            110 #ifdef G4MULTITHREADED
111     G4MUTEXLOCK(&numuNucleusModel);               111     G4MUTEXLOCK(&numuNucleusModel);
112     if(!fData)                                    112     if(!fData)
113     {                                             113     { 
114 #endif                                            114 #endif     
115       fMaster = true;                             115       fMaster = true;
116 #ifdef G4MULTITHREADED                            116 #ifdef G4MULTITHREADED
117     }                                             117     }
118     G4MUTEXUNLOCK(&numuNucleusModel);             118     G4MUTEXUNLOCK(&numuNucleusModel);
119 #endif                                            119 #endif
120   }                                               120   }
121                                                   121 
122   if(fMaster)                                     122   if(fMaster)
123   {                                               123   {  
124     const char* path = G4FindDataDir("G4PARTIC << 124     char* path = getenv("G4PARTICLEXSDATA");
125     std::ostringstream ost1, ost2, ost3, ost4;    125     std::ostringstream ost1, ost2, ost3, ost4;
126     ost1 << path << "/" << "neutrino" << "/" <    126     ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraynckr";
127                                                   127 
128     std::ifstream filein1( ost1.str().c_str()     128     std::ifstream filein1( ost1.str().c_str() );
129                                                   129 
130     // filein.open("$PARTICLEXSDATA/");           130     // filein.open("$PARTICLEXSDATA/");
131                                                   131 
132     filein1>>nSize;                               132     filein1>>nSize;
133                                                   133 
134     for( k = 0; k < fNbin; ++k )                  134     for( k = 0; k < fNbin; ++k )
135     {                                             135     {
136       for( i = 0; i <= fNbin; ++i )               136       for( i = 0; i <= fNbin; ++i )
137       {                                           137       {
138         filein1 >> fNuMuXarrayKR[k][i];           138         filein1 >> fNuMuXarrayKR[k][i];
139         // G4cout<< fNuMuXarrayKR[k][i] << "      139         // G4cout<< fNuMuXarrayKR[k][i] << "  ";
140       }                                           140       }
141     }                                             141     }
142     // G4cout<<G4endl<<G4endl;                    142     // G4cout<<G4endl<<G4endl;
143                                                   143 
144     ost2 << path << "/" << "neutrino" << "/" <    144     ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrnckr";
145     std::ifstream  filein2( ost2.str().c_str()    145     std::ifstream  filein2( ost2.str().c_str() );
146                                                   146 
147     filein2>>nSize;                               147     filein2>>nSize;
148                                                   148 
149     for( k = 0; k < fNbin; ++k )                  149     for( k = 0; k < fNbin; ++k )
150     {                                             150     {
151       for( i = 0; i < fNbin; ++i )                151       for( i = 0; i < fNbin; ++i )
152       {                                           152       {
153         filein2 >> fNuMuXdistrKR[k][i];           153         filein2 >> fNuMuXdistrKR[k][i];
154         // G4cout<< fNuMuXdistrKR[k][i] << "      154         // G4cout<< fNuMuXdistrKR[k][i] << "  ";
155       }                                           155       }
156     }                                             156     }
157     // G4cout<<G4endl<<G4endl;                    157     // G4cout<<G4endl<<G4endl;
158                                                   158 
159     ost3 << path << "/" << "neutrino" << "/" <    159     ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraynckr";
160     std::ifstream  filein3( ost3.str().c_str()    160     std::ifstream  filein3( ost3.str().c_str() );
161                                                   161 
162     filein3>>nSize;                               162     filein3>>nSize;
163                                                   163 
164     for( k = 0; k < fNbin; ++k )                  164     for( k = 0; k < fNbin; ++k )
165     {                                             165     {
166       for( i = 0; i <= fNbin; ++i )               166       for( i = 0; i <= fNbin; ++i )
167       {                                           167       {
168         for( j = 0; j <= fNbin; ++j )             168         for( j = 0; j <= fNbin; ++j )
169         {                                         169         {
170           filein3 >> fNuMuQarrayKR[k][i][j];      170           filein3 >> fNuMuQarrayKR[k][i][j];
171           // G4cout<< fNuMuQarrayKR[k][i][j] <    171           // G4cout<< fNuMuQarrayKR[k][i][j] << "  ";
172         }                                         172         }
173       }                                           173       }
174     }                                             174     }
175     // G4cout<<G4endl<<G4endl;                    175     // G4cout<<G4endl<<G4endl;
176                                                   176 
177     ost4 << path << "/" << "neutrino" << "/" <    177     ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrnckr";
178     std::ifstream  filein4( ost4.str().c_str()    178     std::ifstream  filein4( ost4.str().c_str() );
179                                                   179 
180     filein4>>nSize;                               180     filein4>>nSize;
181                                                   181 
182     for( k = 0; k < fNbin; ++k )                  182     for( k = 0; k < fNbin; ++k )
183     {                                             183     {
184       for( i = 0; i <= fNbin; ++i )               184       for( i = 0; i <= fNbin; ++i )
185       {                                           185       {
186         for( j = 0; j < fNbin; ++j )              186         for( j = 0; j < fNbin; ++j )
187         {                                         187         {
188           filein4 >> fNuMuQdistrKR[k][i][j];      188           filein4 >> fNuMuQdistrKR[k][i][j];
189           // G4cout<< fNuMuQdistrKR[k][i][j] <    189           // G4cout<< fNuMuQdistrKR[k][i][j] << "  ";
190         }                                         190         }
191       }                                           191       }
192     }                                             192     }
193     fData = true;                                 193     fData = true;
194   }                                               194   }
195 }                                                 195 }
196                                                   196 
197 //////////////////////////////////////////////    197 /////////////////////////////////////////////////////////
198                                                   198 
199 G4bool G4NuMuNucleusNcModel::IsApplicable(cons    199 G4bool G4NuMuNucleusNcModel::IsApplicable(const G4HadProjectile & aPart, 
200                   G4Nucleus & )                << 200                  G4Nucleus & targetNucleus)
201 {                                                 201 {
202   G4bool result  = false;                         202   G4bool result  = false;
203   G4String pName = aPart.GetDefinition()->GetP    203   G4String pName = aPart.GetDefinition()->GetParticleName();
204   G4double energy = aPart.GetTotalEnergy();       204   G4double energy = aPart.GetTotalEnergy();
205                                                   205   
206   if(  pName == "nu_mu" // || pName == "anti_n    206   if(  pName == "nu_mu" // || pName == "anti_nu_mu"   ) 
207         &&                                        207         &&
208         energy > fMinNuEnergy                     208         energy > fMinNuEnergy                                )
209   {                                               209   {
210     result = true;                                210     result = true;
211   }                                               211   }
                                                   >> 212   G4int Z = targetNucleus.GetZ_asInt();
                                                   >> 213         Z *= 1;
212                                                   214 
213   return result;                                  215   return result;
214 }                                                 216 }
215                                                   217 
216 /////////////////////////////////////////// Cl    218 /////////////////////////////////////////// ClusterDecay ////////////////////////////////////////////////////////////
217 //                                                219 //
218 //                                                220 //
219                                                   221 
220 G4HadFinalState* G4NuMuNucleusNcModel::ApplyYo    222 G4HadFinalState* G4NuMuNucleusNcModel::ApplyYourself(
221      const G4HadProjectile& aTrack, G4Nucleus&    223      const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
222 {                                                 224 {
223   theParticleChange.Clear();                      225   theParticleChange.Clear();
224   fProton = f2p2h = fBreak = false;               226   fProton = f2p2h = fBreak = false;
225   const G4HadProjectile* aParticle = &aTrack;     227   const G4HadProjectile* aParticle = &aTrack;
226   G4double energy = aParticle->GetTotalEnergy(    228   G4double energy = aParticle->GetTotalEnergy();
227                                                   229 
228   G4String pName  = aParticle->GetDefinition()    230   G4String pName  = aParticle->GetDefinition()->GetParticleName();
229                                                   231 
230   if( energy < fMinNuEnergy )                     232   if( energy < fMinNuEnergy ) 
231   {                                               233   {
232     theParticleChange.SetEnergyChange(energy);    234     theParticleChange.SetEnergyChange(energy);
233     theParticleChange.SetMomentumChange(aTrack    235     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
234     return &theParticleChange;                    236     return &theParticleChange;
235   }                                               237   }
236   SampleLVkr( aTrack, targetNucleus);             238   SampleLVkr( aTrack, targetNucleus);
237                                                   239 
238   if( fBreak == true || fEmu < fMnumu ) // ~5*    240   if( fBreak == true || fEmu < fMnumu ) // ~5*10^-6
239   {                                               241   {
240     // G4cout<<"ni, ";                            242     // G4cout<<"ni, ";
241     theParticleChange.SetEnergyChange(energy);    243     theParticleChange.SetEnergyChange(energy);
242     theParticleChange.SetMomentumChange(aTrack    244     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
243     return &theParticleChange;                    245     return &theParticleChange;
244   }                                               246   }
245                                                   247 
246   // LVs of initial state                         248   // LVs of initial state
247                                                   249 
248   G4LorentzVector lvp1 = aParticle->Get4Moment    250   G4LorentzVector lvp1 = aParticle->Get4Momentum();
249   G4LorentzVector lvt1( 0., 0., 0., fM1 );        251   G4LorentzVector lvt1( 0., 0., 0., fM1 );
250   G4double mPip = G4ParticleTable::GetParticle    252   G4double mPip = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass();
251                                                   253 
252   // 1-pi by fQtransfer && nu-energy              254   // 1-pi by fQtransfer && nu-energy
253   G4LorentzVector lvpip1( 0., 0., 0., mPip );     255   G4LorentzVector lvpip1( 0., 0., 0., mPip );
254   G4LorentzVector lvsum, lv2, lvX;                256   G4LorentzVector lvsum, lv2, lvX;
255   G4ThreeVector eP;                               257   G4ThreeVector eP;
256   G4double cost(1.), sint(0.), phi(0.), muMom(    258   G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.);
257   G4DynamicParticle* aLept = nullptr; // lepto    259   G4DynamicParticle* aLept = nullptr; // lepton lv
258                                                   260 
259   G4int Z = targetNucleus.GetZ_asInt();           261   G4int Z = targetNucleus.GetZ_asInt();
260   G4int A = targetNucleus.GetA_asInt();           262   G4int A = targetNucleus.GetA_asInt();
261   G4double  mTarg = targetNucleus.AtomicMass(A    263   G4double  mTarg = targetNucleus.AtomicMass(A,Z);
262   G4int pdgP(0), qB(0);                           264   G4int pdgP(0), qB(0);
263   // G4double mSum = G4ParticleTable::GetParti    265   // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip;
264                                                   266 
265   G4int iPi     = GetOnePionIndex(energy);        267   G4int iPi     = GetOnePionIndex(energy);
266   G4double p1pi = GetNuMuOnePionProb( iPi, ene    268   G4double p1pi = GetNuMuOnePionProb( iPi, energy);
267                                                   269 
268   if( p1pi > G4UniformRand() && fCosTheta > 0.    270   if( p1pi > G4UniformRand() && fCosTheta > 0.9  ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus
269   {                                               271   {
270     // lvsum = lvp1 + lvpip1;                     272     // lvsum = lvp1 + lvpip1;
271     lvsum = lvp1 + lvt1;                          273     lvsum = lvp1 + lvt1;
272     // cost = fCosThetaPi;                        274     // cost = fCosThetaPi;
273     cost = fCosTheta;                             275     cost = fCosTheta;
274     sint = std::sqrt( (1.0 - cost)*(1.0 + cost    276     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
275     phi  = G4UniformRand()*CLHEP::twopi;          277     phi  = G4UniformRand()*CLHEP::twopi;
276     eP   = G4ThreeVector( sint*std::cos(phi),     278     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
277                                                   279 
278     // muMom = sqrt(fEmuPi*fEmuPi-fMnumu*fMnum    280     // muMom = sqrt(fEmuPi*fEmuPi-fMnumu*fMnumu);
279     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);        281     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
280                                                   282 
281     eP *= muMom;                                  283     eP *= muMom;
282                                                   284 
283     // lv2 = G4LorentzVector( eP, fEmuPi );       285     // lv2 = G4LorentzVector( eP, fEmuPi );
284     lv2 = G4LorentzVector( eP, fEmu );            286     lv2 = G4LorentzVector( eP, fEmu );
285     lv2 = fLVl;                                   287     lv2 = fLVl;
286                                                   288 
287     lvX = lvsum - lv2;                            289     lvX = lvsum - lv2;
288     lvX = fLVh;                                   290     lvX = fLVh;
289     massX2 = lvX.m2();                            291     massX2 = lvX.m2();
290     G4double massX = lvX.m();                     292     G4double massX = lvX.m();
291     G4double massR = fLVt.m();                    293     G4double massR = fLVt.m();
292                                                   294 
293     // if ( massX2 <= 0. ) // vmg: very rarely    295     // if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
294     if ( massX2 <= fM1*fM1 ) // 9-3-20 vmg: ve    296     if ( massX2 <= fM1*fM1 ) // 9-3-20 vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
295       if ( lvX.e() <= fM1 ) // 9-3-20 vmg: ver    297       if ( lvX.e() <= fM1 ) // 9-3-20 vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
296     {                                             298     {
297       theParticleChange.SetEnergyChange(energy    299       theParticleChange.SetEnergyChange(energy);
298       theParticleChange.SetMomentumChange(aTra    300       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
299       return &theParticleChange;                  301       return &theParticleChange;
300     }                                             302     }
301     fW2 = massX2;                                 303     fW2 = massX2;
302                                                   304 
303     if(  pName == "nu_mu" )         aLept = ne    305     if(  pName == "nu_mu" )         aLept = new G4DynamicParticle( theNuMu, lv2 );  
304     else if( pName == "anti_nu_mu") aLept = ne    306     else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theANuMu,  lv2 );
305     else                                          307     else
306     {                                             308     {
307       theParticleChange.SetEnergyChange(energy    309       theParticleChange.SetEnergyChange(energy);
308       theParticleChange.SetMomentumChange(aTra    310       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
309       return &theParticleChange;                  311       return &theParticleChange;
310     }                                             312     } 
311                                                   313  
312     pdgP = 111;                                   314     pdgP = 111;
313                                                   315 
314     G4double eCut; // = fMpi + 0.5*(fMpi*fMpi     316     G4double eCut; // = fMpi + 0.5*(fMpi*fMpi - massX2)/mTarg; // massX -> fMpi
315                                                   317 
316     if( A > 1 )                                   318     if( A > 1 )
317     {                                             319     {
318       eCut = (fMpi + mTarg)*(fMpi + mTarg) - (    320       eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR);
319       eCut /= 2.*massR;                           321       eCut /= 2.*massR;
320       eCut += massX;                              322       eCut += massX;
321     }                                             323     }
322     else  eCut = fM1 + fMpi;                      324     else  eCut = fM1 + fMpi;
323                                                   325 
324     if ( lvX.e() > eCut ) // && sqrt( GetW2()     326     if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) // 
325     {                                             327     {
326       CoherentPion( lvX, pdgP, targetNucleus);    328       CoherentPion( lvX, pdgP, targetNucleus);
327     }                                             329     }
328     else                                          330     else
329     {                                             331     {
330       theParticleChange.SetEnergyChange(energy    332       theParticleChange.SetEnergyChange(energy);
331       theParticleChange.SetMomentumChange(aTra    333       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
332       return &theParticleChange;                  334       return &theParticleChange;
333     }                                             335     } 
334     theParticleChange.AddSecondary( aLept, fSe    336     theParticleChange.AddSecondary( aLept, fSecID );
335                                                   337 
336     return &theParticleChange;                    338     return &theParticleChange;
337   }                                               339   }
338   else // lepton part in lab                      340   else // lepton part in lab
339   {                                               341   { 
340     lvsum = lvp1 + lvt1;                          342     lvsum = lvp1 + lvt1;
341     cost = fCosTheta;                             343     cost = fCosTheta;
342     sint = std::sqrt( (1.0 - cost)*(1.0 + cost    344     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
343     phi  = G4UniformRand()*CLHEP::twopi;          345     phi  = G4UniformRand()*CLHEP::twopi;
344     eP   = G4ThreeVector( sint*std::cos(phi),     346     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
345                                                   347 
346     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);        348     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
347                                                   349 
348     eP *= muMom;                                  350     eP *= muMom;
349                                                   351 
350     lv2 = G4LorentzVector( eP, fEmu );            352     lv2 = G4LorentzVector( eP, fEmu );
351                                                   353 
352     lvX = lvsum - lv2;                            354     lvX = lvsum - lv2;
353                                                   355 
354     massX2 = lvX.m2();                            356     massX2 = lvX.m2();
355                                                   357 
356     if ( massX2 <= 0. ) // vmg: very rarely ~     358     if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved
357     {                                             359     {
358       theParticleChange.SetEnergyChange(energy    360       theParticleChange.SetEnergyChange(energy);
359       theParticleChange.SetMomentumChange(aTra    361       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
360       return &theParticleChange;                  362       return &theParticleChange;
361     }                                             363     }
362     fW2 = massX2;                                 364     fW2 = massX2;
363                                                   365 
364     aLept = new G4DynamicParticle( theNuMu, lv    366     aLept = new G4DynamicParticle( theNuMu, lv2 );  
365                                                   367   
366     theParticleChange.AddSecondary( aLept, fSe    368     theParticleChange.AddSecondary( aLept, fSecID );
367   }                                               369   }
368                                                   370 
369   // hadron part                                  371   // hadron part
370                                                   372 
371   fRecoil  = nullptr;                             373   fRecoil  = nullptr;
372   fCascade = false;                               374   fCascade = false;
373   fString  = false;                               375   fString  = false;
374                                                   376   
375   if( A == 1 )                                    377   if( A == 1 )
376   {                                               378   {
377     qB = 1;                                       379     qB = 1;
378                                                   380 
379     // if( G4UniformRand() > 0.1 ) //  > 0.999    381     // if( G4UniformRand() > 0.1 ) //  > 0.9999 ) // > 0.0001 ) //
380     {                                             382     {
381       ClusterDecay( lvX, qB );                    383       ClusterDecay( lvX, qB );
382     }                                             384     }
383     return &theParticleChange;                    385     return &theParticleChange;
384   }                                               386   }
385   G4Nucleus recoil;                               387   G4Nucleus recoil;
386   G4double rM(0.), ratio = G4double(Z)/G4doubl    388   G4double rM(0.), ratio = G4double(Z)/G4double(A);
387                                                   389 
388   if( ratio > G4UniformRand() ) // proton is e    390   if( ratio > G4UniformRand() ) // proton is excited
389   {                                               391   {
390     fProton = true;                               392     fProton = true;
391     recoil = G4Nucleus(A-1,Z-1);                  393     recoil = G4Nucleus(A-1,Z-1);
392     fRecoil = &recoil;                            394     fRecoil = &recoil;
393     rM = recoil.AtomicMass(A-1,Z-1);              395     rM = recoil.AtomicMass(A-1,Z-1);
394                                                   396 
395     fMt = G4ParticleTable::GetParticleTable()-    397     fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass()
396           + G4ParticleTable::GetParticleTable(    398           + G4ParticleTable::GetParticleTable()->FindParticle(111)->GetPDGMass();
397   }                                               399   }
398   else // excited neutron                         400   else // excited neutron
399   {                                               401   {
400     fProton = false;                              402     fProton = false;
401     recoil = G4Nucleus(A-1,Z);                    403     recoil = G4Nucleus(A-1,Z);
402     fRecoil = &recoil;                            404     fRecoil = &recoil;
403     rM = recoil.AtomicMass(A-1,Z);                405     rM = recoil.AtomicMass(A-1,Z);
404                                                   406 
405     fMt = G4ParticleTable::GetParticleTable()-    407     fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass()
406           + G4ParticleTable::GetParticleTable(    408           + G4ParticleTable::GetParticleTable()->FindParticle(111)->GetPDGMass(); 
407   }                                               409   }
408   // G4int       index = GetEnergyIndex(energy    410   // G4int       index = GetEnergyIndex(energy);
409   G4int nepdg = aParticle->GetDefinition()->Ge    411   G4int nepdg = aParticle->GetDefinition()->GetPDGEncoding();
410                                                   412 
411   G4double qeTotRat; //  = GetNuMuQeTotRat(ind    413   G4double qeTotRat; //  = GetNuMuQeTotRat(index, energy);
412   qeTotRat = CalculateQEratioA( Z, A, energy,     414   qeTotRat = CalculateQEratioA( Z, A, energy, nepdg);
413                                                   415 
414   G4ThreeVector dX = (lvX.vect()).unit();         416   G4ThreeVector dX = (lvX.vect()).unit();
415   G4double eX   = lvX.e();  // excited nucleon    417   G4double eX   = lvX.e();  // excited nucleon
416   G4double mX   = sqrt(massX2);                   418   G4double mX   = sqrt(massX2);
417                                                   419 
418   if( qeTotRat > G4UniformRand() || mX <= fMt     420   if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE
419   {                                               421   {  
420     fString = false;                              422     fString = false;
421                                                   423 
422     if( fProton )                                 424     if( fProton ) 
423     {                                             425     {  
424       fPDGencoding = 2212;                        426       fPDGencoding = 2212;
425       fMr =  proton_mass_c2;                      427       fMr =  proton_mass_c2;
426       recoil = G4Nucleus(A-1,Z-1);                428       recoil = G4Nucleus(A-1,Z-1);
427       fRecoil = &recoil;                          429       fRecoil = &recoil;
428       rM = recoil.AtomicMass(A-1,Z-1);            430       rM = recoil.AtomicMass(A-1,Z-1);
429     }                                             431     } 
430     else                                          432     else
431     {                                             433     {  
432       fPDGencoding = 2112;                        434       fPDGencoding = 2112;
433       fMr =   G4ParticleTable::GetParticleTabl    435       fMr =   G4ParticleTable::GetParticleTable()->
434   FindParticle(fPDGencoding)->GetPDGMass(); //    436   FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV;
435       recoil = G4Nucleus(A-1,Z);                  437       recoil = G4Nucleus(A-1,Z);
436       fRecoil = &recoil;                          438       fRecoil = &recoil;
437       rM = recoil.AtomicMass(A-1,Z);              439       rM = recoil.AtomicMass(A-1,Z);
438     }                                             440     } 
439     G4double eTh = fMr+0.5*(fMr*fMr-mX*mX)/rM;    441     G4double eTh = fMr+0.5*(fMr*fMr-mX*mX)/rM;
440                                                   442 
441     if(eX <= eTh) // vmg, very rarely out of k    443     if(eX <= eTh) // vmg, very rarely out of kinematics
442     {                                             444     {
443       theParticleChange.SetEnergyChange(energy    445       theParticleChange.SetEnergyChange(energy);
444       theParticleChange.SetMomentumChange(aTra    446       theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
445       return &theParticleChange;                  447       return &theParticleChange;
446     }                                             448     } 
447     FinalBarion( lvX, 0, fPDGencoding ); // p(    449     FinalBarion( lvX, 0, fPDGencoding ); // p(n)+deexcited recoil
448   }                                               450   }
449   else // if ( eX < 9500000.*GeV ) // < 25.*Ge    451   else // if ( eX < 9500000.*GeV ) // < 25.*GeV) //  < 95.*GeV ) // < 2.5*GeV ) //cluster decay
450   {                                               452   {  
451     if     (  fProton && pName == "nu_mu" )       453     if     (  fProton && pName == "nu_mu" )      qB =  1;
452     else if( !fProton && pName == "nu_mu" )       454     else if( !fProton && pName == "nu_mu" )      qB =  0;
453                                                   455 
454     ClusterDecay( lvX, qB );                      456     ClusterDecay( lvX, qB );
455   }                                               457   }
456   return &theParticleChange;                      458   return &theParticleChange;
457 }                                                 459 }
458                                                   460 
459                                                   461 
460 //////////////////////////////////////////////    462 /////////////////////////////////////////////////////////////////////
461 //////////////////////////////////////////////    463 ////////////////////////////////////////////////////////////////////
462 //////////////////////////////////////////////    464 ///////////////////////////////////////////////////////////////////
463                                                   465 
464 //////////////////////////////////////////////    466 /////////////////////////////////////////////////
465 //                                                467 //
466 // sample x, then Q2                              468 // sample x, then Q2
467                                                   469 
468 void G4NuMuNucleusNcModel::SampleLVkr(const G4    470 void G4NuMuNucleusNcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
469 {                                                 471 {
470   fBreak = false;                                 472   fBreak = false;
471   G4int A = targetNucleus.GetA_asInt(), iTer(0    473   G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100); 
472   G4int Z = targetNucleus.GetZ_asInt();           474   G4int Z = targetNucleus.GetZ_asInt(); 
473   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.)    475   G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z);
474   G4double cost(1.), sint(0.), phi(0.), muMom(    476   G4double cost(1.), sint(0.), phi(0.), muMom(0.); 
475   G4ThreeVector eP, bst;                          477   G4ThreeVector eP, bst;
476   const G4HadProjectile* aParticle = &aTrack;     478   const G4HadProjectile* aParticle = &aTrack;
477   G4LorentzVector lvp1 = aParticle->Get4Moment    479   G4LorentzVector lvp1 = aParticle->Get4Momentum();
478   nMom = NucleonMomentum( targetNucleus );        480   nMom = NucleonMomentum( targetNucleus );
479                                                   481 
480   if( A == 1 || nMom == 0. ) // hydrogen, no F    482   if( A == 1 || nMom == 0. ) // hydrogen, no Fermi motion ???
481   {                                               483   {
482     fNuEnergy = aParticle->GetTotalEnergy();      484     fNuEnergy = aParticle->GetTotalEnergy();
483     iTer = 0;                                     485     iTer = 0;
484                                                   486 
485     do                                            487     do
486     {                                             488     {
487       fXsample = SampleXkr(fNuEnergy);            489       fXsample = SampleXkr(fNuEnergy);
488       fQtransfer = SampleQkr(fNuEnergy, fXsamp    490       fQtransfer = SampleQkr(fNuEnergy, fXsample);
489       fQ2 = fQtransfer*fQtransfer;                491       fQ2 = fQtransfer*fQtransfer;
490                                                   492 
491      if( fXsample > 0. )                          493      if( fXsample > 0. )
492       {                                           494       {
493         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; //    495         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
494         fEmu = fNuEnergy - fQ2/2./fM1/fXsample    496         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
495       }                                           497       }
496       else                                        498       else
497       {                                           499       {
498         fW2 = fM1*fM1;                            500         fW2 = fM1*fM1;
499         fEmu = fNuEnergy;                         501         fEmu = fNuEnergy;
500       }                                           502       }
501       e3 = fNuEnergy + fM1 - fEmu;                503       e3 = fNuEnergy + fM1 - fEmu;
502                                                   504 
503       // if( e3 < sqrt(fW2) )  G4cout<<"energy    505       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl; // vmg ~10^-5 for NC
504                                                   506     
505       pMu2 = fEmu*fEmu - fMnumu*fMnumu;           507       pMu2 = fEmu*fEmu - fMnumu*fMnumu;
506       pX2  = e3*e3 - fW2;                         508       pX2  = e3*e3 - fW2;
507                                                   509 
508       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2    510       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
509       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);       511       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
510       iTer++;                                     512       iTer++;
511     }                                             513     }
512     while( ( abs(fCosTheta) > 1. || fEmu < fMn    514     while( ( abs(fCosTheta) > 1. || fEmu < fMnumu ) && iTer < iTerMax );
513                                                   515 
514     if( iTer >= iTerMax ) { fBreak = true; ret    516     if( iTer >= iTerMax ) { fBreak = true; return; }
515                                                   517 
516     if( abs(fCosTheta) > 1.) // vmg: due to bi    518     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
517     {                                             519     { 
518       G4cout<<"H2: fCosTheta = "<<fCosTheta<<"    520       G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
519       // fCosTheta = -1. + 2.*G4UniformRand();    521       // fCosTheta = -1. + 2.*G4UniformRand(); 
520       if(fCosTheta < -1.) fCosTheta = -1.;        522       if(fCosTheta < -1.) fCosTheta = -1.;
521       if(fCosTheta >  1.) fCosTheta =  1.;        523       if(fCosTheta >  1.) fCosTheta =  1.;
522     }                                             524     }
523     // LVs                                        525     // LVs
524                                                   526 
525     G4LorentzVector lvt1  = G4LorentzVector( 0    527     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
526     G4LorentzVector lvsum = lvp1 + lvt1;          528     G4LorentzVector lvsum = lvp1 + lvt1;
527                                                   529 
528     cost = fCosTheta;                             530     cost = fCosTheta;
529     sint = std::sqrt( (1.0 - cost)*(1.0 + cost    531     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
530     phi  = G4UniformRand()*CLHEP::twopi;          532     phi  = G4UniformRand()*CLHEP::twopi;
531     eP   = G4ThreeVector( sint*std::cos(phi),     533     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
532     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);        534     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
533     eP *= muMom;                                  535     eP *= muMom;
534     fLVl = G4LorentzVector( eP, fEmu );           536     fLVl = G4LorentzVector( eP, fEmu );
535                                                   537 
536     fLVh = lvsum - fLVl;                          538     fLVh = lvsum - fLVl;
537     fLVt = G4LorentzVector( 0., 0., 0., 0. );     539     fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil
538   }                                               540   }
539   else // Fermi motion, Q2 in nucleon rest fra    541   else // Fermi motion, Q2 in nucleon rest frame
540   {                                               542   {
541     G4ThreeVector nMomDir = nMom*G4RandomDirec    543     G4ThreeVector nMomDir = nMom*G4RandomDirection();
542                                                   544 
543     if( !f2p2h ) // 1p1h                          545     if( !f2p2h ) // 1p1h
544     {                                             546     {
545       G4Nucleus recoil(A-1,Z);                    547       G4Nucleus recoil(A-1,Z);
546       rM = sqrt( recoil.AtomicMass(A-1,Z)*reco    548       rM = sqrt( recoil.AtomicMass(A-1,Z)*recoil.AtomicMass(A-1,Z) + nMom*nMom );
547       hM = tM - rM;                               549       hM = tM - rM;
548                                                   550 
549       fLVt = G4LorentzVector( nMomDir, sqrt( r    551       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
550       fLVh = G4LorentzVector(-nMomDir, sqrt( h    552       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) ); 
551     }                                             553     }
552     else // 2p2h                                  554     else // 2p2h
553     {                                             555     {
554       G4Nucleus recoil(A-2,Z-1);                  556       G4Nucleus recoil(A-2,Z-1);
555       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMo    557       rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1);
556       hM = tM - rM;                               558       hM = tM - rM;
557                                                   559 
558       fLVt = G4LorentzVector( nMomDir, sqrt( r    560       fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) );
559       fLVh = G4LorentzVector(-nMomDir, sqrt( h    561       fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) ); 
560     }                                             562     }
561     // G4cout<<hM<<", ";                          563     // G4cout<<hM<<", ";
562     // bst = fLVh.boostVector(); // 9-3-20        564     // bst = fLVh.boostVector(); // 9-3-20
563                                                   565 
564     // lvp1.boost(-bst); // 9-3-20 -> nucleon     566     // lvp1.boost(-bst); // 9-3-20 -> nucleon rest system, where Q2 transfer is ???
565                                                   567 
566     fNuEnergy  = lvp1.e();                        568     fNuEnergy  = lvp1.e();
567     iTer = 0;                                     569     iTer = 0;
568                                                   570 
569     do                                            571     do
570     {                                             572     {
571       fXsample = SampleXkr(fNuEnergy);            573       fXsample = SampleXkr(fNuEnergy);
572       fQtransfer = SampleQkr(fNuEnergy, fXsamp    574       fQtransfer = SampleQkr(fNuEnergy, fXsample);
573       fQ2 = fQtransfer*fQtransfer;                575       fQ2 = fQtransfer*fQtransfer;
574                                                   576 
575       if( fXsample > 0. )                         577       if( fXsample > 0. )
576       {                                           578       {
577         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; //    579         fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass
578         fEmu = fNuEnergy - fQ2/2./fM1/fXsample    580         fEmu = fNuEnergy - fQ2/2./fM1/fXsample;
579       }                                           581       }
580       else                                        582       else
581       {                                           583       {
582         fW2 = fM1*fM1;                            584         fW2 = fM1*fM1;
583         fEmu = fNuEnergy;                         585         fEmu = fNuEnergy;
584       }                                           586       }
585                                                   587 
586       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu    588       // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl;
587                                                   589 
588       e3 = fNuEnergy + fM1 - fEmu;                590       e3 = fNuEnergy + fM1 - fEmu;
589                                                   591 
590       // if( e3 < sqrt(fW2) )  G4cout<<"energy    592       // if( e3 < sqrt(fW2) )  G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl;
591                                                   593     
592       pMu2 = fEmu*fEmu - fMnumu*fMnumu;           594       pMu2 = fEmu*fEmu - fMnumu*fMnumu;
593       pX2  = e3*e3 - fW2;                         595       pX2  = e3*e3 - fW2;
594                                                   596 
595       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2    597       fCosTheta  = fNuEnergy*fNuEnergy  + pMu2 - pX2;
596       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);       598       fCosTheta /= 2.*fNuEnergy*sqrt(pMu2);
597       iTer++;                                     599       iTer++;
598     }                                             600     }
599     while( ( abs(fCosTheta) > 1. || fEmu < fMn    601     while( ( abs(fCosTheta) > 1. || fEmu < fMnumu ) && iTer < iTerMax );
600                                                   602 
601     if( iTer >= iTerMax ) { fBreak = true; ret    603     if( iTer >= iTerMax ) { fBreak = true; return; }
602                                                   604 
603     if( abs(fCosTheta) > 1.) // vmg: due to bi    605     if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ...
604     {                                             606     { 
605       G4cout<<"FM: fCosTheta = "<<fCosTheta<<"    607       G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl;
606       // fCosTheta = -1. + 2.*G4UniformRand();    608       // fCosTheta = -1. + 2.*G4UniformRand(); 
607       if(fCosTheta < -1.) fCosTheta = -1.;        609       if(fCosTheta < -1.) fCosTheta = -1.;
608       if(fCosTheta >  1.) fCosTheta =  1.;        610       if(fCosTheta >  1.) fCosTheta =  1.;
609     }                                             611     }
610     // LVs                                        612     // LVs
611     G4LorentzVector lvt1  = G4LorentzVector( 0    613     G4LorentzVector lvt1  = G4LorentzVector( 0., 0., 0., fM1 );
612     G4LorentzVector lvsum = lvp1 + lvt1;          614     G4LorentzVector lvsum = lvp1 + lvt1;
613                                                   615 
614     cost = fCosTheta;                             616     cost = fCosTheta;
615     sint = std::sqrt( (1.0 - cost)*(1.0 + cost    617     sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
616     phi  = G4UniformRand()*CLHEP::twopi;          618     phi  = G4UniformRand()*CLHEP::twopi;
617     eP   = G4ThreeVector( sint*std::cos(phi),     619     eP   = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost );
618     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);        620     muMom = sqrt(fEmu*fEmu-fMnumu*fMnumu);
619     eP *= muMom;                                  621     eP *= muMom;
620     fLVl = G4LorentzVector( eP, fEmu );           622     fLVl = G4LorentzVector( eP, fEmu );
621     fLVh = lvsum - fLVl;                          623     fLVh = lvsum - fLVl;
622     // back to lab system                         624     // back to lab system
623     // fLVl.boost(bst); // 9-3-20                 625     // fLVl.boost(bst); // 9-3-20
624     // fLVh.boost(bst); // 9-3-20                 626     // fLVh.boost(bst); // 9-3-20
625   }                                               627   }
626   //G4cout<<iTer<<", "<<fBreak<<"; ";             628   //G4cout<<iTer<<", "<<fBreak<<"; ";
627 }                                                 629 }
628                                                   630 
629 //                                                631 //
630 //                                                632 //
631 ///////////////////////////                       633 ///////////////////////////
632                                                   634