Geant4 Cross Reference |
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: G4NuMuNucleusCcModel.cc 91806 2015-08- 26 // $Id: G4NuMuNucleusCcModel.cc 91806 2015-08-06 12:20:45Z gcosmo $ 27 // 27 // 28 // Geant4 Header : G4NuMuNucleusCcModel 28 // Geant4 Header : G4NuMuNucleusCcModel 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 "G4NuMuNucleusCcModel.hh" 37 #include "G4NuMuNucleusCcModel.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 const G4int G4NuMuNucleusCcModel::fResNumber = 6; >> 90 >> 91 const G4double G4NuMuNucleusCcModel::fResMass[6] = // [fResNumber] = >> 92 {2190., 1920., 1700., 1600., 1440., 1232. }; >> 93 >> 94 const G4int G4NuMuNucleusCcModel::fClustNumber = 4; >> 95 >> 96 const G4double G4NuMuNucleusCcModel::fMesMass[4] = {1260., 980., 770., 139.57}; >> 97 const G4int G4NuMuNucleusCcModel::fMesPDG[4] = {20213, 9000211, 213, 211}; >> 98 >> 99 // const G4double G4NuMuNucleusCcModel::fBarMass[4] = {1905., 1600., 1232., 939.57}; >> 100 // const G4int G4NuMuNucleusCcModel::fBarPDG[4] = {2226, 32224, 2224, 2212}; >> 101 >> 102 const G4double G4NuMuNucleusCcModel::fBarMass[4] = {1700., 1600., 1232., 939.57}; >> 103 const G4int G4NuMuNucleusCcModel::fBarPDG[4] = {12224, 32224, 2224, 2212}; >> 104 >> 105 const G4double G4NuMuNucleusCcModel::fNuMuEnergyLogVector[50] = { >> 106 115.603, 133.424, 153.991, 177.729, 205.126, 236.746, 273.24, 315.361, 363.973, 420.08, 484.836, 559.573, 645.832, >> 107 745.387, 860.289, 992.903, 1145.96, 1322.61, 1526.49, 1761.8, 2033.38, 2346.83, 2708.59, 3126.12, 3608.02, 4164.19, >> 108 4806.1, 5546.97, 6402.04, 7388.91, 8527.92, 9842.5, 11359.7, 13110.8, 15131.9, 17464.5, 20156.6, 23263.8, 26849.9, >> 109 30988.8, 35765.7, 41279, 47642.2, 54986.3, 63462.4, 73245.2, 84536, 97567.2, 112607, 129966 }; >> 110 >> 111 >> 112 G4double G4NuMuNucleusCcModel::fNuMuXarrayKR[50][51] = {{1.0}}; >> 113 G4double G4NuMuNucleusCcModel::fNuMuXdistrKR[50][50] = {{1.0}}; >> 114 G4double G4NuMuNucleusCcModel::fNuMuQarrayKR[50][51][51] = {{{1.0}}}; >> 115 G4double G4NuMuNucleusCcModel::fNuMuQdistrKR[50][51][50] = {{{1.0}}}; >> 116 89 #ifdef G4MULTITHREADED 117 #ifdef G4MULTITHREADED 90 G4Mutex G4NuMuNucleusCcModel::numuNucleusM 118 G4Mutex G4NuMuNucleusCcModel::numuNucleusModel = G4MUTEX_INITIALIZER; 91 #endif 119 #endif 92 120 93 121 94 G4NuMuNucleusCcModel::G4NuMuNucleusCcModel(con 122 G4NuMuNucleusCcModel::G4NuMuNucleusCcModel(const G4String& name) 95 : G4NeutrinoNucleusModel(name) 123 : G4NeutrinoNucleusModel(name) 96 { 124 { 97 fData = fMaster = false; 125 fData = fMaster = false; 98 InitialiseModel(); 126 InitialiseModel(); 99 } 127 } 100 128 101 129 102 G4NuMuNucleusCcModel::~G4NuMuNucleusCcModel() 130 G4NuMuNucleusCcModel::~G4NuMuNucleusCcModel() 103 {} 131 {} 104 132 105 133 106 void G4NuMuNucleusCcModel::ModelDescription(st 134 void G4NuMuNucleusCcModel::ModelDescription(std::ostream& outFile) const 107 { 135 { 108 136 109 outFile << "G4NuMuNucleusCcModel is a neut 137 outFile << "G4NuMuNucleusCcModel is a neutrino-nucleus (charge current) scattering\n" 110 << "model which uses the standard 138 << "model which uses the standard model \n" 111 << "transfer parameterization. Th 139 << "transfer parameterization. The model is fully relativistic\n"; 112 140 113 } 141 } 114 142 115 ////////////////////////////////////////////// 143 ///////////////////////////////////////////////////////// 116 // 144 // 117 // Read data from G4PARTICLEXSDATA (locally PA 145 // Read data from G4PARTICLEXSDATA (locally PARTICLEXSDATA) 118 146 119 void G4NuMuNucleusCcModel::InitialiseModel() 147 void G4NuMuNucleusCcModel::InitialiseModel() 120 { 148 { 121 G4String pName = "nu_mu"; 149 G4String pName = "nu_mu"; 122 // G4String pName = "anti_nu_mu"; << 123 150 124 G4int nSize(0), i(0), j(0), k(0); 151 G4int nSize(0), i(0), j(0), k(0); 125 152 126 if(!fData) 153 if(!fData) 127 { 154 { 128 #ifdef G4MULTITHREADED 155 #ifdef G4MULTITHREADED 129 G4MUTEXLOCK(&numuNucleusModel); 156 G4MUTEXLOCK(&numuNucleusModel); 130 if(!fData) 157 if(!fData) 131 { 158 { 132 #endif 159 #endif 133 fMaster = true; 160 fMaster = true; 134 #ifdef G4MULTITHREADED 161 #ifdef G4MULTITHREADED 135 } 162 } 136 G4MUTEXUNLOCK(&numuNucleusModel); 163 G4MUTEXUNLOCK(&numuNucleusModel); 137 #endif 164 #endif 138 } 165 } 139 166 140 if(fMaster) 167 if(fMaster) 141 { 168 { 142 const char* path = G4FindDataDir("G4PARTIC << 169 char* path = getenv("G4PARTICLEXSDATA"); 143 std::ostringstream ost1, ost2, ost3, ost4; 170 std::ostringstream ost1, ost2, ost3, ost4; 144 ost1 << path << "/" << "neutrino" << "/" < 171 ost1 << path << "/" << "neutrino" << "/" << pName << "/xarraycckr"; 145 172 146 std::ifstream filein1( ost1.str().c_str() 173 std::ifstream filein1( ost1.str().c_str() ); 147 174 148 // filein.open("$PARTICLEXSDATA/"); 175 // filein.open("$PARTICLEXSDATA/"); 149 176 150 filein1>>nSize; 177 filein1>>nSize; 151 178 152 for( k = 0; k < fNbin; ++k ) 179 for( k = 0; k < fNbin; ++k ) 153 { 180 { 154 for( i = 0; i <= fNbin; ++i ) 181 for( i = 0; i <= fNbin; ++i ) 155 { 182 { 156 filein1 >> fNuMuXarrayKR[k][i]; 183 filein1 >> fNuMuXarrayKR[k][i]; 157 // G4cout<< fNuMuXarrayKR[k][i] << " 184 // G4cout<< fNuMuXarrayKR[k][i] << " "; 158 } 185 } 159 } 186 } 160 // G4cout<<G4endl<<G4endl; 187 // G4cout<<G4endl<<G4endl; 161 188 162 ost2 << path << "/" << "neutrino" << "/" < 189 ost2 << path << "/" << "neutrino" << "/" << pName << "/xdistrcckr"; 163 std::ifstream filein2( ost2.str().c_str() 190 std::ifstream filein2( ost2.str().c_str() ); 164 191 165 filein2>>nSize; 192 filein2>>nSize; 166 193 167 for( k = 0; k < fNbin; ++k ) 194 for( k = 0; k < fNbin; ++k ) 168 { 195 { 169 for( i = 0; i < fNbin; ++i ) 196 for( i = 0; i < fNbin; ++i ) 170 { 197 { 171 filein2 >> fNuMuXdistrKR[k][i]; 198 filein2 >> fNuMuXdistrKR[k][i]; 172 // G4cout<< fNuMuXdistrKR[k][i] << " 199 // G4cout<< fNuMuXdistrKR[k][i] << " "; 173 } 200 } 174 } 201 } 175 // G4cout<<G4endl<<G4endl; 202 // G4cout<<G4endl<<G4endl; 176 203 177 ost3 << path << "/" << "neutrino" << "/" < 204 ost3 << path << "/" << "neutrino" << "/" << pName << "/q2arraycckr"; 178 std::ifstream filein3( ost3.str().c_str() 205 std::ifstream filein3( ost3.str().c_str() ); 179 206 180 filein3>>nSize; 207 filein3>>nSize; 181 208 182 for( k = 0; k < fNbin; ++k ) 209 for( k = 0; k < fNbin; ++k ) 183 { 210 { 184 for( i = 0; i <= fNbin; ++i ) 211 for( i = 0; i <= fNbin; ++i ) 185 { 212 { 186 for( j = 0; j <= fNbin; ++j ) 213 for( j = 0; j <= fNbin; ++j ) 187 { 214 { 188 filein3 >> fNuMuQarrayKR[k][i][j]; 215 filein3 >> fNuMuQarrayKR[k][i][j]; 189 // G4cout<< fNuMuQarrayKR[k][i][j] < 216 // G4cout<< fNuMuQarrayKR[k][i][j] << " "; 190 } 217 } 191 } 218 } 192 } 219 } 193 // G4cout<<G4endl<<G4endl; 220 // G4cout<<G4endl<<G4endl; 194 221 195 ost4 << path << "/" << "neutrino" << "/" < 222 ost4 << path << "/" << "neutrino" << "/" << pName << "/q2distrcckr"; 196 std::ifstream filein4( ost4.str().c_str() 223 std::ifstream filein4( ost4.str().c_str() ); 197 224 198 filein4>>nSize; 225 filein4>>nSize; 199 226 200 for( k = 0; k < fNbin; ++k ) 227 for( k = 0; k < fNbin; ++k ) 201 { 228 { 202 for( i = 0; i <= fNbin; ++i ) 229 for( i = 0; i <= fNbin; ++i ) 203 { 230 { 204 for( j = 0; j < fNbin; ++j ) 231 for( j = 0; j < fNbin; ++j ) 205 { 232 { 206 filein4 >> fNuMuQdistrKR[k][i][j]; 233 filein4 >> fNuMuQdistrKR[k][i][j]; 207 // G4cout<< fNuMuQdistrKR[k][i][j] < 234 // G4cout<< fNuMuQdistrKR[k][i][j] << " "; 208 } 235 } 209 } 236 } 210 } 237 } 211 fData = true; 238 fData = true; 212 } 239 } 213 } 240 } 214 241 215 ////////////////////////////////////////////// 242 ///////////////////////////////////////////////////////// 216 243 217 G4bool G4NuMuNucleusCcModel::IsApplicable(cons 244 G4bool G4NuMuNucleusCcModel::IsApplicable(const G4HadProjectile & aPart, 218 G4Nucleus & ) << 245 G4Nucleus & targetNucleus) 219 { 246 { 220 G4bool result = false; 247 G4bool result = false; 221 G4String pName = aPart.GetDefinition()->GetP 248 G4String pName = aPart.GetDefinition()->GetParticleName(); 222 G4double energy = aPart.GetTotalEnergy(); 249 G4double energy = aPart.GetTotalEnergy(); 223 250 224 if( pName == "nu_mu" // || pName == "anti_ << 251 if( pName == "nu_mu" // || pName == "anti_nu_mu" ) 225 && 252 && 226 energy > fMinNuEnergy 253 energy > fMinNuEnergy ) 227 { 254 { 228 result = true; 255 result = true; 229 } 256 } >> 257 G4int Z = targetNucleus.GetZ_asInt(); >> 258 Z *= 1; 230 259 231 return result; 260 return result; 232 } 261 } 233 262 234 /////////////////////////////////////////// Cl 263 /////////////////////////////////////////// ClusterDecay //////////////////////////////////////////////////////////// 235 // 264 // 236 // 265 // 237 266 238 G4HadFinalState* G4NuMuNucleusCcModel::ApplyYo 267 G4HadFinalState* G4NuMuNucleusCcModel::ApplyYourself( 239 const G4HadProjectile& aTrack, G4Nucleus& 268 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus) 240 { 269 { 241 theParticleChange.Clear(); 270 theParticleChange.Clear(); 242 fProton = f2p2h = fBreak = false; 271 fProton = f2p2h = fBreak = false; 243 fCascade = fString = false; << 244 fLVh = fLVl = fLVt = fLVcpi = G4LorentzVecto << 245 << 246 const G4HadProjectile* aParticle = &aTrack; 272 const G4HadProjectile* aParticle = &aTrack; 247 G4double energy = aParticle->GetTotalEnergy( 273 G4double energy = aParticle->GetTotalEnergy(); 248 274 249 G4String pName = aParticle->GetDefinition() 275 G4String pName = aParticle->GetDefinition()->GetParticleName(); 250 276 251 if( energy < fMinNuEnergy ) 277 if( energy < fMinNuEnergy ) 252 { 278 { 253 theParticleChange.SetEnergyChange(energy); 279 theParticleChange.SetEnergyChange(energy); 254 theParticleChange.SetMomentumChange(aTrack 280 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 255 return &theParticleChange; 281 return &theParticleChange; 256 } 282 } 257 << 258 SampleLVkr( aTrack, targetNucleus); 283 SampleLVkr( aTrack, targetNucleus); 259 284 260 if( fBreak == true || fEmu < fMu ) // ~5*10^ 285 if( fBreak == true || fEmu < fMu ) // ~5*10^-6 261 { 286 { 262 // G4cout<<"ni, "; 287 // G4cout<<"ni, "; 263 theParticleChange.SetEnergyChange(energy); 288 theParticleChange.SetEnergyChange(energy); 264 theParticleChange.SetMomentumChange(aTrack 289 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 265 return &theParticleChange; 290 return &theParticleChange; 266 } 291 } 267 292 268 // LVs of initial state 293 // LVs of initial state 269 294 270 G4LorentzVector lvp1 = aParticle->Get4Moment 295 G4LorentzVector lvp1 = aParticle->Get4Momentum(); 271 G4LorentzVector lvt1( 0., 0., 0., fM1 ); 296 G4LorentzVector lvt1( 0., 0., 0., fM1 ); 272 G4double mPip = G4ParticleTable::GetParticle 297 G4double mPip = G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass(); 273 298 274 // 1-pi by fQtransfer && nu-energy 299 // 1-pi by fQtransfer && nu-energy 275 G4LorentzVector lvpip1( 0., 0., 0., mPip ); 300 G4LorentzVector lvpip1( 0., 0., 0., mPip ); 276 G4LorentzVector lvsum, lv2, lvX; 301 G4LorentzVector lvsum, lv2, lvX; 277 G4ThreeVector eP; 302 G4ThreeVector eP; 278 G4double cost(1.), sint(0.), phi(0.), muMom( 303 G4double cost(1.), sint(0.), phi(0.), muMom(0.), massX2(0.), massX(0.), massR(0.), eCut(0.); 279 G4DynamicParticle* aLept = nullptr; // lepto 304 G4DynamicParticle* aLept = nullptr; // lepton lv 280 305 281 G4int Z = targetNucleus.GetZ_asInt(); 306 G4int Z = targetNucleus.GetZ_asInt(); 282 G4int A = targetNucleus.GetA_asInt(); 307 G4int A = targetNucleus.GetA_asInt(); 283 G4double mTarg = targetNucleus.AtomicMass(A 308 G4double mTarg = targetNucleus.AtomicMass(A,Z); 284 G4int pdgP(0), qB(0); 309 G4int pdgP(0), qB(0); 285 // G4double mSum = G4ParticleTable::GetParti 310 // G4double mSum = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() + mPip; 286 311 287 G4int iPi = GetOnePionIndex(energy); 312 G4int iPi = GetOnePionIndex(energy); 288 G4double p1pi = GetNuMuOnePionProb( iPi, ene 313 G4double p1pi = GetNuMuOnePionProb( iPi, energy); 289 314 290 if( p1pi > G4UniformRand() && fCosTheta > 0 << 315 if( p1pi > G4UniformRand() ) // && fQtransfer < 0.95*GeV ) // mu- & coherent pion + nucleus 291 { 316 { 292 // lvsum = lvp1 + lvpip1; 317 // lvsum = lvp1 + lvpip1; 293 lvsum = lvp1 + lvt1; 318 lvsum = lvp1 + lvt1; 294 // cost = fCosThetaPi; 319 // cost = fCosThetaPi; 295 cost = fCosTheta; 320 cost = fCosTheta; 296 sint = std::sqrt( (1.0 - cost)*(1.0 + cost 321 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) ); 297 phi = G4UniformRand()*CLHEP::twopi; 322 phi = G4UniformRand()*CLHEP::twopi; 298 eP = G4ThreeVector( sint*std::cos(phi), 323 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost ); 299 324 300 // muMom = sqrt(fEmuPi*fEmuPi-fMu*fMu); 325 // muMom = sqrt(fEmuPi*fEmuPi-fMu*fMu); 301 muMom = sqrt(fEmu*fEmu-fMu*fMu); 326 muMom = sqrt(fEmu*fEmu-fMu*fMu); 302 327 303 eP *= muMom; 328 eP *= muMom; 304 329 305 // lv2 = G4LorentzVector( eP, fEmuPi ); 330 // lv2 = G4LorentzVector( eP, fEmuPi ); 306 // lv2 = G4LorentzVector( eP, fEmu ); 331 // lv2 = G4LorentzVector( eP, fEmu ); 307 lv2 = fLVl; 332 lv2 = fLVl; 308 333 309 // lvX = lvsum - lv2; 334 // lvX = lvsum - lv2; 310 lvX = fLVh; 335 lvX = fLVh; 311 massX2 = lvX.m2(); 336 massX2 = lvX.m2(); 312 massX = lvX.m(); 337 massX = lvX.m(); 313 massR = fLVt.m(); 338 massR = fLVt.m(); 314 339 315 if ( massX2 <= 0. ) // vmg: very rarely ~ 340 if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved 316 { 341 { 317 fCascade = true; << 318 theParticleChange.SetEnergyChange(energy 342 theParticleChange.SetEnergyChange(energy); 319 theParticleChange.SetMomentumChange(aTra 343 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 320 return &theParticleChange; 344 return &theParticleChange; 321 } 345 } 322 fW2 = massX2; 346 fW2 = massX2; 323 347 324 if( pName == "nu_mu" ) aLept = ne 348 if( pName == "nu_mu" ) aLept = new G4DynamicParticle( theMuonMinus, lv2 ); 325 // else if( pName == "anti_nu_mu") aLept = << 349 else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus, lv2 ); 326 else 350 else 327 { 351 { 328 theParticleChange.SetEnergyChange(energy 352 theParticleChange.SetEnergyChange(energy); 329 theParticleChange.SetMomentumChange(aTra 353 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 330 return &theParticleChange; 354 return &theParticleChange; 331 } 355 } 332 if( pName == "nu_mu" ) pdgP = 211; 356 if( pName == "nu_mu" ) pdgP = 211; 333 // else pdgP = -211; << 357 else pdgP = -211; >> 358 334 // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mT 359 // eCut = fMpi + 0.5*(fMpi*fMpi-massX2)/mTarg; // massX -> fMpi 335 360 336 if( A > 1 ) << 361 eCut = (fMpi + mTarg)*(fMpi + mTarg) - (massX + massR)*(massX + massR); 337 { << 362 eCut /= 2.*massR; 338 eCut = (fMpi + mTarg)*(fMpi + mTarg) - ( << 363 eCut += massX; 339 eCut /= 2.*massR; << 340 eCut += massX; << 341 } << 342 else eCut = fM1 + fMpi; << 343 364 344 if ( lvX.e() > eCut ) // && sqrt( GetW2() 365 if ( lvX.e() > eCut ) // && sqrt( GetW2() ) < 1.4*GeV ) // 345 { 366 { 346 CoherentPion( lvX, pdgP, targetNucleus); 367 CoherentPion( lvX, pdgP, targetNucleus); 347 } 368 } 348 else 369 else 349 { 370 { 350 fCascade = true; << 351 theParticleChange.SetEnergyChange(energy 371 theParticleChange.SetEnergyChange(energy); 352 theParticleChange.SetMomentumChange(aTra 372 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 353 return &theParticleChange; 373 return &theParticleChange; 354 } 374 } 355 theParticleChange.AddSecondary( aLept, fSe << 375 theParticleChange.AddSecondary( aLept ); 356 376 357 return &theParticleChange; 377 return &theParticleChange; 358 } 378 } 359 else // lepton part in lab 379 else // lepton part in lab 360 { 380 { 361 lvsum = lvp1 + lvt1; 381 lvsum = lvp1 + lvt1; 362 cost = fCosTheta; 382 cost = fCosTheta; 363 sint = std::sqrt( (1.0 - cost)*(1.0 + cost 383 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) ); 364 phi = G4UniformRand()*CLHEP::twopi; 384 phi = G4UniformRand()*CLHEP::twopi; 365 eP = G4ThreeVector( sint*std::cos(phi), 385 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost ); 366 386 367 muMom = sqrt(fEmu*fEmu-fMu*fMu); 387 muMom = sqrt(fEmu*fEmu-fMu*fMu); 368 388 369 eP *= muMom; 389 eP *= muMom; 370 390 371 lv2 = G4LorentzVector( eP, fEmu ); 391 lv2 = G4LorentzVector( eP, fEmu ); 372 lv2 = fLVl; 392 lv2 = fLVl; 373 lvX = lvsum - lv2; 393 lvX = lvsum - lv2; 374 lvX = fLVh; 394 lvX = fLVh; 375 massX2 = lvX.m2(); 395 massX2 = lvX.m2(); 376 396 377 if ( massX2 <= 0. ) // vmg: very rarely ~ 397 if ( massX2 <= 0. ) // vmg: very rarely ~ (1-4)e-6 due to big Q2/x, to be improved 378 { 398 { 379 fCascade = true; << 380 theParticleChange.SetEnergyChange(energy 399 theParticleChange.SetEnergyChange(energy); 381 theParticleChange.SetMomentumChange(aTra 400 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 382 return &theParticleChange; 401 return &theParticleChange; 383 } 402 } 384 fW2 = massX2; 403 fW2 = massX2; 385 404 386 if( pName == "nu_mu" ) aLept = ne 405 if( pName == "nu_mu" ) aLept = new G4DynamicParticle( theMuonMinus, lv2 ); 387 // else if( pName == "anti_nu_mu") aLept = << 406 else if( pName == "anti_nu_mu") aLept = new G4DynamicParticle( theMuonPlus, lv2 ); 388 else 407 else 389 { 408 { 390 theParticleChange.SetEnergyChange(energy 409 theParticleChange.SetEnergyChange(energy); 391 theParticleChange.SetMomentumChange(aTra 410 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 392 return &theParticleChange; 411 return &theParticleChange; 393 } 412 } 394 theParticleChange.AddSecondary( aLept, fSe << 413 theParticleChange.AddSecondary( aLept ); 395 } 414 } 396 415 397 // hadron part 416 // hadron part 398 417 399 fRecoil = nullptr; 418 fRecoil = nullptr; >> 419 fCascade = false; >> 420 fString = false; 400 421 401 if( A == 1 ) 422 if( A == 1 ) 402 { 423 { 403 if( pName == "nu_mu" ) qB = 2; 424 if( pName == "nu_mu" ) qB = 2; 404 // else qB = 0; << 425 else qB = 0; 405 426 406 // if( G4UniformRand() > 0.1 ) // > 0.999 427 // if( G4UniformRand() > 0.1 ) // > 0.9999 ) // > 0.0001 ) // 407 { 428 { 408 ClusterDecay( lvX, qB ); 429 ClusterDecay( lvX, qB ); 409 } 430 } 410 return &theParticleChange; 431 return &theParticleChange; 411 } 432 } 412 /* 433 /* 413 // else 434 // else 414 { 435 { 415 if( pName == "nu_mu" ) pdgP = 211; 436 if( pName == "nu_mu" ) pdgP = 211; 416 else pdgP = -211; 437 else pdgP = -211; 417 438 418 439 419 if ( fQtransfer < 0.95*GeV ) // < 0.35*G 440 if ( fQtransfer < 0.95*GeV ) // < 0.35*GeV ) // 420 { 441 { 421 if( lvX.m() > mSum ) CoherentPion( lvX, pdgP 442 if( lvX.m() > mSum ) CoherentPion( lvX, pdgP, targetNucleus); 422 } 443 } 423 } 444 } 424 return &theParticleChange; 445 return &theParticleChange; 425 } 446 } 426 */ 447 */ 427 G4Nucleus recoil; 448 G4Nucleus recoil; 428 G4double rM(0.), ratio = G4double(Z)/G4doubl 449 G4double rM(0.), ratio = G4double(Z)/G4double(A); 429 450 430 if( ratio > G4UniformRand() ) // proton is e 451 if( ratio > G4UniformRand() ) // proton is excited 431 { 452 { 432 fProton = true; 453 fProton = true; 433 recoil = G4Nucleus(A-1,Z-1); 454 recoil = G4Nucleus(A-1,Z-1); 434 fRecoil = &recoil; 455 fRecoil = &recoil; 435 rM = recoil.AtomicMass(A-1,Z-1); 456 rM = recoil.AtomicMass(A-1,Z-1); 436 457 437 if( pName == "nu_mu" ) // (++) state -> p 458 if( pName == "nu_mu" ) // (++) state -> p + pi+ 438 { 459 { 439 fMt = G4ParticleTable::GetParticleTable( 460 fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() 440 + G4ParticleTable::GetParticleTable( 461 + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass(); 441 } 462 } 442 else // (0) state -> p + pi-, n + pi0 463 else // (0) state -> p + pi-, n + pi0 443 { 464 { 444 // fMt = G4ParticleTable::GetParticleTab << 465 fMt = G4ParticleTable::GetParticleTable()->FindParticle(2212)->GetPDGMass() 445 // + G4ParticleTable::GetParticleTab << 466 + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass(); 446 } 467 } 447 } 468 } 448 else // excited neutron 469 else // excited neutron 449 { 470 { 450 fProton = false; 471 fProton = false; 451 recoil = G4Nucleus(A-1,Z); 472 recoil = G4Nucleus(A-1,Z); 452 fRecoil = &recoil; 473 fRecoil = &recoil; 453 rM = recoil.AtomicMass(A-1,Z); 474 rM = recoil.AtomicMass(A-1,Z); 454 475 455 if( pName == "nu_mu" ) // (+) state -> n + 476 if( pName == "nu_mu" ) // (+) state -> n + pi+ 456 { 477 { 457 fMt = G4ParticleTable::GetParticleTable( 478 fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass() 458 + G4ParticleTable::GetParticleTable( 479 + G4ParticleTable::GetParticleTable()->FindParticle(211)->GetPDGMass(); 459 } 480 } 460 else // (-) state -> n + pi-, // n + pi0 481 else // (-) state -> n + pi-, // n + pi0 461 { 482 { 462 // fMt = G4ParticleTable::GetParticleTab << 483 fMt = G4ParticleTable::GetParticleTable()->FindParticle(2112)->GetPDGMass() 463 // + G4ParticleTable::GetParticleTab << 484 + G4ParticleTable::GetParticleTable()->FindParticle(-211)->GetPDGMass(); 464 } 485 } 465 } 486 } 466 // G4int index = GetEnergyIndex(energy << 487 G4int index = GetEnergyIndex(energy); 467 G4int nepdg = aParticle->GetDefinition()->Ge << 488 G4double qeTotRat = GetNuMuQeTotRat(index, energy); 468 << 469 G4double qeTotRat; // = GetNuMuQeTotRat(inde << 470 qeTotRat = CalculateQEratioA( Z, A, energy, << 471 489 472 G4ThreeVector dX = (lvX.vect()).unit(); 490 G4ThreeVector dX = (lvX.vect()).unit(); 473 G4double eX = lvX.e(); // excited nucleon 491 G4double eX = lvX.e(); // excited nucleon 474 G4double mX = sqrt(massX2); 492 G4double mX = sqrt(massX2); 475 // G4double pX = sqrt( eX*eX - mX*mX ); 493 // G4double pX = sqrt( eX*eX - mX*mX ); 476 // G4double sumE = eX + rM; 494 // G4double sumE = eX + rM; 477 495 478 if( qeTotRat > G4UniformRand() || mX <= fMt 496 if( qeTotRat > G4UniformRand() || mX <= fMt ) // || eX <= 1232.*MeV) // QE 479 { 497 { 480 fString = false; 498 fString = false; 481 499 482 if( fProton ) << 500 if( pName == "nu_mu" ) 483 { 501 { 484 fPDGencoding = 2212; 502 fPDGencoding = 2212; 485 fMr = proton_mass_c2; 503 fMr = proton_mass_c2; 486 recoil = G4Nucleus(A-1,Z-1); << 504 recoil = G4Nucleus(A-1,Z); 487 fRecoil = &recoil; 505 fRecoil = &recoil; 488 rM = recoil.AtomicMass(A-1,Z-1); << 506 rM = recoil.AtomicMass(A-1,Z); 489 } 507 } 490 else // if( pName == "anti_nu_mu" ) 508 else // if( pName == "anti_nu_mu" ) 491 { 509 { 492 fPDGencoding = 2112; 510 fPDGencoding = 2112; 493 fMr = G4ParticleTable::GetParticleTabl 511 fMr = G4ParticleTable::GetParticleTable()-> 494 FindParticle(fPDGencoding)->GetPDGMass(); // 512 FindParticle(fPDGencoding)->GetPDGMass(); // 939.5654133*MeV; 495 recoil = G4Nucleus(A-1,Z); << 513 recoil = G4Nucleus(A-1,Z-1); 496 fRecoil = &recoil; 514 fRecoil = &recoil; 497 rM = recoil.AtomicMass(A-1,Z); << 515 rM = recoil.AtomicMass(A-1,Z-1); 498 } 516 } 499 // sumE = eX + rM; 517 // sumE = eX + rM; 500 G4double eTh = fMr + 0.5*(fMr*fMr - mX*mX) 518 G4double eTh = fMr + 0.5*(fMr*fMr - mX*mX)/rM; 501 519 502 if( eX <= eTh ) // vmg, very rarely out of 520 if( eX <= eTh ) // vmg, very rarely out of kinematics 503 { 521 { 504 fString = true; << 505 theParticleChange.SetEnergyChange(energy 522 theParticleChange.SetEnergyChange(energy); 506 theParticleChange.SetMomentumChange(aTra 523 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 507 return &theParticleChange; 524 return &theParticleChange; 508 } 525 } 509 // FinalBarion( fLVh, 0, fPDGencoding ); / << 526 FinalBarion( fLVh, 0, fPDGencoding ); // p(n)+deexcited recoil 510 FinalBarion( lvX, 0, fPDGencoding ); // p( << 511 } 527 } 512 else // if ( eX < 9500000.*GeV ) // < 25.*G << 528 else if ( eX < 95000.*GeV ) // < 25.*GeV) // < 95.*GeV ) // < 2.5*GeV ) //cluster decay 513 { 529 { 514 if ( fProton && pName == "nu_mu" ) 530 if ( fProton && pName == "nu_mu" ) qB = 2; 515 // else if( fProton && pName == "anti_nu_ << 531 else if( fProton && pName == "anti_nu_mu" ) qB = 0; 516 else if( !fProton && pName == "nu_mu" ) 532 else if( !fProton && pName == "nu_mu" ) qB = 1; 517 // else if( !fProton && pName == "anti_nu_ << 533 else if( !fProton && pName == "anti_nu_mu" ) qB = -1; 518 << 519 534 >> 535 // if( G4UniformRand() > 0.1 ) >> 536 { 520 ClusterDecay( lvX, qB ); 537 ClusterDecay( lvX, qB ); >> 538 } >> 539 // else >> 540 { >> 541 if( pName == "nu_mu" ) pdgP = 211; >> 542 else pdgP = -211; >> 543 >> 544 if ( fQtransfer < 0.95*GeV ) // < 0.99*GeV ) // >> 545 { >> 546 // if( lvX.m() > mSum ) CoherentPion( lvX, pdgP, targetNucleus); >> 547 } >> 548 } 521 } 549 } >> 550 else // string >> 551 { >> 552 return &theParticleChange; >> 553 } 522 return &theParticleChange; 554 return &theParticleChange; 523 } 555 } 524 556 525 557 526 ////////////////////////////////////////////// 558 ///////////////////////////////////////////////////////////////////// 527 ////////////////////////////////////////////// 559 //////////////////////////////////////////////////////////////////// 528 ////////////////////////////////////////////// 560 /////////////////////////////////////////////////////////////////// 529 561 530 ////////////////////////////////////////////// 562 ///////////////////////////////////////////////// 531 // 563 // 532 // sample x, then Q2 564 // sample x, then Q2 533 565 534 void G4NuMuNucleusCcModel::SampleLVkr(const G4 566 void G4NuMuNucleusCcModel::SampleLVkr(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus) 535 { 567 { 536 fBreak = false; 568 fBreak = false; 537 G4int A = targetNucleus.GetA_asInt(), iTer(0 569 G4int A = targetNucleus.GetA_asInt(), iTer(0), iTerMax(100); 538 G4int Z = targetNucleus.GetZ_asInt(); 570 G4int Z = targetNucleus.GetZ_asInt(); 539 G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.) 571 G4double e3(0.), pMu2(0.), pX2(0.), nMom(0.), rM(0.), hM(0.), tM = targetNucleus.AtomicMass(A,Z); 540 G4double Ex(0.), ei(0.), nm2(0.); 572 G4double Ex(0.), ei(0.), nm2(0.); 541 G4double cost(1.), sint(0.), phi(0.), muMom( 573 G4double cost(1.), sint(0.), phi(0.), muMom(0.); 542 G4ThreeVector eP, bst; 574 G4ThreeVector eP, bst; 543 const G4HadProjectile* aParticle = &aTrack; 575 const G4HadProjectile* aParticle = &aTrack; 544 G4LorentzVector lvp1 = aParticle->Get4Moment 576 G4LorentzVector lvp1 = aParticle->Get4Momentum(); 545 577 546 if( A == 1 ) // hydrogen, no Fermi motion ?? 578 if( A == 1 ) // hydrogen, no Fermi motion ??? 547 { 579 { 548 fNuEnergy = aParticle->GetTotalEnergy(); 580 fNuEnergy = aParticle->GetTotalEnergy(); 549 iTer = 0; 581 iTer = 0; 550 582 551 do 583 do 552 { 584 { 553 fXsample = SampleXkr(fNuEnergy); 585 fXsample = SampleXkr(fNuEnergy); 554 fQtransfer = SampleQkr(fNuEnergy, fXsamp 586 fQtransfer = SampleQkr(fNuEnergy, fXsample); 555 fQ2 = fQtransfer*fQtransfer; 587 fQ2 = fQtransfer*fQtransfer; 556 588 557 if( fXsample > 0. ) 589 if( fXsample > 0. ) 558 { 590 { 559 fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // 591 fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass 560 fEmu = fNuEnergy - fQ2/2./fM1/fXsample 592 fEmu = fNuEnergy - fQ2/2./fM1/fXsample; 561 } 593 } 562 else 594 else 563 { 595 { 564 fW2 = fM1*fM1; 596 fW2 = fM1*fM1; 565 fEmu = fNuEnergy; 597 fEmu = fNuEnergy; 566 } 598 } 567 e3 = fNuEnergy + fM1 - fEmu; 599 e3 = fNuEnergy + fM1 - fEmu; 568 600 569 if( e3 < sqrt(fW2) ) G4cout<<"energyX = 601 if( e3 < sqrt(fW2) ) G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl; 570 602 571 pMu2 = fEmu*fEmu - fMu*fMu; 603 pMu2 = fEmu*fEmu - fMu*fMu; 572 << 573 if(pMu2 < 0.) { fBreak = true; return; } << 574 << 575 pX2 = e3*e3 - fW2; 604 pX2 = e3*e3 - fW2; 576 605 577 fCosTheta = fNuEnergy*fNuEnergy + pMu2 606 fCosTheta = fNuEnergy*fNuEnergy + pMu2 - pX2; 578 fCosTheta /= 2.*fNuEnergy*sqrt(pMu2); 607 fCosTheta /= 2.*fNuEnergy*sqrt(pMu2); 579 iTer++; 608 iTer++; 580 } 609 } 581 while( ( abs(fCosTheta) > 1. || fEmu < fMu 610 while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax ); 582 611 583 if( iTer >= iTerMax ) { fBreak = true; ret 612 if( iTer >= iTerMax ) { fBreak = true; return; } 584 613 585 if( abs(fCosTheta) > 1.) // vmg: due to bi 614 if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ... 586 { 615 { 587 G4cout<<"H2: fCosTheta = "<<fCosTheta<<" 616 G4cout<<"H2: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl; 588 // fCosTheta = -1. + 2.*G4UniformRand(); 617 // fCosTheta = -1. + 2.*G4UniformRand(); 589 if(fCosTheta < -1.) fCosTheta = -1.; 618 if(fCosTheta < -1.) fCosTheta = -1.; 590 if(fCosTheta > 1.) fCosTheta = 1.; 619 if(fCosTheta > 1.) fCosTheta = 1.; 591 } 620 } 592 // LVs 621 // LVs 593 622 594 G4LorentzVector lvt1 = G4LorentzVector( 0 623 G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., fM1 ); 595 G4LorentzVector lvsum = lvp1 + lvt1; 624 G4LorentzVector lvsum = lvp1 + lvt1; 596 625 597 cost = fCosTheta; 626 cost = fCosTheta; 598 sint = std::sqrt( (1.0 - cost)*(1.0 + cost 627 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) ); 599 phi = G4UniformRand()*CLHEP::twopi; 628 phi = G4UniformRand()*CLHEP::twopi; 600 eP = G4ThreeVector( sint*std::cos(phi), 629 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost ); 601 muMom = sqrt(fEmu*fEmu-fMu*fMu); 630 muMom = sqrt(fEmu*fEmu-fMu*fMu); 602 eP *= muMom; 631 eP *= muMom; 603 fLVl = G4LorentzVector( eP, fEmu ); 632 fLVl = G4LorentzVector( eP, fEmu ); 604 633 605 fLVh = lvsum - fLVl; 634 fLVh = lvsum - fLVl; 606 fLVt = G4LorentzVector( 0., 0., 0., 0. ); 635 fLVt = G4LorentzVector( 0., 0., 0., 0. ); // no recoil 607 } 636 } 608 else // Fermi motion, Q2 in nucleon rest fra 637 else // Fermi motion, Q2 in nucleon rest frame 609 { 638 { 610 G4Nucleus recoil1( A-1, Z ); 639 G4Nucleus recoil1( A-1, Z ); 611 rM = recoil1.AtomicMass(A-1,Z); 640 rM = recoil1.AtomicMass(A-1,Z); 612 do 641 do 613 { 642 { 614 // nMom = NucleonMomentumBR( targetNucle 643 // nMom = NucleonMomentumBR( targetNucleus ); // BR 615 nMom = GgSampleNM( targetNucleus ); // G 644 nMom = GgSampleNM( targetNucleus ); // Gg 616 Ex = GetEx(A-1, fProton); 645 Ex = GetEx(A-1, fProton); 617 ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nM 646 ei = tM - sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ); 618 // ei = 0.5*( tM - s2M - 2*eX ); 647 // ei = 0.5*( tM - s2M - 2*eX ); 619 648 620 nm2 = ei*ei - nMom*nMom; 649 nm2 = ei*ei - nMom*nMom; 621 iTer++; 650 iTer++; 622 } 651 } 623 while( nm2 < 0. && iTer < iTerMax ); 652 while( nm2 < 0. && iTer < iTerMax ); 624 653 625 if( iTer >= iTerMax ) { fBreak = true; ret 654 if( iTer >= iTerMax ) { fBreak = true; return; } 626 655 627 G4ThreeVector nMomDir = nMom*G4RandomDirec 656 G4ThreeVector nMomDir = nMom*G4RandomDirection(); 628 657 629 if( !f2p2h || A < 3 ) // 1p1h << 658 if( !f2p2h ) // 1p1h 630 { 659 { 631 // hM = tM - rM; 660 // hM = tM - rM; 632 661 633 fLVt = G4LorentzVector( -nMomDir, sqrt( 662 fLVt = G4LorentzVector( -nMomDir, sqrt( (rM + Ex)*(rM + Ex) + nMom*nMom ) ); // rM ); // 634 fLVh = G4LorentzVector( nMomDir, ei ); 663 fLVh = G4LorentzVector( nMomDir, ei ); // hM); // 635 } 664 } 636 else // 2p2h 665 else // 2p2h 637 { 666 { 638 G4Nucleus recoil(A-2,Z-1); 667 G4Nucleus recoil(A-2,Z-1); 639 rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMo 668 rM = recoil.AtomicMass(A-2,Z-1)+sqrt(nMom*nMom+fM1*fM1); 640 hM = tM - rM; 669 hM = tM - rM; 641 670 642 fLVt = G4LorentzVector( nMomDir, sqrt( r 671 fLVt = G4LorentzVector( nMomDir, sqrt( rM*rM+nMom*nMom ) ); 643 fLVh = G4LorentzVector(-nMomDir, sqrt( h 672 fLVh = G4LorentzVector(-nMomDir, sqrt( hM*hM+nMom*nMom ) ); 644 } 673 } 645 // G4cout<<hM<<", "; 674 // G4cout<<hM<<", "; 646 // bst = fLVh.boostVector(); << 675 bst = fLVh.boostVector(); 647 676 648 // lvp1.boost(-bst); // -> nucleon rest sy << 677 lvp1.boost(-bst); // -> nucleon rest system, where Q2 transfer is ??? 649 678 650 fNuEnergy = lvp1.e(); 679 fNuEnergy = lvp1.e(); 651 // G4double mN = fLVh.m(); // better mN = << 680 G4double mN = fLVh.m(); 652 iTer = 0; 681 iTer = 0; 653 682 654 do // no FM!?, 5.4.20 vmg << 683 do 655 { 684 { 656 fXsample = SampleXkr(fNuEnergy); 685 fXsample = SampleXkr(fNuEnergy); 657 fQtransfer = SampleQkr(fNuEnergy, fXsamp 686 fQtransfer = SampleQkr(fNuEnergy, fXsample); 658 fQ2 = fQtransfer*fQtransfer; 687 fQ2 = fQtransfer*fQtransfer; 659 688 660 // G4double mR = mN + fM1*(A-1.)*std::ex << 661 << 662 if( fXsample > 0. ) 689 if( fXsample > 0. ) 663 { 690 { 664 fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // << 691 // fW2 = fM1*fM1 - fQ2 + fQ2/fXsample; // sample excited hadron mass 665 << 692 fW2 = mN*mN - fQ2 + fQ2/fXsample; // sample excited hadron mass 666 // fW2 = mN*mN - fQ2 + fQ2/fXsample; / << 693 fEmu = fNuEnergy - fQ2/2./fM1/fXsample; 667 // fEmu = fNuEnergy - fQ2/2./mR/fXsamp << 668 << 669 fEmu = fNuEnergy - fQ2/2./fM1/fXsample << 670 } 694 } 671 else 695 else 672 { 696 { 673 // fW2 = mN*mN; << 697 fW2 = fM1*fM1; 674 << 675 fW2 = fM1*fM1; << 676 fEmu = fNuEnergy; 698 fEmu = fNuEnergy; 677 } 699 } >> 700 678 // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu 701 // if(fEmu < 0.) G4cout<<"fEmu = "<<fEmu<<" hM = "<<hM<<G4endl; 679 // e3 = fNuEnergy + mR - fEmu; << 680 702 681 e3 = fNuEnergy + fM1 - fEmu; 703 e3 = fNuEnergy + fM1 - fEmu; 682 704 683 // if( e3 < sqrt(fW2) ) G4cout<<"energy 705 // if( e3 < sqrt(fW2) ) G4cout<<"energyX = "<<e3/GeV<<", fW = "<<sqrt(fW2)/GeV<<G4endl; 684 706 685 pMu2 = fEmu*fEmu - fMu*fMu; 707 pMu2 = fEmu*fEmu - fMu*fMu; 686 pX2 = e3*e3 - fW2; 708 pX2 = e3*e3 - fW2; 687 709 688 if(pMu2 < 0.) { fBreak = true; return; } << 710 if(pMu2 < 0.) 689 << 711 { >> 712 fBreak = true; >> 713 return; >> 714 } 690 fCosTheta = fNuEnergy*fNuEnergy + pMu2 715 fCosTheta = fNuEnergy*fNuEnergy + pMu2 - pX2; 691 fCosTheta /= 2.*fNuEnergy*sqrt(pMu2); 716 fCosTheta /= 2.*fNuEnergy*sqrt(pMu2); 692 iTer++; 717 iTer++; 693 } 718 } 694 while( ( abs(fCosTheta) > 1. || fEmu < fMu 719 while( ( abs(fCosTheta) > 1. || fEmu < fMu ) && iTer < iTerMax ); 695 720 696 if( iTer >= iTerMax ) { fBreak = true; ret 721 if( iTer >= iTerMax ) { fBreak = true; return; } 697 722 698 if( abs(fCosTheta) > 1.) // vmg: due to bi 723 if( abs(fCosTheta) > 1.) // vmg: due to big Q2/x values. To be improved ... 699 { 724 { 700 G4cout<<"FM: fCosTheta = "<<fCosTheta<<" 725 G4cout<<"FM: fCosTheta = "<<fCosTheta<<", fEmu = "<<fEmu<<G4endl; 701 // fCosTheta = -1. + 2.*G4UniformRand(); 726 // fCosTheta = -1. + 2.*G4UniformRand(); 702 if( fCosTheta < -1.) fCosTheta = -1.; << 727 if(fCosTheta < -1.) fCosTheta = -1.; 703 if( fCosTheta > 1.) fCosTheta = 1.; << 728 if(fCosTheta > 1.) fCosTheta = 1.; 704 } 729 } 705 // LVs 730 // LVs 706 // G4LorentzVector lvt1 = G4LorentzVector << 731 G4LorentzVector lvt1 = G4LorentzVector( 0., 0., 0., mN ); // fM1 ); 707 << 708 G4LorentzVector lvt1 = G4LorentzVector( 0 << 709 G4LorentzVector lvsum = lvp1 + lvt1; 732 G4LorentzVector lvsum = lvp1 + lvt1; 710 733 711 cost = fCosTheta; 734 cost = fCosTheta; 712 sint = std::sqrt( (1.0 - cost)*(1.0 + cost 735 sint = std::sqrt( (1.0 - cost)*(1.0 + cost) ); 713 phi = G4UniformRand()*CLHEP::twopi; 736 phi = G4UniformRand()*CLHEP::twopi; 714 eP = G4ThreeVector( sint*std::cos(phi), 737 eP = G4ThreeVector( sint*std::cos(phi), sint*std::sin(phi), cost ); 715 muMom = sqrt(fEmu*fEmu-fMu*fMu); 738 muMom = sqrt(fEmu*fEmu-fMu*fMu); 716 eP *= muMom; 739 eP *= muMom; 717 fLVl = G4LorentzVector( eP, fEmu ); 740 fLVl = G4LorentzVector( eP, fEmu ); 718 fLVh = lvsum - fLVl; 741 fLVh = lvsum - fLVl; >> 742 // back to lab system >> 743 fLVl.boost(bst); >> 744 fLVh.boost(bst); >> 745 } >> 746 //G4cout<<iTer<<", "<<fBreak<<"; "; >> 747 } 719 748 720 // if( fLVh.e() < mN || fLVh.m2() < 0.) { << 749 ////////////////////////////////////// 721 750 722 if( fLVh.e() < fM1 || fLVh.m2() < 0.) { fB << 751 G4double G4NuMuNucleusCcModel::SampleXkr(G4double energy) >> 752 { >> 753 G4int i(0), nBin(50); >> 754 G4double xx(0.), prob = G4UniformRand(); 723 755 724 // back to lab system << 756 for( i = 0; i < nBin; ++i ) >> 757 { >> 758 if( energy <= fNuMuEnergyLogVector[i] ) break; >> 759 } >> 760 if( i <= 0) // E-edge >> 761 { >> 762 fEindex = 0; >> 763 xx = GetXkr( 0, prob); >> 764 } >> 765 else if ( i >= nBin) >> 766 { >> 767 fEindex = nBin-1; >> 768 xx = GetXkr( nBin-1, prob); >> 769 } >> 770 else >> 771 { >> 772 fEindex = i; >> 773 G4double x1 = GetXkr(i-1,prob); >> 774 G4double x2 = GetXkr(i,prob); >> 775 >> 776 G4double e1 = G4Log(fNuMuEnergyLogVector[i-1]); >> 777 G4double e2 = G4Log(fNuMuEnergyLogVector[i]); >> 778 G4double e = G4Log(energy); 725 779 726 // fLVl.boost(bst); << 780 if( e2 <= e1) xx = x1 + G4UniformRand()*(x2-x1); 727 // fLVh.boost(bst); << 781 else xx = x1 + (e-e1)*(x2-x1)/(e2-e1); // lin in energy log-scale 728 } 782 } 729 //G4cout<<iTer<<", "<<fBreak<<"; "; << 783 return xx; >> 784 } >> 785 >> 786 ////////////////////////////////////////////// >> 787 // >> 788 // sample X according to prob (xmin,1) at a given energy index iEnergy >> 789 >> 790 G4double G4NuMuNucleusCcModel::GetXkr(G4int iEnergy, G4double prob) >> 791 { >> 792 G4int i(0), nBin=50; >> 793 G4double xx(0.); >> 794 >> 795 for( i = 0; i < nBin; ++i ) >> 796 { >> 797 if( prob <= fNuMuXdistrKR[iEnergy][i] ) >> 798 break; >> 799 } >> 800 if(i <= 0 ) // X-edge >> 801 { >> 802 fXindex = 0; >> 803 xx = fNuMuXarrayKR[iEnergy][0]; >> 804 } >> 805 if ( i >= nBin ) >> 806 { >> 807 fXindex = nBin; >> 808 xx = fNuMuXarrayKR[iEnergy][nBin]; >> 809 } >> 810 else >> 811 { >> 812 fXindex = i; >> 813 G4double x1 = fNuMuXarrayKR[iEnergy][i]; >> 814 G4double x2 = fNuMuXarrayKR[iEnergy][i+1]; >> 815 >> 816 G4double p1 = 0.; >> 817 >> 818 if( i > 0 ) p1 = fNuMuXdistrKR[iEnergy][i-1]; >> 819 >> 820 G4double p2 = fNuMuXdistrKR[iEnergy][i]; >> 821 >> 822 if( p2 <= p1 ) xx = x1 + G4UniformRand()*(x2-x1); >> 823 else xx = x1 + (prob-p1)*(x2-x1)/(p2-p1); >> 824 } >> 825 return xx; >> 826 } >> 827 >> 828 ////////////////////////////////////// >> 829 // >> 830 // Sample fQtransfer at a given Enu and fX >> 831 >> 832 G4double G4NuMuNucleusCcModel::SampleQkr( G4double energy, G4double xx) >> 833 { >> 834 G4int nBin(50), iE=fEindex, jX=fXindex; >> 835 G4double qq(0.), qq1(0.), qq2(0.); >> 836 G4double prob = G4UniformRand(); >> 837 >> 838 // first E >> 839 >> 840 if( iE <= 0 ) >> 841 { >> 842 qq1 = GetQkr( 0, jX, prob); >> 843 } >> 844 else if ( iE >= nBin-1) >> 845 { >> 846 qq1 = GetQkr( nBin-1, jX, prob); >> 847 } >> 848 else >> 849 { >> 850 G4double q1 = GetQkr(iE-1,jX, prob); >> 851 G4double q2 = GetQkr(iE,jX, prob); >> 852 >> 853 G4double e1 = G4Log(fNuMuEnergyLogVector[iE-1]); >> 854 G4double e2 = G4Log(fNuMuEnergyLogVector[iE]); >> 855 G4double e = G4Log(energy); >> 856 >> 857 if( e2 <= e1) qq1 = q1 + G4UniformRand()*(q2-q1); >> 858 else qq1 = q1 + (e-e1)*(q2-q1)/(e2-e1); // lin in energy log-scale >> 859 } >> 860 >> 861 // then X >> 862 >> 863 if( jX <= 0 ) >> 864 { >> 865 qq2 = GetQkr( iE, 0, prob); >> 866 } >> 867 else if ( jX >= nBin) >> 868 { >> 869 qq2 = GetQkr( iE, nBin, prob); >> 870 } >> 871 else >> 872 { >> 873 G4double q1 = GetQkr(iE,jX-1, prob); >> 874 G4double q2 = GetQkr(iE,jX, prob); >> 875 >> 876 G4double e1 = G4Log(fNuMuXarrayKR[iE][jX-1]); >> 877 G4double e2 = G4Log(fNuMuXarrayKR[iE][jX]); >> 878 G4double e = G4Log(xx); >> 879 >> 880 if( e2 <= e1) qq2 = q1 + G4UniformRand()*(q2-q1); >> 881 else qq2 = q1 + (e-e1)*(q2-q1)/(e2-e1); // lin in energy log-scale >> 882 } >> 883 qq = 0.5*(qq1+qq2); >> 884 >> 885 return qq; >> 886 } >> 887 >> 888 ////////////////////////////////////////////// >> 889 // >> 890 // sample Q according to prob (qmin,qmax) at a given energy index iE and X index jX >> 891 >> 892 G4double G4NuMuNucleusCcModel::GetQkr( G4int iE, G4int jX, G4double prob ) >> 893 { >> 894 G4int i(0), nBin=50; >> 895 G4double qq(0.); >> 896 >> 897 for( i = 0; i < nBin; ++i ) >> 898 { >> 899 if( prob <= fNuMuQdistrKR[iE][jX][i] ) >> 900 break; >> 901 } >> 902 if(i <= 0 ) // Q-edge >> 903 { >> 904 fQindex = 0; >> 905 qq = fNuMuQarrayKR[iE][jX][0]; >> 906 } >> 907 if ( i >= nBin ) >> 908 { >> 909 fQindex = nBin; >> 910 qq = fNuMuQarrayKR[iE][jX][nBin]; >> 911 } >> 912 else >> 913 { >> 914 fQindex = i; >> 915 G4double q1 = fNuMuQarrayKR[iE][jX][i]; >> 916 G4double q2 = fNuMuQarrayKR[iE][jX][i+1]; >> 917 >> 918 G4double p1 = 0.; >> 919 >> 920 if( i > 0 ) p1 = fNuMuQdistrKR[iE][jX][i-1]; >> 921 >> 922 G4double p2 = fNuMuQdistrKR[iE][jX][i]; >> 923 >> 924 if( p2 <= p1 ) qq = q1 + G4UniformRand()*(q2-q1); >> 925 else qq = q1 + (prob-p1)*(q2-q1)/(p2-p1); >> 926 } >> 927 return qq; 730 } 928 } 731 929 732 // 930 // 733 // 931 // 734 /////////////////////////// 932 /////////////////////////// 735 933