[819] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
[1340] | 26 | // $Id: G4MuNuclearInteraction.cc,v 1.14 2010/09/08 08:59:29 vnivanch Exp $ |
---|
| 27 | // GEANT4 tag $Name: geant4-09-03-ref-09 $ |
---|
[819] | 28 | // |
---|
| 29 | // $Id: |
---|
| 30 | // -------------------------------------------------------------- |
---|
| 31 | // GEANT 4 class implementation file |
---|
| 32 | // |
---|
| 33 | // History: first implementation, based on object model of |
---|
| 34 | // 2nd December 1995, G.Cosmo |
---|
| 35 | // -------- G4MuNuclearInteraction physics process --------- |
---|
| 36 | // by Laszlo Urban, May 1998 |
---|
| 37 | // added simple model for hadronic vertex, J.P. Wellisch, November 1998 |
---|
| 38 | // -------------------------------------------------------------- |
---|
[1055] | 39 | // 26/10/1998: new corr.s from R.Kokoulin + cleanup , L.Urban |
---|
| 40 | // 23/01/2009 V.Ivanchenko Add deregistration |
---|
[819] | 41 | // |
---|
[1055] | 42 | |
---|
[819] | 43 | #include "G4MuNuclearInteraction.hh" |
---|
| 44 | #include "G4UnitsTable.hh" |
---|
[962] | 45 | #include "G4HadronicProcessStore.hh" |
---|
[819] | 46 | |
---|
| 47 | // static members ........ |
---|
| 48 | G4int G4MuNuclearInteraction::nzdat = 5 ; |
---|
| 49 | G4double G4MuNuclearInteraction::zdat[]={1.,4.,13.,29.,92.}; |
---|
| 50 | G4double G4MuNuclearInteraction::adat[]={1.01,9.01,26.98,63.55,238.03}; |
---|
| 51 | G4int G4MuNuclearInteraction::ntdat = 8 ; |
---|
| 52 | G4double G4MuNuclearInteraction::tdat[]={1.e3,1.e4,1.e5,1.e6,1.e7,1.e8, |
---|
| 53 | 1.e9,1.e10}; |
---|
| 54 | G4int G4MuNuclearInteraction::NBIN = 1000 ; |
---|
| 55 | G4double G4MuNuclearInteraction::ya[1001]={0.}; |
---|
| 56 | G4double G4MuNuclearInteraction::proba[5][8][1001]={{{0.}}}; |
---|
| 57 | |
---|
| 58 | G4MuNuclearInteraction::G4MuNuclearInteraction(const G4String& processName) |
---|
| 59 | : G4VDiscreteProcess(processName), |
---|
| 60 | theMeanFreePathTable(0), |
---|
| 61 | theCrossSectionTable(0), |
---|
| 62 | LowestKineticEnergy (1.*GeV), |
---|
| 63 | HighestKineticEnergy (1000000.*TeV), |
---|
| 64 | TotBin(100), |
---|
| 65 | CutFixed ( 0.200*GeV), |
---|
| 66 | GramPerMole(g/mole), |
---|
| 67 | theMuonMinus ( G4MuonMinus::MuonMinus() ), |
---|
| 68 | theMuonPlus ( G4MuonPlus::MuonPlus() ), |
---|
| 69 | thePionZero (G4PionZero::PionZero() ) |
---|
[962] | 70 | { |
---|
| 71 | SetProcessSubType(fHadronInelastic); |
---|
| 72 | G4HadronicProcessStore::Instance()->RegisterExtraProcess(this); |
---|
| 73 | } |
---|
[819] | 74 | |
---|
| 75 | G4MuNuclearInteraction::~G4MuNuclearInteraction() |
---|
| 76 | { |
---|
[1055] | 77 | G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this); |
---|
| 78 | |
---|
[962] | 79 | if (theMeanFreePathTable) { |
---|
| 80 | theMeanFreePathTable->clearAndDestroy(); |
---|
| 81 | delete theMeanFreePathTable; |
---|
| 82 | } |
---|
| 83 | if (theCrossSectionTable) { |
---|
| 84 | theCrossSectionTable->clearAndDestroy(); |
---|
| 85 | delete theCrossSectionTable; |
---|
| 86 | } |
---|
| 87 | |
---|
| 88 | if (&PartialSumSigma) { |
---|
| 89 | PartialSumSigma.clearAndDestroy(); |
---|
| 90 | } |
---|
[819] | 91 | } |
---|
| 92 | |
---|
| 93 | void G4MuNuclearInteraction::SetPhysicsTableBining(G4double lowE, |
---|
[962] | 94 | G4double highE, G4int nBins) |
---|
[819] | 95 | { |
---|
| 96 | LowestKineticEnergy = lowE; HighestKineticEnergy = highE ; TotBin = nBins ; |
---|
| 97 | } |
---|
| 98 | |
---|
[1228] | 99 | |
---|
| 100 | G4bool G4MuNuclearInteraction::IsApplicable(const G4ParticleDefinition& particle) |
---|
| 101 | { |
---|
| 102 | return( (&particle == theMuonMinus)||(&particle == theMuonPlus)) ; |
---|
| 103 | } |
---|
| 104 | |
---|
[962] | 105 | void G4MuNuclearInteraction::PreparePhysicsTable( |
---|
| 106 | const G4ParticleDefinition& aParticleType) |
---|
| 107 | { |
---|
| 108 | G4HadronicProcessStore::Instance() |
---|
| 109 | ->RegisterParticleForExtraProcess(this, &aParticleType); |
---|
| 110 | } |
---|
| 111 | |
---|
[819] | 112 | void G4MuNuclearInteraction::BuildPhysicsTable( |
---|
| 113 | const G4ParticleDefinition& aParticleType) |
---|
| 114 | { |
---|
[962] | 115 | G4HadronicProcessStore::Instance()->PrintInfo(&aParticleType); |
---|
| 116 | |
---|
[819] | 117 | G4double LowEdgeEnergy , Value; |
---|
| 118 | G4PhysicsLogVector* ptrVector; |
---|
| 119 | |
---|
| 120 | if (theCrossSectionTable) { |
---|
[962] | 121 | theCrossSectionTable->clearAndDestroy() ; |
---|
| 122 | delete theCrossSectionTable ; |
---|
[819] | 123 | } |
---|
| 124 | |
---|
| 125 | // make tables for the sampling at initialization |
---|
| 126 | if (theMeanFreePathTable == 0) MakeSamplingTables(&aParticleType); |
---|
| 127 | |
---|
| 128 | theCrossSectionTable = new G4PhysicsTable (G4Element::GetNumberOfElements()); |
---|
| 129 | const G4ElementTable* theElementTable = G4Element::GetElementTable() ; |
---|
| 130 | G4double AtomicNumber,AtomicWeight ; |
---|
| 131 | |
---|
| 132 | for (size_t J=0; J < G4Element::GetNumberOfElements(); J++ ) |
---|
| 133 | { |
---|
| 134 | ptrVector = new G4PhysicsLogVector(LowestKineticEnergy, |
---|
[962] | 135 | HighestKineticEnergy,TotBin) ; |
---|
[819] | 136 | AtomicNumber = (*theElementTable )[J]->GetZ() ; |
---|
| 137 | AtomicWeight = (*theElementTable )[J]->GetA() ; |
---|
| 138 | |
---|
| 139 | for ( G4int i = 0 ; i < TotBin ; i++) |
---|
| 140 | { |
---|
| 141 | LowEdgeEnergy = ptrVector->GetLowEdgeEnergy(i) ; |
---|
| 142 | Value = ComputeMicroscopicCrossSection(&aParticleType, |
---|
[962] | 143 | LowEdgeEnergy, |
---|
| 144 | AtomicNumber,AtomicWeight) ; |
---|
[819] | 145 | ptrVector->PutValue(i,Value) ; |
---|
| 146 | } |
---|
| 147 | |
---|
| 148 | theCrossSectionTable->insertAt( J , ptrVector ) ; |
---|
| 149 | } |
---|
| 150 | |
---|
| 151 | G4double FixedEnergy = (LowestKineticEnergy + HighestKineticEnergy)/2. ; |
---|
| 152 | const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ; |
---|
| 153 | if (theMeanFreePathTable) { |
---|
[962] | 154 | theMeanFreePathTable->clearAndDestroy(); |
---|
| 155 | delete theMeanFreePathTable; |
---|
[819] | 156 | } |
---|
| 157 | theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials()); |
---|
| 158 | |
---|
| 159 | PartialSumSigma.clearAndDestroy(); |
---|
| 160 | PartialSumSigma.resize(G4Material::GetNumberOfMaterials()); |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | for (size_t K=0 ; K < G4Material::GetNumberOfMaterials(); K++ ) |
---|
| 164 | { |
---|
| 165 | ptrVector = new G4PhysicsLogVector(LowestKineticEnergy, |
---|
[962] | 166 | HighestKineticEnergy, |
---|
| 167 | TotBin ) ; |
---|
[819] | 168 | |
---|
| 169 | const G4Material* material= (*theMaterialTable)[K]; |
---|
| 170 | |
---|
| 171 | for ( G4int i = 0 ; i < TotBin ; i++ ) |
---|
| 172 | { |
---|
| 173 | LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ; |
---|
| 174 | Value = ComputeMeanFreePath( &aParticleType, LowEdgeEnergy, |
---|
| 175 | material ); |
---|
| 176 | |
---|
| 177 | ptrVector->PutValue( i , Value ) ; |
---|
| 178 | } |
---|
| 179 | |
---|
| 180 | theMeanFreePathTable->insertAt( K , ptrVector ); |
---|
| 181 | |
---|
| 182 | // Compute the PartialSumSigma table at a given fixed energy |
---|
| 183 | ComputePartialSumSigma( &aParticleType, FixedEnergy, material); |
---|
| 184 | } |
---|
| 185 | |
---|
| 186 | if (&aParticleType == theMuonPlus) PrintInfoDefinition(); |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | void G4MuNuclearInteraction::ComputePartialSumSigma( |
---|
| 190 | const G4ParticleDefinition* ParticleType, |
---|
[962] | 191 | G4double KineticEnergy, |
---|
| 192 | const G4Material* aMaterial) |
---|
[819] | 193 | |
---|
| 194 | // Build the table of cross section per element. The table is built for MATERIALS. |
---|
| 195 | // This table is used by DoIt to select randomly an element in the material. |
---|
| 196 | { |
---|
| 197 | G4int Imate = aMaterial->GetIndex(); |
---|
| 198 | G4int NbOfElements = aMaterial->GetNumberOfElements(); |
---|
| 199 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
| 200 | const G4double* theAtomNumDensityVector = aMaterial->GetAtomicNumDensityVector(); |
---|
| 201 | |
---|
| 202 | PartialSumSigma[Imate] = new G4DataVector(); |
---|
| 203 | |
---|
| 204 | G4double SIGMA = 0. ; |
---|
| 205 | |
---|
| 206 | for ( G4int Ielem=0 ; Ielem < NbOfElements ; Ielem++ ) |
---|
[962] | 207 | { |
---|
| 208 | SIGMA += theAtomNumDensityVector[Ielem] * |
---|
| 209 | ComputeMicroscopicCrossSection( ParticleType, KineticEnergy, |
---|
| 210 | (*theElementVector)[Ielem]->GetZ(), |
---|
| 211 | (*theElementVector)[Ielem]->GetA()) ; |
---|
| 212 | PartialSumSigma[Imate]->push_back(SIGMA); |
---|
| 213 | } |
---|
[819] | 214 | } |
---|
| 215 | |
---|
| 216 | G4double G4MuNuclearInteraction::ComputeMicroscopicCrossSection( |
---|
| 217 | const G4ParticleDefinition* ParticleType, |
---|
| 218 | G4double KineticEnergy, |
---|
| 219 | G4double AtomicNumber,G4double AtomicWeight) |
---|
| 220 | { |
---|
| 221 | static const G4double |
---|
| 222 | xgi[] ={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 }; |
---|
| 223 | static const G4double |
---|
| 224 | wgi[] ={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 }; |
---|
| 225 | static const G4double ak1=6.9 ; |
---|
| 226 | static const G4double ak2=1.0 ; |
---|
| 227 | |
---|
| 228 | G4double Mass,epmin,epmax,epln,ep,aaa,bbb,hhh,x ; |
---|
| 229 | G4int kkk ; |
---|
| 230 | |
---|
| 231 | Mass = ParticleType->GetPDGMass() ; |
---|
| 232 | |
---|
| 233 | G4double CrossSection = 0.0 ; |
---|
| 234 | |
---|
| 235 | if ( AtomicNumber < 1. ) return CrossSection; |
---|
| 236 | if ( KineticEnergy <= CutFixed ) return CrossSection; |
---|
| 237 | |
---|
| 238 | epmin = CutFixed ; |
---|
| 239 | epmax = KineticEnergy + Mass - 0.5*proton_mass_c2 ; |
---|
| 240 | if (epmax <= epmin ) return CrossSection; // NaN bug correction |
---|
| 241 | |
---|
| 242 | aaa = std::log(epmin) ; |
---|
| 243 | bbb = std::log(epmax) ; |
---|
| 244 | kkk = G4int((bbb-aaa)/ak1 +ak2) ; |
---|
| 245 | hhh = (bbb-aaa)/kkk ; |
---|
| 246 | |
---|
| 247 | for (G4int l=0 ; l<kkk; l++) |
---|
| 248 | { |
---|
| 249 | x = aaa + hhh*l ; |
---|
| 250 | for (G4int ll=0; ll<8; ll++) |
---|
| 251 | { |
---|
| 252 | epln=x+xgi[ll]*hhh ; |
---|
| 253 | ep = std::exp(epln) ; |
---|
| 254 | CrossSection += ep*wgi[ll]* ComputeDMicroscopicCrossSection(ParticleType, |
---|
| 255 | KineticEnergy, |
---|
| 256 | AtomicNumber,AtomicWeight, |
---|
| 257 | ep) ; |
---|
| 258 | } |
---|
| 259 | } |
---|
| 260 | CrossSection *= hhh ; |
---|
| 261 | |
---|
| 262 | if (CrossSection < 0.) CrossSection = 0.; |
---|
| 263 | |
---|
| 264 | return CrossSection; |
---|
| 265 | } |
---|
| 266 | |
---|
| 267 | G4double G4MuNuclearInteraction::ComputeDMicroscopicCrossSection( |
---|
| 268 | const G4ParticleDefinition* ParticleType, |
---|
| 269 | G4double KineticEnergy, |
---|
| 270 | G4double, G4double AtomicWeight, |
---|
| 271 | G4double epsilon) |
---|
| 272 | // Calculates the differential (D) microscopic cross section |
---|
| 273 | // using the cross section formula of R.P. Kokoulin (18/01/98) |
---|
| 274 | { |
---|
| 275 | const G4double alam2 = 0.400*GeV*GeV ; |
---|
| 276 | const G4double alam = 0.632456*GeV ; |
---|
| 277 | const G4double coeffn = fine_structure_const/pi ; |
---|
| 278 | |
---|
| 279 | G4double ep,a,aeff,sigph,v,v1,v2,mass2,up,down ; |
---|
| 280 | G4double ParticleMass = ParticleType->GetPDGMass() ; |
---|
| 281 | G4double TotalEnergy = KineticEnergy + ParticleMass ; |
---|
| 282 | |
---|
| 283 | G4double DCrossSection = 0. ; |
---|
| 284 | |
---|
| 285 | if((epsilon >= TotalEnergy - 0.5*proton_mass_c2) |
---|
| 286 | || |
---|
| 287 | (epsilon <= CutFixed)) |
---|
| 288 | return DCrossSection ; |
---|
| 289 | |
---|
| 290 | ep = epsilon/GeV ; |
---|
| 291 | a = AtomicWeight/(g/mole) ; |
---|
| 292 | |
---|
| 293 | aeff = 0.22*a+0.78*std::exp(0.89*std::log(a)) ; //shadowing |
---|
| 294 | sigph = (49.2+11.1*std::log(ep)+151.8/std::sqrt(ep))*microbarn ; //!!!!!!!!!!! |
---|
| 295 | |
---|
| 296 | v=epsilon/TotalEnergy ; |
---|
| 297 | v1 = 1.-v ; |
---|
| 298 | v2 = v*v ; |
---|
| 299 | mass2 = ParticleMass*ParticleMass ; |
---|
| 300 | |
---|
| 301 | up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1)) ; |
---|
| 302 | down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam) ; |
---|
| 303 | |
---|
| 304 | DCrossSection = coeffn*aeff*sigph/epsilon* |
---|
| 305 | (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*std::log(up/down)) ; |
---|
| 306 | |
---|
| 307 | if( DCrossSection < 0.) |
---|
| 308 | DCrossSection = 0. ; |
---|
| 309 | |
---|
| 310 | return DCrossSection ; |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | void G4MuNuclearInteraction::MakeSamplingTables( |
---|
| 314 | const G4ParticleDefinition* ParticleType) |
---|
| 315 | { |
---|
| 316 | |
---|
| 317 | G4int nbin; |
---|
| 318 | G4double AtomicNumber,AtomicWeight,KineticEnergy, |
---|
| 319 | TotalEnergy,Maxep ; |
---|
| 320 | |
---|
| 321 | G4double ParticleMass = ParticleType->GetPDGMass() ; |
---|
| 322 | |
---|
| 323 | for (G4int iz=0; iz<nzdat; iz++) |
---|
| 324 | { |
---|
| 325 | AtomicNumber = zdat[iz]; |
---|
| 326 | AtomicWeight = adat[iz]*GramPerMole ; |
---|
| 327 | |
---|
| 328 | for (G4int it=0; it<ntdat; it++) |
---|
| 329 | { |
---|
| 330 | KineticEnergy = tdat[it]; |
---|
| 331 | TotalEnergy = KineticEnergy + ParticleMass; |
---|
| 332 | Maxep = TotalEnergy - 0.5*proton_mass_c2 ; |
---|
| 333 | |
---|
| 334 | G4double CrossSection = 0.0 ; |
---|
| 335 | |
---|
| 336 | G4double c,y,ymin,ymax,dy,yy,dx,x,ep ; |
---|
| 337 | |
---|
| 338 | // calculate the differential cross section |
---|
| 339 | // numerical integration in |
---|
| 340 | // log ............... |
---|
| 341 | c = std::log(Maxep/CutFixed) ; |
---|
| 342 | ymin = -5. ; |
---|
| 343 | ymax = 0. ; |
---|
| 344 | dy = (ymax-ymin)/NBIN ; |
---|
| 345 | |
---|
| 346 | nbin=-1; |
---|
| 347 | |
---|
| 348 | y = ymin - 0.5*dy ; |
---|
| 349 | yy = ymin - dy ; |
---|
| 350 | for (G4int i=0 ; i<NBIN; i++) |
---|
| 351 | { |
---|
| 352 | y += dy ; |
---|
| 353 | x = std::exp(y) ; |
---|
| 354 | yy += dy ; |
---|
| 355 | dx = std::exp(yy+dy)-std::exp(yy) ; |
---|
| 356 | |
---|
| 357 | ep = CutFixed*std::exp(c*x) ; |
---|
| 358 | |
---|
| 359 | CrossSection += ep*dx*ComputeDMicroscopicCrossSection(ParticleType, |
---|
| 360 | KineticEnergy,AtomicNumber, |
---|
| 361 | AtomicWeight,ep) ; |
---|
| 362 | if(nbin<NBIN) |
---|
| 363 | { |
---|
| 364 | nbin += 1 ; |
---|
| 365 | ya[nbin]=y ; |
---|
| 366 | proba[iz][it][nbin] = CrossSection ; |
---|
| 367 | } |
---|
| 368 | } |
---|
| 369 | ya[NBIN]=0. ; |
---|
| 370 | |
---|
| 371 | if(CrossSection > 0.) |
---|
| 372 | { |
---|
| 373 | for(G4int ib=0; ib<=nbin; ib++) |
---|
| 374 | { |
---|
| 375 | proba[iz][it][ib] /= CrossSection ; |
---|
| 376 | } |
---|
| 377 | } |
---|
| 378 | } |
---|
| 379 | } |
---|
| 380 | } |
---|
| 381 | |
---|
| 382 | G4VParticleChange* G4MuNuclearInteraction::PostStepDoIt( |
---|
| 383 | const G4Track& trackData, |
---|
| 384 | const G4Step& stepData) |
---|
| 385 | { |
---|
| 386 | static const G4double Mass=theMuonPlus->GetPDGMass() ; |
---|
| 387 | static const G4double m0=0.2*GeV ; |
---|
| 388 | |
---|
| 389 | aParticleChange.Initialize(trackData); |
---|
| 390 | G4Material* aMaterial=trackData.GetMaterial() ; |
---|
| 391 | |
---|
| 392 | const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle(); |
---|
| 393 | |
---|
| 394 | G4double KineticEnergy = aDynamicParticle->GetKineticEnergy(); |
---|
| 395 | G4ParticleMomentum ParticleDirection = |
---|
| 396 | aDynamicParticle->GetMomentumDirection(); |
---|
| 397 | |
---|
| 398 | // limits of the energy sampling |
---|
| 399 | G4double TotalEnergy = KineticEnergy + Mass ; |
---|
| 400 | G4double epmin = CutFixed ; |
---|
| 401 | G4double epmax = TotalEnergy - 0.5*proton_mass_c2 ; |
---|
| 402 | // check against insufficient energy |
---|
| 403 | if (epmax <= epmin ) |
---|
| 404 | { |
---|
| 405 | aParticleChange.ProposeMomentumDirection( ParticleDirection ); |
---|
| 406 | aParticleChange.ProposeEnergy( KineticEnergy ); |
---|
| 407 | aParticleChange.ProposeLocalEnergyDeposit (0.); |
---|
| 408 | aParticleChange.SetNumberOfSecondaries(0); |
---|
| 409 | return G4VDiscreteProcess::PostStepDoIt(trackData,stepData); |
---|
| 410 | } |
---|
| 411 | |
---|
| 412 | // select randomly one element constituing the material |
---|
| 413 | G4Element* anElement = SelectRandomAtom(aMaterial); |
---|
| 414 | |
---|
| 415 | // sample energy of the secondary ( pi0) |
---|
| 416 | // sampling using tables |
---|
| 417 | G4double ep,x,y ; |
---|
| 418 | G4int iy ; |
---|
| 419 | |
---|
| 420 | // select sampling table ; |
---|
| 421 | G4double lnZ = std::log(anElement->GetZ()) ; |
---|
| 422 | G4double delmin = 1.e10 ; |
---|
| 423 | G4double del ; |
---|
| 424 | G4int izz=0,itt=0,NBINminus1 ; |
---|
| 425 | NBINminus1 = NBIN-1 ; |
---|
| 426 | for (G4int iz=0; iz<nzdat; iz++) |
---|
| 427 | { |
---|
| 428 | del = std::abs(lnZ-std::log(zdat[iz])) ; |
---|
| 429 | if(del<delmin) |
---|
| 430 | { |
---|
| 431 | delmin=del ; |
---|
| 432 | izz=iz ; |
---|
| 433 | } |
---|
| 434 | } |
---|
| 435 | delmin = 1.e10 ; |
---|
| 436 | for (G4int it=0; it<ntdat; it++) |
---|
| 437 | { |
---|
| 438 | del = std::abs(std::log(KineticEnergy)-std::log(tdat[it])) ; |
---|
| 439 | if(del<delmin) |
---|
| 440 | { |
---|
| 441 | delmin=del; |
---|
| 442 | itt=it ; |
---|
| 443 | } |
---|
| 444 | } |
---|
| 445 | |
---|
| 446 | //sample energy transfer according to the sampling table |
---|
| 447 | |
---|
| 448 | G4double r = G4UniformRand() ; |
---|
| 449 | |
---|
| 450 | iy = -1 ; |
---|
| 451 | do { |
---|
| 452 | iy += 1 ; |
---|
| 453 | } while (((proba[izz][itt][iy]) < r)&&(iy < NBINminus1)) ; |
---|
| 454 | |
---|
| 455 | //sampling is Done uniformly in y in the bin |
---|
| 456 | if( iy < NBIN ) |
---|
| 457 | y = ya[iy] + G4UniformRand() * ( ya[iy+1] - ya[iy] ) ; |
---|
| 458 | else |
---|
| 459 | y = ya[iy] ; |
---|
| 460 | |
---|
| 461 | x = std::exp(y) ; |
---|
| 462 | ep = epmin*std::exp(x*std::log(epmax/epmin)) ; |
---|
| 463 | |
---|
| 464 | // sample scattering angle of mu, but first t should be sampled. |
---|
| 465 | G4double yy = ep/TotalEnergy ; |
---|
| 466 | G4double tmin=Mass*Mass*yy*yy/(1.-yy) ; |
---|
| 467 | G4double tmax=2.*proton_mass_c2*ep ; |
---|
| 468 | G4double t1,t2,t,w1,w2,w3,y1,y2,y3,rej ; |
---|
| 469 | if(m0<ep) |
---|
| 470 | { |
---|
| 471 | t1=m0*m0; |
---|
| 472 | t2=ep*ep; |
---|
| 473 | } |
---|
| 474 | else |
---|
| 475 | { |
---|
| 476 | t1=ep*ep; |
---|
| 477 | t2=m0*m0; |
---|
| 478 | } |
---|
| 479 | |
---|
| 480 | w1=tmax*t1; |
---|
| 481 | w2=tmax+t1 ; |
---|
| 482 | w3=tmax*(tmin+t1)/(tmin*w2); |
---|
| 483 | y1=1.-yy; |
---|
| 484 | y2=0.5*yy*yy; |
---|
| 485 | y3=y1+y2; |
---|
| 486 | |
---|
| 487 | // now the sampling of t |
---|
| 488 | G4int ntry=0; |
---|
| 489 | do |
---|
| 490 | { |
---|
| 491 | ntry += 1 ; |
---|
| 492 | t=w1/(w2*std::exp(G4UniformRand()*std::log(w3))-tmax) ; |
---|
| 493 | rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2)); |
---|
| 494 | } while (G4UniformRand() > rej) ; |
---|
| 495 | |
---|
| 496 | // compute angle from t |
---|
| 497 | G4double sinth2,theta ; // sinth2 = std::sin(theta/2)*std::sin(theta/2) ! |
---|
| 498 | sinth2 = 0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin) ; |
---|
| 499 | theta = std::acos(1.-2.*sinth2) ; |
---|
| 500 | |
---|
| 501 | G4double phi=twopi*G4UniformRand() ; |
---|
| 502 | G4double sinth=std::sin(theta) ; |
---|
| 503 | G4double dirx=sinth*std::cos(phi) , diry=sinth*std::sin(phi) , dirz=std::cos(theta); |
---|
| 504 | |
---|
| 505 | G4ThreeVector finalDirection(dirx,diry,dirz) ; |
---|
| 506 | finalDirection.rotateUz(ParticleDirection) ; |
---|
| 507 | |
---|
| 508 | G4double NewKinEnergy = KineticEnergy - ep ; |
---|
| 509 | G4double finalMomentum=std::sqrt(NewKinEnergy* |
---|
| 510 | (NewKinEnergy+2.*Mass)) ; |
---|
| 511 | |
---|
| 512 | G4double Ef=NewKinEnergy+Mass ; |
---|
| 513 | G4double initMomentum=std::sqrt(KineticEnergy*(TotalEnergy+Mass)) ; |
---|
| 514 | |
---|
| 515 | // G4double Q2=2.*(TotalEnergy*Ef-initMomentum*finalMomentum*std::cos(theta)-Mass*Mass) ; |
---|
| 516 | |
---|
| 517 | aParticleChange.ProposeMomentumDirection( finalDirection ); |
---|
| 518 | aParticleChange.ProposeEnergy( NewKinEnergy ); |
---|
| 519 | |
---|
| 520 | G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy); |
---|
| 521 | G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef); |
---|
| 522 | G4LorentzVector momentumTransfere = primaryMomentum-fsMomentum; |
---|
| 523 | |
---|
| 524 | G4DynamicParticle* aGamma = |
---|
| 525 | new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfere); |
---|
| 526 | G4Track gammaTrack(aGamma, trackData.GetGlobalTime(), trackData.GetPosition() ); |
---|
| 527 | gammaTrack.SetStep(trackData.GetStep()); |
---|
| 528 | G4Nucleus theTarget(aMaterial); |
---|
| 529 | |
---|
[1228] | 530 | G4VParticleChange* aHadronicFS = theHadronicVertex.ApplyYourself(theTarget, gammaTrack); |
---|
| 531 | //delete aGamma; |
---|
[819] | 532 | |
---|
| 533 | G4int numSecondaries = aHadronicFS->GetNumberOfSecondaries(); |
---|
| 534 | aParticleChange.SetNumberOfSecondaries(numSecondaries); |
---|
| 535 | |
---|
[1228] | 536 | // G4ParticleMomentum secondaryMomentum = G4ThreeVector(0.,0.,0.); |
---|
[819] | 537 | for(G4int iSec=0; iSec<numSecondaries; iSec++) |
---|
| 538 | { |
---|
[1228] | 539 | //secondaryMomentum |
---|
| 540 | // = secondaryMomentum + aHadronicFS->GetSecondary(iSec)->GetMomentum(); |
---|
[819] | 541 | aParticleChange.AddSecondary(aHadronicFS->GetSecondary(iSec)); |
---|
| 542 | } |
---|
[1340] | 543 | aHadronicFS->Clear(); |
---|
[819] | 544 | |
---|
| 545 | return G4VDiscreteProcess::PostStepDoIt(trackData,stepData); |
---|
| 546 | } |
---|
| 547 | |
---|
| 548 | G4Element* G4MuNuclearInteraction::SelectRandomAtom(G4Material* aMaterial) const |
---|
| 549 | { |
---|
| 550 | // select randomly 1 element within the material |
---|
| 551 | |
---|
[962] | 552 | G4int Index = aMaterial->GetIndex(); |
---|
| 553 | G4int NumberOfElements = aMaterial->GetNumberOfElements(); |
---|
[819] | 554 | const G4ElementVector* theElementVector = aMaterial->GetElementVector(); |
---|
[962] | 555 | if(1 == NumberOfElements) return ((*theElementVector)[0]); |
---|
[819] | 556 | |
---|
| 557 | G4double rval = G4UniformRand()*((*PartialSumSigma[Index])[NumberOfElements-1]); |
---|
[962] | 558 | for ( G4int i=0; i < NumberOfElements; i++ ) { |
---|
[819] | 559 | if (rval <= (*PartialSumSigma[Index])[i]) return ((*theElementVector)[i]); |
---|
[962] | 560 | } |
---|
| 561 | G4cout << "G4MuNuclearInteraction WARNING !!! no element selected for '" |
---|
| 562 | << aMaterial->GetName() |
---|
| 563 | << " 1st element returned." << G4endl; |
---|
| 564 | return ((*theElementVector)[0]); |
---|
[819] | 565 | } |
---|
[962] | 566 | |
---|
[819] | 567 | void G4MuNuclearInteraction::PrintInfoDefinition() |
---|
| 568 | { |
---|
| 569 | G4String comments = "cross sections from R. Kokoulin \n "; |
---|
| 570 | comments += " Good description up to 1000 PeV."; |
---|
| 571 | |
---|
| 572 | G4cout << G4endl << GetProcessName() << ": " << comments |
---|
| 573 | << "\n PhysicsTables from " << G4BestUnit(LowestKineticEnergy, |
---|
| 574 | "Energy") |
---|
| 575 | << " to " << G4BestUnit(HighestKineticEnergy,"Energy") |
---|
| 576 | << " in " << TotBin << " bins. \n"; |
---|
| 577 | |
---|
| 578 | G4cout << G4endl; |
---|
| 579 | } |
---|
| 580 | |
---|
[1228] | 581 | G4double G4MuNuclearInteraction::GetMeanFreePath(const G4Track& trackData, |
---|
| 582 | G4double, |
---|
| 583 | G4ForceCondition* condition) |
---|
| 584 | { |
---|
| 585 | const G4DynamicParticle* aDynamicParticle; |
---|
| 586 | G4Material* aMaterial; |
---|
| 587 | G4double MeanFreePath; |
---|
| 588 | G4bool isOutRange ; |
---|
| 589 | |
---|
| 590 | *condition = NotForced ; |
---|
| 591 | |
---|
| 592 | aDynamicParticle = trackData.GetDynamicParticle(); |
---|
| 593 | aMaterial = trackData.GetMaterial(); |
---|
| 594 | |
---|
| 595 | G4double KineticEnergy = aDynamicParticle->GetKineticEnergy(); |
---|
| 596 | |
---|
| 597 | if (KineticEnergy < LowestKineticEnergy) |
---|
| 598 | MeanFreePath = DBL_MAX ; |
---|
| 599 | else { |
---|
| 600 | if (KineticEnergy > HighestKineticEnergy) |
---|
| 601 | KineticEnergy = 0.99*HighestKineticEnergy ; |
---|
| 602 | MeanFreePath = (*theMeanFreePathTable)(aMaterial->GetIndex())-> |
---|
| 603 | GetValue( KineticEnergy, isOutRange ); |
---|
| 604 | } |
---|
| 605 | return MeanFreePath; |
---|
| 606 | } |
---|
| 607 | |
---|
| 608 | G4double G4MuNuclearInteraction::ComputeMeanFreePath( |
---|
| 609 | const G4ParticleDefinition* ParticleType, |
---|
| 610 | G4double KineticEnergy, |
---|
| 611 | const G4Material* aMaterial) |
---|
| 612 | { |
---|
| 613 | const G4ElementVector* theElementVector = aMaterial->GetElementVector() ; |
---|
| 614 | const G4double* theAtomNumDensityVector = |
---|
| 615 | aMaterial->GetAtomicNumDensityVector(); |
---|
| 616 | |
---|
| 617 | G4double SIGMA = 0 ; |
---|
| 618 | |
---|
| 619 | for ( size_t i=0 ; i < aMaterial->GetNumberOfElements() ; i++ ) |
---|
| 620 | { |
---|
| 621 | SIGMA += theAtomNumDensityVector[i] * |
---|
| 622 | ComputeMicroscopicCrossSection( ParticleType, KineticEnergy, |
---|
| 623 | (*theElementVector)[i]->GetZ(), |
---|
| 624 | (*theElementVector)[i]->GetA()) ; |
---|
| 625 | } |
---|
| 626 | |
---|
| 627 | return SIGMA<=0.0 ? DBL_MAX : 1./SIGMA ; |
---|
| 628 | } |
---|
| 629 | |
---|
| 630 | |
---|