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 ]

Diff markup

Differences between /processes/hadronic/models/lepto_nuclear/src/G4ANuMuNucleusCcModel.cc (Version 11.3.0) and /processes/hadronic/models/lepto_nuclear/src/G4ANuMuNucleusCcModel.cc (Version 10.7.p2)


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