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: G4hhElastic.cc,v 1.5 2010-11-09 09:04:29 grichine Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 // 29 // 28 // Physics model class G4hhElastic 30 // Physics model class G4hhElastic 29 // 31 // 30 // 32 // 31 // G4 Model: qQ hadron hadron elastic scatteri 33 // G4 Model: qQ hadron hadron elastic scattering with 4-momentum balance 32 // 34 // 33 // 02.05.2014 V. Grichine 1-st version 35 // 02.05.2014 V. Grichine 1-st version 34 // 36 // 35 37 36 #include "G4hhElastic.hh" 38 #include "G4hhElastic.hh" 37 #include "G4ParticleTable.hh" 39 #include "G4ParticleTable.hh" 38 #include "G4ParticleDefinition.hh" 40 #include "G4ParticleDefinition.hh" 39 #include "G4IonTable.hh" 41 #include "G4IonTable.hh" 40 #include "G4NucleiProperties.hh" 42 #include "G4NucleiProperties.hh" 41 43 42 #include "Randomize.hh" 44 #include "Randomize.hh" 43 #include "G4Integrator.hh" 45 #include "G4Integrator.hh" 44 #include "globals.hh" 46 #include "globals.hh" 45 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh" 46 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh" 47 49 48 #include "G4Proton.hh" 50 #include "G4Proton.hh" 49 #include "G4Neutron.hh" 51 #include "G4Neutron.hh" 50 #include "G4PionPlus.hh" 52 #include "G4PionPlus.hh" 51 #include "G4PionMinus.hh" 53 #include "G4PionMinus.hh" 52 54 53 #include "G4Element.hh" 55 #include "G4Element.hh" 54 #include "G4ElementTable.hh" 56 #include "G4ElementTable.hh" 55 #include "G4PhysicsTable.hh" 57 #include "G4PhysicsTable.hh" 56 #include "G4PhysicsLogVector.hh" 58 #include "G4PhysicsLogVector.hh" 57 #include "G4PhysicsFreeVector.hh" 59 #include "G4PhysicsFreeVector.hh" 58 60 59 #include "G4HadronNucleonXsc.hh" 61 #include "G4HadronNucleonXsc.hh" 60 62 61 #include "G4Pow.hh" 63 #include "G4Pow.hh" 62 64 63 #include "G4HadronicParameters.hh" << 64 << 65 using namespace std; 65 using namespace std; 66 66 67 67 68 ////////////////////////////////////////////// 68 ///////////////////////////////////////////////////////////////////////// 69 // 69 // 70 // Tracking constructor. Target is proton 70 // Tracking constructor. Target is proton 71 71 72 72 73 G4hhElastic::G4hhElastic() 73 G4hhElastic::G4hhElastic() 74 : G4HadronElastic("HadrHadrElastic") 74 : G4HadronElastic("HadrHadrElastic") 75 { 75 { 76 SetMinEnergy( 1.*GeV ); 76 SetMinEnergy( 1.*GeV ); 77 SetMaxEnergy( G4HadronicParameters::Instance << 77 SetMaxEnergy( 10000.*TeV ); 78 verboseLevel = 0; 78 verboseLevel = 0; 79 lowEnergyRecoilLimit = 100.*keV; 79 lowEnergyRecoilLimit = 100.*keV; 80 lowEnergyLimitQ = 0.0*GeV; 80 lowEnergyLimitQ = 0.0*GeV; 81 lowEnergyLimitHE = 0.0*GeV; 81 lowEnergyLimitHE = 0.0*GeV; 82 lowestEnergyLimit= 0.0*keV; 82 lowestEnergyLimit= 0.0*keV; 83 plabLowLimit = 20.0*MeV; 83 plabLowLimit = 20.0*MeV; 84 84 85 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 85 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 86 fInTkin=0; 86 fInTkin=0; 87 theProton = G4Proton::Proton(); 87 theProton = G4Proton::Proton(); 88 theNeutron = G4Neutron::Neutron(); 88 theNeutron = G4Neutron::Neutron(); 89 thePionPlus = G4PionPlus::PionPlus(); 89 thePionPlus = G4PionPlus::PionPlus(); 90 thePionMinus= G4PionMinus::PionMinus(); 90 thePionMinus= G4PionMinus::PionMinus(); 91 91 92 fTarget = G4Proton::Proton(); 92 fTarget = G4Proton::Proton(); 93 fProjectile = 0; 93 fProjectile = 0; 94 fHadrNuclXsc = new G4HadronNucleonXsc(); 94 fHadrNuclXsc = new G4HadronNucleonXsc(); 95 95 96 fEnergyBin = 200; 96 fEnergyBin = 200; 97 fBinT = 514; // 514; // 500; // 200; 97 fBinT = 514; // 514; // 500; // 200; 98 98 99 fEnergyVector = new G4PhysicsLogVector( the 99 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 100 100 101 fTableT = 0; 101 fTableT = 0; 102 fOldTkin = 0.; 102 fOldTkin = 0.; 103 SetParameters(); 103 SetParameters(); 104 104 105 Initialise(); 105 Initialise(); 106 } 106 } 107 107 108 108 109 ////////////////////////////////////////////// 109 ///////////////////////////////////////////////////////////////////////// 110 // 110 // 111 // test constructor 111 // test constructor 112 112 113 113 114 G4hhElastic::G4hhElastic( G4ParticleDefinition 114 G4hhElastic::G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile, G4double plab) 115 : G4HadronElastic("HadrHadrElastic") 115 : G4HadronElastic("HadrHadrElastic") 116 { 116 { 117 SetMinEnergy( 1.*GeV ); 117 SetMinEnergy( 1.*GeV ); 118 SetMaxEnergy( G4HadronicParameters::Instance << 118 SetMaxEnergy( 10000.*TeV ); 119 verboseLevel = 0; 119 verboseLevel = 0; 120 lowEnergyRecoilLimit = 100.*keV; 120 lowEnergyRecoilLimit = 100.*keV; 121 lowEnergyLimitQ = 0.0*GeV; 121 lowEnergyLimitQ = 0.0*GeV; 122 lowEnergyLimitHE = 0.0*GeV; 122 lowEnergyLimitHE = 0.0*GeV; 123 lowestEnergyLimit = 0.0*keV; 123 lowestEnergyLimit = 0.0*keV; 124 plabLowLimit = 20.0*MeV; 124 plabLowLimit = 20.0*MeV; 125 125 126 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 126 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 127 fInTkin=0; 127 fInTkin=0; 128 theProton = G4Proton::Proton(); 128 theProton = G4Proton::Proton(); 129 theNeutron = G4Neutron::Neutron(); 129 theNeutron = G4Neutron::Neutron(); 130 thePionPlus = G4PionPlus::PionPlus(); 130 thePionPlus = G4PionPlus::PionPlus(); 131 thePionMinus= G4PionMinus::PionMinus(); 131 thePionMinus= G4PionMinus::PionMinus(); 132 132 133 fTarget = target; 133 fTarget = target; 134 fProjectile = projectile; 134 fProjectile = projectile; 135 fMassTarg = fTarget->GetPDGMass(); 135 fMassTarg = fTarget->GetPDGMass(); 136 fMassProj = fProjectile->GetPDGMass(); 136 fMassProj = fProjectile->GetPDGMass(); 137 fMassSum2 = (fMassTarg+fMassProj)*(fMassTa 137 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj); 138 fMassDif2 = (fMassTarg-fMassProj)*(fMassTa 138 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj); 139 fHadrNuclXsc = new G4HadronNucleonXsc(); 139 fHadrNuclXsc = new G4HadronNucleonXsc(); 140 140 141 fEnergyBin = 200; 141 fEnergyBin = 200; 142 fBinT = 514; // 200; 142 fBinT = 514; // 200; 143 143 144 fEnergyVector = new G4PhysicsLogVector( the 144 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 145 fTableT = 0; 145 fTableT = 0; 146 fOldTkin = 0.; 146 fOldTkin = 0.; 147 147 148 148 149 SetParameters(); 149 SetParameters(); 150 SetParametersCMS( plab); 150 SetParametersCMS( plab); 151 } 151 } 152 152 153 153 154 ////////////////////////////////////////////// 154 ///////////////////////////////////////////////////////////////////////// 155 // 155 // 156 // constructor used for low mass diffraction 156 // constructor used for low mass diffraction 157 157 158 158 159 G4hhElastic::G4hhElastic( G4ParticleDefinition 159 G4hhElastic::G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile) 160 : G4HadronElastic("HadrHadrElastic") 160 : G4HadronElastic("HadrHadrElastic") 161 { 161 { 162 SetMinEnergy( 1.*GeV ); 162 SetMinEnergy( 1.*GeV ); 163 SetMaxEnergy( G4HadronicParameters::Instance << 163 SetMaxEnergy( 10000.*TeV ); 164 verboseLevel = 0; 164 verboseLevel = 0; 165 lowEnergyRecoilLimit = 100.*keV; 165 lowEnergyRecoilLimit = 100.*keV; 166 lowEnergyLimitQ = 0.0*GeV; 166 lowEnergyLimitQ = 0.0*GeV; 167 lowEnergyLimitHE = 0.0*GeV; 167 lowEnergyLimitHE = 0.0*GeV; 168 lowestEnergyLimit= 0.0*keV; 168 lowestEnergyLimit= 0.0*keV; 169 plabLowLimit = 20.0*MeV; 169 plabLowLimit = 20.0*MeV; 170 170 171 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 171 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 172 fInTkin=0; 172 fInTkin=0; 173 173 174 fTarget = target; // later vmg 174 fTarget = target; // later vmg 175 fProjectile = projectile; 175 fProjectile = projectile; 176 theProton = G4Proton::Proton(); 176 theProton = G4Proton::Proton(); 177 theNeutron = G4Neutron::Neutron(); 177 theNeutron = G4Neutron::Neutron(); 178 thePionPlus = G4PionPlus::PionPlus(); 178 thePionPlus = G4PionPlus::PionPlus(); 179 thePionMinus= G4PionMinus::PionMinus(); 179 thePionMinus= G4PionMinus::PionMinus(); 180 180 181 fTarget = G4Proton::Proton(); // later vmg 181 fTarget = G4Proton::Proton(); // later vmg >> 182 fProjectile = 0; 182 fMassTarg = fTarget->GetPDGMass(); 183 fMassTarg = fTarget->GetPDGMass(); 183 fMassProj = fProjectile->GetPDGMass(); 184 fMassProj = fProjectile->GetPDGMass(); 184 fMassSum2 = (fMassTarg+fMassProj)*(fMassTa 185 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj); 185 fMassDif2 = (fMassTarg-fMassProj)*(fMassTa 186 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj); 186 fHadrNuclXsc = new G4HadronNucleonXsc(); 187 fHadrNuclXsc = new G4HadronNucleonXsc(); 187 188 188 fEnergyBin = 200; 189 fEnergyBin = 200; 189 fBinT = 514; // 514; // 500; // 200; 190 fBinT = 514; // 514; // 500; // 200; 190 191 191 fEnergyVector = new G4PhysicsLogVector( the 192 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 192 193 193 fTableT = 0; 194 fTableT = 0; 194 fOldTkin = 0.; 195 fOldTkin = 0.; 195 196 196 SetParameters(); 197 SetParameters(); 197 } 198 } 198 199 199 200 200 201 201 ////////////////////////////////////////////// 202 ////////////////////////////////////////////////////////////////////////////// 202 // 203 // 203 // Destructor 204 // Destructor 204 205 205 G4hhElastic::~G4hhElastic() 206 G4hhElastic::~G4hhElastic() 206 { 207 { 207 if ( fEnergyVector ) { 208 if ( fEnergyVector ) { 208 delete fEnergyVector; 209 delete fEnergyVector; 209 fEnergyVector = 0; 210 fEnergyVector = 0; 210 } 211 } 211 212 212 for ( std::vector<G4PhysicsTable*>::iterator 213 for ( std::vector<G4PhysicsTable*>::iterator it = fBankT.begin(); 213 it != fBankT.end(); ++it ) { 214 it != fBankT.end(); ++it ) { 214 if ( (*it) ) (*it)->clearAndDestroy(); 215 if ( (*it) ) (*it)->clearAndDestroy(); 215 delete *it; 216 delete *it; 216 *it = 0; 217 *it = 0; 217 } 218 } 218 fTableT = 0; 219 fTableT = 0; 219 if(fHadrNuclXsc) delete fHadrNuclXsc; 220 if(fHadrNuclXsc) delete fHadrNuclXsc; 220 } 221 } 221 222 222 ////////////////////////////////////////////// 223 ///////////////////////////////////////////////////////////////////////////// 223 ///////////////////// Table preparation and r 224 ///////////////////// Table preparation and reading //////////////////////// 224 225 225 226 226 ////////////////////////////////////////////// 227 ////////////////////////////////////////////////////////////////////////////// 227 // 228 // 228 // Initialisation for given particle on the pr 229 // Initialisation for given particle on the proton target 229 230 230 void G4hhElastic::Initialise() 231 void G4hhElastic::Initialise() 231 { 232 { 232 // pp,pn 233 // pp,pn 233 234 234 fProjectile = G4Proton::Proton(); 235 fProjectile = G4Proton::Proton(); 235 BuildTableT(fTarget, fProjectile); 236 BuildTableT(fTarget, fProjectile); 236 fBankT.push_back(fTableT); // 0 237 fBankT.push_back(fTableT); // 0 237 238 238 // pi+-p 239 // pi+-p 239 240 240 fProjectile = G4PionPlus::PionPlus(); 241 fProjectile = G4PionPlus::PionPlus(); 241 BuildTableT(fTarget, fProjectile); 242 BuildTableT(fTarget, fProjectile); 242 fBankT.push_back(fTableT); // 1 243 fBankT.push_back(fTableT); // 1 243 //K+-p 244 //K+-p 244 fProjectile = G4KaonPlus::KaonPlus(); 245 fProjectile = G4KaonPlus::KaonPlus(); 245 BuildTableT(fTarget, fProjectile); 246 BuildTableT(fTarget, fProjectile); 246 fBankT.push_back(fTableT); // 2 247 fBankT.push_back(fTableT); // 2 247 248 248 } 249 } 249 250 250 ////////////////////////////////////////////// 251 /////////////////////////////////////////////////////////////////////////////// 251 // 252 // 252 // Build for given particle and proton table o 253 // Build for given particle and proton table of momentum transfers. 253 254 254 void G4hhElastic::BuildTableT( G4ParticleDefin 255 void G4hhElastic::BuildTableT( G4ParticleDefinition* target, G4ParticleDefinition* projectile) // , G4double plab) 255 { 256 { 256 G4int iTkin, jTransfer; 257 G4int iTkin, jTransfer; 257 G4double plab, Tkin, tMax; 258 G4double plab, Tkin, tMax; 258 G4double t1, t2, dt, delta = 0., sum = 0.; 259 G4double t1, t2, dt, delta = 0., sum = 0.; 259 260 260 fTarget = target; 261 fTarget = target; 261 fProjectile = projectile; 262 fProjectile = projectile; 262 fMassTarg = fTarget->GetPDGMass(); 263 fMassTarg = fTarget->GetPDGMass(); 263 fMassProj = fProjectile->GetPDGMass(); 264 fMassProj = fProjectile->GetPDGMass(); 264 fMassSum2 = (fMassTarg+fMassProj)*(fMassTa 265 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj); 265 fMassDif2 = (fMassTarg-fMassProj)*(fMassTa 266 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj); 266 267 267 G4Integrator<G4hhElastic,G4double(G4hhElasti 268 G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral; 268 // G4HadronNucleonXsc* hnXsc = new 269 // G4HadronNucleonXsc* hnXsc = new G4HadronNucleonXsc(); 269 fTableT = new G4PhysicsTable(fEnergyBin); 270 fTableT = new G4PhysicsTable(fEnergyBin); 270 271 271 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 272 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 272 { 273 { 273 Tkin = fEnergyVector->GetLowEdgeEnergy(iT 274 Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin); 274 plab = std::sqrt( Tkin*( Tkin + 2*fMassPr 275 plab = std::sqrt( Tkin*( Tkin + 2*fMassProj ) ); 275 // G4DynamicParticle* theDynamicParticle 276 // G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(projectile, 276 // 277 // G4ParticleMomentum(0.,0.,1.), 277 // 278 // Tkin); 278 // fSigmaTot = fHadrNuclXsc->GetHadronNucl 279 // fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, target ); 279 280 280 SetParametersCMS( plab ); 281 SetParametersCMS( plab ); 281 282 282 tMax = 4.*fPcms*fPcms; 283 tMax = 4.*fPcms*fPcms; 283 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*Ge 284 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ??? 284 285 285 G4PhysicsFreeVector* vectorT = new G4Physi 286 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1); 286 sum = 0.; 287 sum = 0.; 287 dt = tMax/fBinT; 288 dt = tMax/fBinT; 288 289 289 // for(j = 1; j < fBinT; j++) 290 // for(j = 1; j < fBinT; j++) 290 291 291 for( jTransfer = fBinT-1; jTransfer >= 1; 292 for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer--) 292 { 293 { 293 t1 = dt*(jTransfer-1); 294 t1 = dt*(jTransfer-1); 294 t2 = t1 + dt; 295 t2 = t1 + dt; 295 296 296 if( fMassProj > 900.*MeV ) // pp, pn 297 if( fMassProj > 900.*MeV ) // pp, pn 297 { 298 { 298 delta = integral.Legendre10(this, &G4h 299 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2); 299 // delta = integral.Legendre96(this, & 300 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2); 300 } 301 } 301 else // pi+-p, K+-p 302 else // pi+-p, K+-p 302 { 303 { 303 delta = integral.Legendre10(this, &G4h 304 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2); 304 // delta = integral.Legendre96(this, & 305 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2); 305 } 306 } 306 sum += delta; 307 sum += delta; 307 vectorT->PutValue( jTransfer-1, t1, sum 308 vectorT->PutValue( jTransfer-1, t1, sum ); // t2 308 } 309 } 309 // vectorT->PutValue( fBinT-1, dt*(fBinT-1 310 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2 310 fTableT->insertAt( iTkin, vectorT ); 311 fTableT->insertAt( iTkin, vectorT ); 311 // delete theDynamicParticle; 312 // delete theDynamicParticle; 312 } 313 } 313 // delete hnXsc; 314 // delete hnXsc; 314 315 315 return; 316 return; 316 } 317 } 317 318 318 ////////////////////////////////////////////// 319 //////////////////////////////////////////////////////////////////////////// 319 // 320 // 320 // Return inv momentum transfer -t > 0 from in 321 // Return inv momentum transfer -t > 0 from initialisation table 321 322 322 G4double G4hhElastic::SampleInvariantT( const 323 G4double G4hhElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p, 323 324 G4int, G4int ) 324 { 325 { 325 G4int iTkin, iTransfer; 326 G4int iTkin, iTransfer; 326 G4double t, t2, position, m1 = aParticle->Ge 327 G4double t, t2, position, m1 = aParticle->GetPDGMass(); 327 G4double Tkin = std::sqrt(m1*m1+p*p) - m1; 328 G4double Tkin = std::sqrt(m1*m1+p*p) - m1; 328 329 329 if( aParticle == G4Proton::Proton() || aPart 330 if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() ) 330 { 331 { 331 fTableT = fBankT[0]; 332 fTableT = fBankT[0]; 332 } 333 } 333 if( aParticle == G4PionPlus::PionPlus() || a 334 if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() ) 334 { 335 { 335 fTableT = fBankT[1]; 336 fTableT = fBankT[1]; 336 } 337 } 337 if( aParticle == G4KaonPlus::KaonPlus() || a 338 if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() ) 338 { 339 { 339 fTableT = fBankT[2]; 340 fTableT = fBankT[2]; 340 } 341 } 341 342 342 G4double delta = std::abs(Tkin - fOldTkin 343 G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin); 343 G4double deltaMax = 1.e-2; 344 G4double deltaMax = 1.e-2; 344 345 345 if ( delta < deltaMax ) iTkin = fInTkin; 346 if ( delta < deltaMax ) iTkin = fInTkin; 346 else 347 else 347 { 348 { 348 for( iTkin = 0; iTkin < fEnergyBin; iTkin+ 349 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 349 { 350 { 350 if( Tkin < fEnergyVector->GetLowEdgeEner 351 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break; 351 } 352 } 352 } 353 } 353 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi 354 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy 354 if ( iTkin < 0 ) iTkin = 0; // aga 355 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy 355 356 356 fOldTkin = Tkin; 357 fOldTkin = Tkin; 357 fInTkin = iTkin; 358 fInTkin = iTkin; 358 359 359 if (iTkin == fEnergyBin -1 || iTkin == 0 ) 360 if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges 360 { 361 { 361 position = (*(*fTableT)(iTkin))(0)*G4Unifo 362 position = (*(*fTableT)(iTkin))(0)*G4UniformRand(); 362 363 363 // G4cout<<"position = "<<position<<G4endl 364 // G4cout<<"position = "<<position<<G4endl; 364 365 365 for(iTransfer = 0; iTransfer < fBinT-1; iT 366 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++) 366 { 367 { 367 if( position >= (*(*fTableT)(iTkin))(iTr 368 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break; 368 } 369 } 369 if (iTransfer >= fBinT-1) iTransfer = fBin 370 if (iTransfer >= fBinT-1) iTransfer = fBinT-2; 370 371 371 // G4cout<<"iTransfer = "<<iTransfer<<G4en 372 // G4cout<<"iTransfer = "<<iTransfer<<G4endl; 372 373 373 t = GetTransfer(iTkin, iTransfer, position 374 t = GetTransfer(iTkin, iTransfer, position); 374 375 375 // G4cout<<"t = "<<t<<G4endl; 376 // G4cout<<"t = "<<t<<G4endl; 376 } 377 } 377 else // Tkin inside between energy table ed 378 else // Tkin inside between energy table edges 378 { 379 { 379 // position = (*(*fTableT)(iTkin))(fBinT-2 380 // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand(); 380 position = (*(*fTableT)(iTkin))(0)*G4Unifo 381 position = (*(*fTableT)(iTkin))(0)*G4UniformRand(); 381 382 382 // G4cout<<"position = "<<position<<G4endl 383 // G4cout<<"position = "<<position<<G4endl; 383 384 384 for(iTransfer = 0; iTransfer < fBinT-1; iT 385 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++) 385 { 386 { 386 // if( position < (*(*fTableT)(iTkin))(i 387 // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break; 387 if( position >= (*(*fTableT)(iTkin))(iTr 388 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break; 388 } 389 } 389 if (iTransfer >= fBinT-1) iTransfer = fBin 390 if (iTransfer >= fBinT-1) iTransfer = fBinT-2; 390 391 391 // G4cout<<"iTransfer = "<<iTransfer<<G4en 392 // G4cout<<"iTransfer = "<<iTransfer<<G4endl; 392 393 393 t2 = GetTransfer(iTkin, iTransfer, positi 394 t2 = GetTransfer(iTkin, iTransfer, position); 394 return t2; 395 return t2; 395 /* 396 /* 396 G4double t1, E1, E2, W, W1, W2; 397 G4double t1, E1, E2, W, W1, W2; 397 // G4cout<<"t2 = "<<t2<<G4endl; 398 // G4cout<<"t2 = "<<t2<<G4endl; 398 399 399 E2 = fEnergyVector->GetLowEdgeEnergy(iTkin 400 E2 = fEnergyVector->GetLowEdgeEnergy(iTkin); 400 401 401 // G4cout<<"E2 = "<<E2<<G4endl; 402 // G4cout<<"E2 = "<<E2<<G4endl; 402 403 403 iTkin--; 404 iTkin--; 404 405 405 // position = (*(*fTableT)(iTkin))(fBinT-2 406 // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand(); 406 407 407 // G4cout<<"position = "<<position<<G4endl 408 // G4cout<<"position = "<<position<<G4endl; 408 409 409 for(iTransfer = 0; iTransfer < fBinT-1; iT 410 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++) 410 { 411 { 411 // if( position < (*(*fTableT)(iTkin))(i 412 // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break; 412 if( position >= (*(*fTableT)(iTkin))(iTr 413 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break; 413 } 414 } 414 if (iTransfer >= fBinT-1) iTransfer = fBin 415 if (iTransfer >= fBinT-1) iTransfer = fBinT-2; 415 416 416 t1 = GetTransfer(iTkin, iTransfer, positi 417 t1 = GetTransfer(iTkin, iTransfer, position); 417 418 418 // G4cout<<"t1 = "<<t1<<G4endl; 419 // G4cout<<"t1 = "<<t1<<G4endl; 419 420 420 E1 = fEnergyVector->GetLowEdgeEnergy(iTkin 421 E1 = fEnergyVector->GetLowEdgeEnergy(iTkin); 421 422 422 // G4cout<<"E1 = "<<E1<<G4endl; 423 // G4cout<<"E1 = "<<E1<<G4endl; 423 424 424 W = 1.0/(E2 - E1); 425 W = 1.0/(E2 - E1); 425 W1 = (E2 - Tkin)*W; 426 W1 = (E2 - Tkin)*W; 426 W2 = (Tkin - E1)*W; 427 W2 = (Tkin - E1)*W; 427 428 428 t = W1*t1 + W2*t2; 429 t = W1*t1 + W2*t2; 429 */ 430 */ 430 } 431 } 431 return t; 432 return t; 432 } 433 } 433 434 434 435 435 ////////////////////////////////////////////// 436 //////////////////////////////////////////////////////////////////////////// 436 // 437 // 437 // Return inv momentum transfer -t > 0 from in 438 // Return inv momentum transfer -t > 0 from initialisation table 438 439 439 G4double G4hhElastic::SampleBisectionalT( cons 440 G4double G4hhElastic::SampleBisectionalT( const G4ParticleDefinition* aParticle, G4double p) 440 { 441 { 441 G4int iTkin, iTransfer; 442 G4int iTkin, iTransfer; 442 G4double t, position, m1 = aParticle->GetPD 443 G4double t, position, m1 = aParticle->GetPDGMass(); 443 G4double Tkin = std::sqrt(m1*m1+p*p) - m1; 444 G4double Tkin = std::sqrt(m1*m1+p*p) - m1; 444 445 445 if( aParticle == G4Proton::Proton() || aPart 446 if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() ) 446 { 447 { 447 fTableT = fBankT[0]; 448 fTableT = fBankT[0]; 448 } 449 } 449 if( aParticle == G4PionPlus::PionPlus() || a 450 if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() ) 450 { 451 { 451 fTableT = fBankT[1]; 452 fTableT = fBankT[1]; 452 } 453 } 453 if( aParticle == G4KaonPlus::KaonPlus() || a 454 if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() ) 454 { 455 { 455 fTableT = fBankT[2]; 456 fTableT = fBankT[2]; 456 } 457 } 457 G4double delta = std::abs(Tkin - fOldTkin 458 G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin); 458 G4double deltaMax = 1.e-2; 459 G4double deltaMax = 1.e-2; 459 460 460 if ( delta < deltaMax ) iTkin = fInTkin; 461 if ( delta < deltaMax ) iTkin = fInTkin; 461 else 462 else 462 { 463 { 463 for( iTkin = 0; iTkin < fEnergyBin; iTkin+ 464 for( iTkin = 0; iTkin < fEnergyBin; iTkin++ ) 464 { 465 { 465 if( Tkin < fEnergyVector->GetLowEdgeEner 466 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break; 466 } 467 } 467 } 468 } 468 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi 469 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy 469 if ( iTkin < 0 ) iTkin = 0; // aga 470 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy 470 471 471 fOldTkin = Tkin; 472 fOldTkin = Tkin; 472 fInTkin = iTkin; 473 fInTkin = iTkin; 473 474 474 if (iTkin == fEnergyBin -1 || iTkin == 0 ) 475 if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges 475 { 476 { 476 position = (*(*fTableT)(iTkin))(0)*G4Unifo 477 position = (*(*fTableT)(iTkin))(0)*G4UniformRand(); 477 478 478 for(iTransfer = 0; iTransfer < fBinT-1; iT 479 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++) 479 { 480 { 480 if( position >= (*(*fTableT)(iTkin))(iTr 481 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break; 481 } 482 } 482 if (iTransfer >= fBinT-1) iTransfer = fBin 483 if (iTransfer >= fBinT-1) iTransfer = fBinT-2; 483 484 484 t = GetTransfer(iTkin, iTransfer, position 485 t = GetTransfer(iTkin, iTransfer, position); 485 486 486 487 487 } 488 } 488 else // Tkin inside between energy table ed 489 else // Tkin inside between energy table edges 489 { 490 { 490 G4double rand = G4UniformRand(); 491 G4double rand = G4UniformRand(); 491 position = (*(*fTableT)(iTkin))(0)*rand; 492 position = (*(*fTableT)(iTkin))(0)*rand; 492 493 493 // 494 // 494 // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBi 495 // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBinT-2); 495 G4int sTransfer = 0, fTransfer = fBinT - 496 G4int sTransfer = 0, fTransfer = fBinT - 2, dTransfer = fTransfer - sTransfer; 496 G4double y2; 497 G4double y2; 497 498 498 for( iTransfer = 0; iTransfer < fBinT - 1; 499 for( iTransfer = 0; iTransfer < fBinT - 1; iTransfer++ ) 499 { 500 { 500 // dTransfer %= 2; 501 // dTransfer %= 2; 501 dTransfer /= 2; 502 dTransfer /= 2; 502 // dTransfer *= 0.5; 503 // dTransfer *= 0.5; 503 y2 = (*(*fTableT)(iTkin))( sTransfer + d 504 y2 = (*(*fTableT)(iTkin))( sTransfer + dTransfer ); 504 505 505 if( y2 > position ) sTransfer += dTransf 506 if( y2 > position ) sTransfer += dTransfer; 506 507 507 // if( dTransfer <= 1 ) break; 508 // if( dTransfer <= 1 ) break; 508 if( dTransfer < 1 ) break; 509 if( dTransfer < 1 ) break; 509 } 510 } 510 t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sT 511 t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sTransfer); // +(-0.5+rand)*(*fTableT)(iTkin)->GetLowEdgeEnergy(3); 511 } 512 } 512 return t; 513 return t; 513 } 514 } 514 515 515 516 516 ////////////////////////////////////////////// 517 /////////////////////////////////////////////////////////////////////////////// 517 // 518 // 518 // Build for given particle and proton table o 519 // Build for given particle and proton table of momentum transfers. 519 520 520 void G4hhElastic::BuildTableTest( G4ParticleDe 521 void G4hhElastic::BuildTableTest( G4ParticleDefinition* target, G4ParticleDefinition* projectile, G4double plab) 521 { 522 { 522 G4int jTransfer; 523 G4int jTransfer; 523 G4double tMax; // , sQq, sQG; 524 G4double tMax; // , sQq, sQG; 524 G4double t1, t2, dt, delta = 0., sum = 0. ; 525 G4double t1, t2, dt, delta = 0., sum = 0. ; // , threshold; 525 526 526 fTarget = target; 527 fTarget = target; 527 fProjectile = projectile; 528 fProjectile = projectile; 528 fMassTarg = fTarget->GetPDGMass(); 529 fMassTarg = fTarget->GetPDGMass(); 529 fMassProj = fProjectile->GetPDGMass(); 530 fMassProj = fProjectile->GetPDGMass(); 530 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg 531 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj); 531 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg 532 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj); 532 fSpp = fMassProj*fMassProj + fMassTarg*fMas 533 fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj); 533 fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp 534 fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp); 534 535 535 G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fM 536 G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fMassProj = "<<fMassProj<<" MeV"<<G4endl; 536 tMax = 4.*fPcms*fPcms; 537 tMax = 4.*fPcms*fPcms; 537 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; 538 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ??? 538 539 539 540 540 G4Integrator<G4hhElastic,G4double(G4hhElasti 541 G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral; 541 fTableT = new G4PhysicsTable(1); 542 fTableT = new G4PhysicsTable(1); 542 G4PhysicsFreeVector* vectorT = new G4Physics 543 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1); 543 544 544 sum = 0.; 545 sum = 0.; 545 dt = tMax/G4double(fBinT); 546 dt = tMax/G4double(fBinT); 546 G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; 547 G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; fPcms = "<<fPcms/GeV 547 <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt 548 <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt = "<<dt/GeV/GeV<<" GeV2"<<G4endl; 548 549 549 // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fR 550 // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fRB*GeV<<G4endl; 550 551 551 // for(jTransfer = 1; jTransfer < fBinT; j 552 // for(jTransfer = 1; jTransfer < fBinT; jTransfer++) 552 for( jTransfer = fBinT-1; jTransfer >= 1; jT 553 for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer-- ) 553 { 554 { 554 t1 = dt*(jTransfer-1); 555 t1 = dt*(jTransfer-1); 555 t2 = t1 + dt; 556 t2 = t1 + dt; 556 557 557 if( fMassProj > 900.*MeV ) // pp, pn 558 if( fMassProj > 900.*MeV ) // pp, pn 558 { 559 { 559 delta = integral.Legendre10(this, &G4hhE 560 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2); 560 // threshold = integral.Legendre96(this, 561 // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, tMax); 561 } 562 } 562 else // pi+-p, K+-p 563 else // pi+-p, K+-p 563 { 564 { 564 delta = integral.Legendre10(this, &G4hhE 565 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2); 565 // threshold = integral.Legendre96(this, 566 // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, tMax); 566 // delta = integral.Legendre96(this, &G4 567 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2); 567 } 568 } 568 sum += delta; 569 sum += delta; 569 // G4cout<<delta<<"\t"<<sum<<"\t"<<thresho 570 // G4cout<<delta<<"\t"<<sum<<"\t"<<threshold<<G4endl; 570 571 571 // sQq = GetdsdtF123(q1); 572 // sQq = GetdsdtF123(q1); 572 // sQG = GetdsdtF123qQgG(q1); 573 // sQG = GetdsdtF123qQgG(q1); 573 // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/milli 574 // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/millibarn<<"\t"<<sQq*GeV*GeV/millibarn<<G4endl; 574 // G4cout<<"sum = "<<sum<<", "; 575 // G4cout<<"sum = "<<sum<<", "; 575 576 576 vectorT->PutValue( jTransfer-1, t1, sum ); 577 vectorT->PutValue( jTransfer-1, t1, sum ); // t2 577 } 578 } 578 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 579 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2 579 fTableT->insertAt( 0, vectorT ); 580 fTableT->insertAt( 0, vectorT ); 580 fBankT.push_back( fTableT ); // 0 581 fBankT.push_back( fTableT ); // 0 581 582 582 // for(jTransfer = 0; jTransfer < fBinT-1; j 583 // for(jTransfer = 0; jTransfer < fBinT-1; jTransfer++) 583 // G4cout<<(*(*fTableT)(0))(jTransfer)/sum 584 // G4cout<<(*(*fTableT)(0))(jTransfer)/sum<<"\t\t"<<G4Pow::GetInstance()->powN(2.,-jTransfer)<<G4endl; 584 585 585 return; 586 return; 586 } 587 } 587 588 588 589 589 ////////////////////////////////////////////// 590 //////////////////////////////////////////////////////////////////////////// 590 // 591 // 591 // Return inv momentum transfer -t > 0 from in 592 // Return inv momentum transfer -t > 0 from initialisation table 592 593 593 G4double G4hhElastic::SampleTest(G4double tMin 594 G4double G4hhElastic::SampleTest(G4double tMin ) // const G4ParticleDefinition* aParticle, ) 594 { 595 { 595 G4int iTkin, iTransfer, iTmin; 596 G4int iTkin, iTransfer, iTmin; 596 G4double t, position; 597 G4double t, position; 597 // G4double qMin = std::sqrt(tMin); 598 // G4double qMin = std::sqrt(tMin); 598 599 599 fTableT = fBankT[0]; 600 fTableT = fBankT[0]; 600 iTkin = 0; 601 iTkin = 0; 601 602 602 for(iTransfer = 0; iTransfer < fBinT-1; iTra 603 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++) 603 { 604 { 604 // if( qMin <= (*fTableT)(iTkin)->GetLowEd 605 // if( qMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break; 605 if( tMin <= (*fTableT)(iTkin)->GetLowEdgeE 606 if( tMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break; 606 } 607 } 607 iTmin = iTransfer-1; 608 iTmin = iTransfer-1; 608 if(iTmin < 0 ) iTmin = 0; 609 if(iTmin < 0 ) iTmin = 0; 609 610 610 position = (*(*fTableT)(iTkin))(iTmin)*G4Uni 611 position = (*(*fTableT)(iTkin))(iTmin)*G4UniformRand(); 611 612 612 for( iTmin = 0; iTransfer < fBinT-1; iTransf 613 for( iTmin = 0; iTransfer < fBinT-1; iTransfer++) 613 { 614 { 614 if( position > (*(*fTableT)(iTkin))(iTrans 615 if( position > (*(*fTableT)(iTkin))(iTransfer) ) break; 615 } 616 } 616 if (iTransfer >= fBinT-1) iTransfer = fBinT- 617 if (iTransfer >= fBinT-1) iTransfer = fBinT-2; 617 618 618 t = GetTransfer(iTkin, iTransfer, position); 619 t = GetTransfer(iTkin, iTransfer, position); 619 620 620 return t; 621 return t; 621 } 622 } 622 623 623 624 624 ////////////////////////////////////////////// 625 ///////////////////////////////////////////////////////////////////////////////// 625 // 626 // 626 // Check with PAI sampling 627 // Check with PAI sampling 627 628 628 G4double 629 G4double 629 G4hhElastic:: GetTransfer( G4int iTkin, G4int 630 G4hhElastic:: GetTransfer( G4int iTkin, G4int iTransfer, G4double position ) 630 { 631 { 631 G4double x1, x2, y1, y2, randTransfer, delta 632 G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6; 632 633 633 if( iTransfer == 0 ) 634 if( iTransfer == 0 ) 634 { 635 { 635 randTransfer = (*fTableT)(iTkin)->GetLowEd 636 randTransfer = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer); 636 // iTransfer++; 637 // iTransfer++; 637 } 638 } 638 else 639 else 639 { 640 { 640 if ( iTransfer >= G4int((*fTableT)(iTkin)- 641 if ( iTransfer >= G4int((*fTableT)(iTkin)->GetVectorLength()) ) 641 { 642 { 642 iTransfer = G4int((*fTableT)(iTkin)->Get << 643 iTransfer = (*fTableT)(iTkin)->GetVectorLength() - 1; 643 } 644 } 644 y1 = (*(*fTableT)(iTkin))(iTransfer-1); 645 y1 = (*(*fTableT)(iTkin))(iTransfer-1); 645 y2 = (*(*fTableT)(iTkin))(iTransfer); 646 y2 = (*(*fTableT)(iTkin))(iTransfer); 646 647 647 x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(i 648 x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer-1); 648 x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(i 649 x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer); 649 650 650 delta = y2 - y1; 651 delta = y2 - y1; 651 mean = y2 + y1; 652 mean = y2 + y1; 652 653 653 if ( x1 == x2 ) randTransfer = x2; 654 if ( x1 == x2 ) randTransfer = x2; 654 else 655 else 655 { 656 { 656 // if ( y1 == y2 ) 657 // if ( y1 == y2 ) 657 if ( delta < epsilon*mean ) 658 if ( delta < epsilon*mean ) 658 randTransfer = x1 + ( x2 - x1 )*G4U 659 randTransfer = x1 + ( x2 - x1 )*G4UniformRand(); 659 else randTransfer = x1 + ( position - y1 660 else randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 ); 660 } 661 } 661 } 662 } 662 return randTransfer; 663 return randTransfer; 663 } 664 } 664 665 665 const G4double G4hhElastic::theNuclNuclData[19 << 666 const G4double G4hhElastic::theNuclNuclData[18][6] = 666 { 667 { 667 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1 668 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof 668 669 669 { 2.76754, 4.8, 4.8, 0.05, 0.742441, 10.5 670 { 2.76754, 4.8, 4.8, 0.05, 0.742441, 10.5 }, // pp 3GeV/c 670 { 3.07744, 5.4, 5.4, 0.02, 0.83818, 6.5 671 { 3.07744, 5.4, 5.4, 0.02, 0.83818, 6.5 }, // pp 4GeV/c 671 { 3.36305, 5.2, 5.2, 0.02, 0.838893, 7.5 672 { 3.36305, 5.2, 5.2, 0.02, 0.838893, 7.5 }, // np 5GeV/c 672 { 4.32941, 6, 6, 0.03, 0.769389, 7.5 }, 673 { 4.32941, 6, 6, 0.03, 0.769389, 7.5 }, // np 9 GeV/c 673 { 4.62126, 6, 6, 0.03, 0.770111, 6.5 }, 674 { 4.62126, 6, 6, 0.03, 0.770111, 6.5 }, // pp 10.4 GeV/c 674 675 675 { 5.47416, 4.5, 4.5, 0.03, 0.813185, 7.5 676 { 5.47416, 4.5, 4.5, 0.03, 0.813185, 7.5 }, // np 15 GeV/c 676 { 6.15088, 6.5, 6.5, 0.02, 0.799539, 6.5 677 { 6.15088, 6.5, 6.5, 0.02, 0.799539, 6.5 }, // pp 19.2 GeV/c 677 { 6.77474, 5.2, 5.2, 0.03, 0.784901, 7.5 678 { 6.77474, 5.2, 5.2, 0.03, 0.784901, 7.5 }, // np 23.5 GeV/c 678 { 9.77775, 7, 7, 0.03, 0.742531, 6.5 }, 679 { 9.77775, 7, 7, 0.03, 0.742531, 6.5 }, // pp 50 GeV/c 679 // {9.77775, 7, 7, 0.011, 680 // {9.77775, 7, 7, 0.011, 0.84419, 4.5 }, // pp 50 GeV/c 680 { 10.4728, 5.2, 5.2, 0.03, 0.780439, 7.5 681 { 10.4728, 5.2, 5.2, 0.03, 0.780439, 7.5 }, // np 57.5 GeV/c 681 682 682 { 13.7631, 7, 7, 0.008, 0.8664, 683 { 13.7631, 7, 7, 0.008, 0.8664, 5.0 }, // pp 100 GeV/c 683 { 19.4184, 6.8, 6.8, 0.009, 0.861337, 2. 684 { 19.4184, 6.8, 6.8, 0.009, 0.861337, 2.5 }, // pp 200 GeV/c 684 { 23.5, 6.8, 6.8, 0.007, 0.878112, 1.5 } 685 { 23.5, 6.8, 6.8, 0.007, 0.878112, 1.5 }, // pp 23.5 GeV 685 // {24.1362, 6.4, 6.4, 0.09 686 // {24.1362, 6.4, 6.4, 0.09, 0.576215, 7.5 }, // np 309.5 GeV/c 686 { 24.1362, 7.2, 7.2, 0.008, 0.864745, 5. 687 { 24.1362, 7.2, 7.2, 0.008, 0.864745, 5.5 }, 687 { 52.8, 6.8, 6.8, 0.008, 0.871929, 1.5 } 688 { 52.8, 6.8, 6.8, 0.008, 0.871929, 1.5 }, // pp 58.2 GeV 688 689 689 { 546, 7.4, 7.4, 0.013, 0.845877, 690 { 546, 7.4, 7.4, 0.013, 0.845877, 5.5 }, // pb-p 546 GeV 690 { 1960, 7.8, 7.8, 0.022, 0.809062, 7.5 } 691 { 1960, 7.8, 7.8, 0.022, 0.809062, 7.5 }, // pb-p 1960 GeV 691 { 7000, 8, 8, 0.024, 0.820441, 5.5 }, / << 692 { 7000, 8, 8, 0.024, 0.820441, 5.5 } // pp TOTEM 692 { 13000, 8.5, 8.5, 0.03, 0.796721, 10.5 << 693 693 694 }; 694 }; 695 695 696 ////////////////////////////////////////////// 696 ////////////////////////////////////////////////////////////////////////////////// 697 697 698 const G4double G4hhElastic::thePiKaNuclData[8] 698 const G4double G4hhElastic::thePiKaNuclData[8][6] = 699 { 699 { 700 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1 700 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof 701 701 702 { 2.5627, 3.8, 3.3, 0.22, 0.222, 702 { 2.5627, 3.8, 3.3, 0.22, 0.222, 1.5 }, // pipp 3.017 GeV/c 703 { 2.93928, 4.3, 3.8, 0.2, 0.2506 703 { 2.93928, 4.3, 3.8, 0.2, 0.250601, 1.3 }, // pipp 4.122 GeV/c 704 { 3.22326, 4.8, 4.3, 0.13, 0.32751, 2.5 704 { 3.22326, 4.8, 4.3, 0.13, 0.32751, 2.5 }, // pipp 5.055 GeV/c 705 { 7.80704, 5.5, 5, 0.13, 0.340631, 2.5 }, 705 { 7.80704, 5.5, 5, 0.13, 0.340631, 2.5 }, // pipp 32 GeV/c 706 { 9.7328, 5, 4.5, 0.05, 0.416319, 5.5 }, 706 { 9.7328, 5, 4.5, 0.05, 0.416319, 5.5 }, // pipp 50 GeV/c 707 707 708 { 13.7315, 5.3, 4.8, 0.05, 0.418426, 5.5 708 { 13.7315, 5.3, 4.8, 0.05, 0.418426, 5.5 }, // pipp 100 GeV/c 709 { 16.6359, 6.3, 5.8, 0.05, 0.423817, 5.5 709 { 16.6359, 6.3, 5.8, 0.05, 0.423817, 5.5 }, // pipp 147 GeV/c 710 { 19.3961, 5, 4.5, 0.05, 0.413477, 3.5 } 710 { 19.3961, 5, 4.5, 0.05, 0.413477, 3.5 } // pimp 200 GeV/c 711 711 712 }; 712 }; 713 713 714 // 714 // 715 // 715 // 716 ////////////////////////////////////////////// 716 ///////////////////////////////////////////////////////////////////////////////// 717 717