Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of >> 2 // the GEANT4 collaboration. 1 // 3 // 2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work 3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement, 4 // * << 6 // and all its terms. 5 // * The Geant4 software is copyright of th << 7 // 6 // * the Geant4 Collaboration. It is provided << 8 // $Id: G4MuPairProduction.cc,v 1.15 2000/05/23 09:58:49 urban Exp $ 7 // * conditions of the Geant4 Software License << 9 // GEANT4 tag $Name: geant4-02-00 $ 8 // * LICENSE and available at http://cern.ch/ << 10 // 9 // * include a list of copyright holders. << 11 // -------------------------------------------------------------- 10 // * << 12 // GEANT 4 class implementation file 11 // * Neither the authors of this software syst << 13 // 12 // * institutes,nor the agencies providing fin << 14 // For information related to this code contact: 13 // * work make any representation or warran << 15 // CERN, CN Division, ASD group 14 // * regarding this software system or assum << 16 // History: first implementation, based on object model of 15 // * use. Please see the license in the file << 17 // 2nd December 1995, G.Cosmo 16 // * for the full disclaimer and the limitatio << 18 // -------- G4MuPairProduction physics process --------- 17 // * << 19 // by Laszlo Urban, May 1998 18 // * This code implementation is the result << 20 // ************************************************************** 19 // * technical work of the GEANT4 collaboratio << 21 // 04-06-98, in DoIt,secondary production condition:range>G4std::min(threshold,safety) 20 // * By using, copying, modifying or distri << 22 // 26/10/98, new stuff from R. Kokoulin + cleanup , L.Urban 21 // * any work based on the software) you ag << 23 // 06/05/99 , bug fixed , L.Urban 22 // * use in resulting scientific publicati << 24 // 10/02/00 modifications+bug fix , new e.m. structure, L.Urban 23 // * acceptance of all terms of the Geant4 Sof << 25 // -------------------------------------------------------------- 24 // ******************************************* << 25 // << 26 // << 27 // ------------------------------------------- << 28 // << 29 // GEANT4 Class file << 30 // << 31 // << 32 // File name: G4MuPairProduction << 33 // << 34 // Author: Laszlo Urban << 35 // << 36 // Creation date: 30.05.1998 << 37 // << 38 // Modifications: << 39 // << 40 // 04-06-98 in DoIt,secondary production condi << 41 // range>std::min(threshold,safety) << 42 // 26-10-98 new stuff from R. Kokoulin + clean << 43 // 06-05-99 bug fixed , L.Urban << 44 // 10-02-00 modifications+bug fix , new e.m. s << 45 // 29-05-01 V.Ivanchenko minor changes to prov << 46 // 10-08-01 new methods Store/Retrieve Physics << 47 // 17-09-01 migration of Materials to pure STL << 48 // 20-09-01 (L.Urban) in ComputeMicroscopicCro << 49 // if(MaxPairEnergy<CutInPairEnergy) << 50 // 26-09-01 completion of store/retrieve Physi << 51 // 28-09-01 suppression of theMuonPlus ..etc.. << 52 // 29-10-01 all static functions no more inlin << 53 // 07-11-01 particleMass becomes a local varia << 54 // 19-08-02 V.Ivanchenko update to new design << 55 // 23-12-02 Change interface in order to move << 56 // 26-12-02 Secondary production moved to deri << 57 // 13-02-03 SubCutoff regime is assigned to a << 58 // 08-08-03 STD substitute standard (V.Ivanch << 59 // 27-09-03 e+ set to be a secondary particle << 60 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossPro << 61 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko) << 62 // 17-08-04 Utilise mu+ tables for mu- (V.Ivan << 63 // 08-11-04 Migration to new interface of Stor << 64 // 08-04-05 Major optimisation of internal int << 65 // << 66 // ------------------------------------------- << 67 // << 68 //....oooOO0OOooo........oooOO0OOooo........oo << 69 //....oooOO0OOooo........oooOO0OOooo........oo << 70 26 71 #include "G4MuPairProduction.hh" 27 #include "G4MuPairProduction.hh" 72 #include "G4SystemOfUnits.hh" << 28 #include "G4EnergyLossTables.hh" 73 #include "G4Positron.hh" << 29 #include "G4UnitsTable.hh" 74 #include "G4VEmModel.hh" << 75 #include "G4MuPairProductionModel.hh" << 76 #include "G4ElementData.hh" << 77 #include "G4EmParameters.hh" << 78 30 79 //....oooOO0OOooo........oooOO0OOooo........oo << 31 // static members ........ 80 32 81 G4MuPairProduction::G4MuPairProduction(const G << 33 G4int G4MuPairProduction::nzdat = 5 ; 82 : G4VEnergyLossProcess(name), << 34 G4double G4MuPairProduction::zdat[]={1.,4.,13.,26.,92.}; 83 lowestKinEnergy(0.85*CLHEP::GeV) << 35 G4int G4MuPairProduction::ntdat = 8 ; >> 36 G4double G4MuPairProduction::tdat[]={1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10}; >> 37 G4int G4MuPairProduction::NBIN = 1000 ; //100 ; >> 38 G4double G4MuPairProduction::ya[1001]={0.}; >> 39 G4double G4MuPairProduction::proba[5][8][1001]={0.}; >> 40 G4double G4MuPairProduction::MinPairEnergy = 4.*electron_mass_c2 ; >> 41 >> 42 G4double G4MuPairProduction::LowerBoundLambda = 1.*keV ; >> 43 G4double G4MuPairProduction::UpperBoundLambda = 1000000.*TeV ; >> 44 G4int G4MuPairProduction::NbinLambda = 150 ; >> 45 >> 46 >> 47 G4MuPairProduction::G4MuPairProduction(const G4String& processName) >> 48 : G4VMuEnergyLoss(processName), >> 49 theMeanFreePathTable(NULL) >> 50 { } >> 51 >> 52 >> 53 G4MuPairProduction::~G4MuPairProduction() 84 { 54 { 85 SetProcessSubType(fPairProdByCharged); << 55 if (theMeanFreePathTable) { 86 SetSecondaryParticle(G4Positron::Positron()) << 56 theMeanFreePathTable->clearAndDestroy(); 87 SetIonisation(false); << 57 delete theMeanFreePathTable; >> 58 } >> 59 >> 60 if (&PartialSumSigma) { >> 61 PartialSumSigma.clearAndDestroy(); >> 62 } 88 } 63 } 89 64 90 //....oooOO0OOooo........oooOO0OOooo........oo << 65 void G4MuPairProduction::BuildPhysicsTable( 91 << 66 const G4ParticleDefinition& aParticleType) 92 G4bool G4MuPairProduction::IsApplicable(const << 67 // just call BuildLossTable+BuildLambdaTable 93 { 68 { 94 return (p.GetPDGCharge() != 0.0); << 69 >> 70 // get bining from EnergyLoss >> 71 LowestKineticEnergy = GetLowerBoundEloss() ; >> 72 HighestKineticEnergy = GetUpperBoundEloss() ; >> 73 TotBin = GetNbinEloss() ; >> 74 >> 75 BuildLossTable(aParticleType) ; >> 76 >> 77 if(&aParticleType==theMuonMinus) >> 78 { >> 79 RecorderOfmuminusProcess[CounterOfmuminusProcess] = (*this).theLossTable ; >> 80 CounterOfmuminusProcess++; >> 81 } >> 82 else >> 83 { >> 84 RecorderOfmuplusProcess[CounterOfmuplusProcess] = (*this).theLossTable ; >> 85 CounterOfmuplusProcess++; >> 86 } >> 87 >> 88 // sampling table should be made only once ! >> 89 if(theMeanFreePathTable == NULL) >> 90 MakeSamplingTables(&aParticleType) ; >> 91 >> 92 G4double electronCutInRange = G4Electron::Electron()->GetCuts(); >> 93 if(electronCutInRange != lastelectronCutInRange) >> 94 BuildLambdaTable(aParticleType) ; >> 95 >> 96 G4VMuEnergyLoss::BuildDEDXTable(aParticleType) ; >> 97 >> 98 if(&aParticleType==theMuonPlus) >> 99 PrintInfoDefinition() ; 95 } 100 } 96 101 97 //....oooOO0OOooo........oooOO0OOooo........oo << 102 void G4MuPairProduction::BuildLossTable( >> 103 const G4ParticleDefinition& aParticleType) >> 104 { >> 105 G4double KineticEnergy,TotalEnergy,pairloss,Z, >> 106 loss,natom,eCut,pCut ; >> 107 >> 108 const G4MaterialTable* theMaterialTable = >> 109 G4Material::GetMaterialTable(); >> 110 ParticleMass = aParticleType.GetPDGMass() ; >> 111 ElectronCutInKineticEnergy = (*theElectron).GetEnergyCuts() ; >> 112 PositronCutInKineticEnergy = (*thePositron).GetEnergyCuts() ; >> 113 >> 114 G4int numOfMaterials = theMaterialTable->length() ; >> 115 >> 116 if (theLossTable) { >> 117 theLossTable->clearAndDestroy(); >> 118 delete theLossTable; >> 119 } >> 120 >> 121 theLossTable = new G4PhysicsTable(numOfMaterials) ; >> 122 >> 123 for (G4int J=0; J<numOfMaterials; J++) >> 124 { >> 125 G4PhysicsLogVector* aVector = new G4PhysicsLogVector( >> 126 LowestKineticEnergy,HighestKineticEnergy,TotBin); >> 127 ElectronCutInKineticEnergyNow = ElectronCutInKineticEnergy[J] ; >> 128 PositronCutInKineticEnergyNow = PositronCutInKineticEnergy[J] ; >> 129 const G4Material* material = (*theMaterialTable)[J] ; >> 130 const G4ElementVector* theElementVector = >> 131 material->GetElementVector() ; >> 132 const G4double* theAtomicNumDensityVector = >> 133 material->GetAtomicNumDensityVector() ; >> 134 const G4int NumberOfElements = >> 135 material->GetNumberOfElements() ; >> 136 >> 137 for (G4int i=0; i<TotBin; i++) >> 138 { >> 139 KineticEnergy = aVector->GetLowEdgeEnergy(i) ; >> 140 TotalEnergy = KineticEnergy+ParticleMass ; >> 141 eCut = ElectronCutInKineticEnergyNow ; >> 142 pCut = PositronCutInKineticEnergyNow ; >> 143 >> 144 if(eCut>KineticEnergy) >> 145 eCut = KineticEnergy ; >> 146 if(pCut>KineticEnergy) >> 147 pCut = KineticEnergy ; >> 148 >> 149 pairloss = 0.; >> 150 for (G4int iel=0; iel<NumberOfElements; iel++) >> 151 { >> 152 Z=(*theElementVector)(iel)->GetZ(); >> 153 natom = theAtomicNumDensityVector[iel] ; >> 154 loss = ComputePairLoss(&aParticleType, >> 155 Z,KineticEnergy,eCut,pCut) ; >> 156 pairloss += natom*loss ; >> 157 } >> 158 if(pairloss<0.) >> 159 pairloss = 0. ; >> 160 aVector->PutValue(i,pairloss); >> 161 } 98 162 99 G4double G4MuPairProduction::MinPrimaryEnergy( << 163 theLossTable->insert(aVector); 100 const G4Material*, << 164 } 101 G4double) << 165 } >> 166 >> 167 G4double G4MuPairProduction::ComputePairLoss( >> 168 const G4ParticleDefinition* ParticleType, >> 169 G4double AtomicNumber, >> 170 G4double KineticEnergy, >> 171 G4double ElectronEnergyCut, >> 172 G4double PositronEnergyCut) 102 { 173 { 103 return lowestKinEnergy; << 174 static const G4double >> 175 xgi[] ={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 }; >> 176 static const G4double >> 177 wgi[] ={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 }; >> 178 static const G4double ak1=6.9 ; >> 179 static const G4double ak2=1.0 ; >> 180 static const G4double sqrte = sqrt(exp(1.)) ; >> 181 G4double z13 = exp(log(AtomicNumber)/3.) ; >> 182 >> 183 G4double loss = 0.0 ; >> 184 >> 185 if ( AtomicNumber < 1. ) return loss; >> 186 >> 187 G4double CutInPairEnergy = ElectronEnergyCut+PositronEnergyCut >> 188 +2.*electron_mass_c2 ; >> 189 if( CutInPairEnergy <= MinPairEnergy ) return loss ; >> 190 >> 191 G4double MaxPairEnergy = KineticEnergy+ParticleMass*(1.-0.75*sqrte*z13) ; >> 192 if(MaxPairEnergy < MinPairEnergy) >> 193 MaxPairEnergy = MinPairEnergy ; >> 194 >> 195 if( CutInPairEnergy >= MaxPairEnergy ) >> 196 CutInPairEnergy = MaxPairEnergy ; >> 197 >> 198 if(CutInPairEnergy <= MinPairEnergy) return loss ; >> 199 >> 200 G4double aaa,bbb,hhh,x,epln,ep ; >> 201 G4int kkk ; >> 202 // calculate the rectricted loss >> 203 // numerical integration in log(PairEnergy) >> 204 aaa = log(MinPairEnergy) ; >> 205 bbb = log(CutInPairEnergy) ; >> 206 kkk = int((bbb-aaa)/ak1+ak2) ; >> 207 hhh = (bbb-aaa)/kkk ; >> 208 >> 209 for (G4int l=0 ; l<kkk; l++) >> 210 { >> 211 x = aaa+hhh*l ; >> 212 for (G4int ll=0; ll<8; ll++) >> 213 { >> 214 epln=x+xgi[ll]*hhh ; >> 215 ep = exp(epln) ; >> 216 loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(ParticleType, >> 217 KineticEnergy,AtomicNumber, >> 218 ep) ; >> 219 } >> 220 } >> 221 loss *= hhh ; >> 222 if (loss < 0.) loss = 0.; >> 223 return loss ; 104 } 224 } >> 225 105 226 106 //....oooOO0OOooo........oooOO0OOooo........oo << 227 void G4MuPairProduction::BuildLambdaTable( >> 228 const G4ParticleDefinition& ParticleType) >> 229 { >> 230 G4double LowEdgeEnergy , Value; >> 231 G4double FixedEnergy = (LowestKineticEnergy + HighestKineticEnergy)/2. ; >> 232 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ; >> 233 >> 234 if (theMeanFreePathTable) { >> 235 theMeanFreePathTable->clearAndDestroy(); >> 236 delete theMeanFreePathTable; >> 237 } >> 238 theMeanFreePathTable = new >> 239 G4PhysicsTable( G4Material::GetNumberOfMaterials() ) ; >> 240 PartialSumSigma.resize(G4Material::GetNumberOfMaterials()); >> 241 >> 242 G4PhysicsLogVector* ptrVector; >> 243 for ( G4int J=0 ; J < G4Material::GetNumberOfMaterials(); J++ ) >> 244 { >> 245 ptrVector = new G4PhysicsLogVector( >> 246 LowerBoundLambda,UpperBoundLambda,NbinLambda); >> 247 >> 248 >> 249 const G4Material* material= (*theMaterialTable)[J]; >> 250 >> 251 for ( G4int i = 0 ; i < NbinLambda ; i++ ) >> 252 { >> 253 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ; >> 254 Value = ComputeMeanFreePath( &ParticleType, LowEdgeEnergy, >> 255 material ); >> 256 ptrVector->PutValue( i , Value ) ; >> 257 } >> 258 >> 259 theMeanFreePathTable->insertAt( J , ptrVector ); >> 260 >> 261 // Compute the PartialSumSigma table at a given fixed energy >> 262 ComputePartialSumSigma( &ParticleType, FixedEnergy, material) ; >> 263 } >> 264 } 107 265 108 void G4MuPairProduction::InitialiseEnergyLossP << 266 void G4MuPairProduction::ComputePartialSumSigma( 109 const G4ParticleDefin << 267 const G4ParticleDefinition* ParticleType, 110 const G4ParticleDefinition*) << 268 G4double KineticEnergy, >> 269 const G4Material* aMaterial) 111 { 270 { 112 if (isInitialised) { return; } << 271 G4int Imate = aMaterial->GetIndex(); 113 isInitialised = true; << 272 G4int NbOfElements = aMaterial->GetNumberOfElements(); >> 273 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); >> 274 const G4double* theAtomNumDensityVector = aMaterial-> >> 275 GetAtomicNumDensityVector(); >> 276 G4double ElectronEnergyCut = (G4Electron::GetCutsInEnergy())[Imate]; >> 277 G4double PositronEnergyCut = (G4Positron::GetCutsInEnergy())[Imate]; >> 278 >> 279 PartialSumSigma(Imate) = new G4ValVector(NbOfElements); >> 280 >> 281 G4double SIGMA = 0. ; >> 282 >> 283 for ( G4int Ielem=0 ; Ielem < NbOfElements ; Ielem++ ) >> 284 { >> 285 SIGMA += theAtomNumDensityVector[Ielem] * >> 286 ComputeMicroscopicCrossSection( ParticleType, KineticEnergy, >> 287 (*theElementVector)(Ielem)->GetZ(), >> 288 ElectronEnergyCut,PositronEnergyCut ); 114 289 115 theParticle = part; << 290 PartialSumSigma(Imate)->insertAt(Ielem, SIGMA); >> 291 } >> 292 } 116 293 117 G4VEmModel* mod = EmModel(0); << 294 G4double G4MuPairProduction::ComputeMicroscopicCrossSection( 118 if(nullptr == mod) { << 295 const G4ParticleDefinition* ParticleType, 119 lowestKinEnergy = std::max(lowestKinEnergy << 296 G4double KineticEnergy, 120 auto ptr = new G4MuPairProductionModel(par << 297 G4double AtomicNumber, 121 ptr->SetLowestKineticEnergy(lowestKinEnerg << 298 G4double ElectronEnergyCut, 122 mod = ptr; << 299 G4double PositronEnergyCut) 123 SetEmModel(mod); << 300 >> 301 { >> 302 static const G4double >> 303 xgi[] ={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 }; >> 304 static const G4double >> 305 wgi[] ={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 }; >> 306 static const G4double ak1=6.9 ; >> 307 static const G4double ak2=1.0 ; >> 308 >> 309 static const G4double sqrte = sqrt(exp(1.)) ; >> 310 G4double z13 = exp(log(AtomicNumber)/3.) ; >> 311 >> 312 G4double CrossSection = 0.0 ; >> 313 >> 314 if ( AtomicNumber < 1. ) return CrossSection; >> 315 >> 316 G4double CutInPairEnergy = ElectronEnergyCut+PositronEnergyCut >> 317 +2.*electron_mass_c2 ; >> 318 >> 319 if( CutInPairEnergy < 4.*electron_mass_c2 ) >> 320 CutInPairEnergy = 4.*electron_mass_c2 ; >> 321 >> 322 G4double MaxPairEnergy = KineticEnergy+ParticleMass*(1.-0.75*sqrte*z13) ; >> 323 if(MaxPairEnergy < CutInPairEnergy) >> 324 MaxPairEnergy = CutInPairEnergy ; >> 325 if( CutInPairEnergy >= MaxPairEnergy ) return CrossSection ; >> 326 >> 327 G4double aaa,bbb,hhh,x,epln,ep ; >> 328 G4int kkk ; >> 329 // calculate the total cross section >> 330 // numerical integration in log(PairEnergy) >> 331 aaa = log(CutInPairEnergy) ; >> 332 bbb = log(MaxPairEnergy) ; >> 333 kkk = int((bbb-aaa)/ak1+ak2) ; >> 334 hhh = (bbb-aaa)/kkk ; >> 335 for (G4int l=0 ; l<kkk; l++) >> 336 { >> 337 x = aaa+hhh*l ; >> 338 for (G4int ll=0; ll<8; ll++) >> 339 { >> 340 epln = x+xgi[ll]*hhh; >> 341 ep = exp(epln) ; >> 342 CrossSection += wgi[ll]*ep*ComputeDMicroscopicCrossSection(ParticleType, >> 343 KineticEnergy,AtomicNumber, >> 344 ep) ; >> 345 } 124 } 346 } >> 347 CrossSection *= hhh ; 125 348 126 G4VEmFluctuationModel* fm = nullptr; << 349 if (CrossSection < 0.) CrossSection = 0.; 127 G4EmParameters* param = G4EmParameters::Inst << 128 mod->SetLowEnergyLimit(param->MinKinEnergy() << 129 mod->SetHighEnergyLimit(param->MaxKinEnergy( << 130 mod->SetSecondaryThreshold(param->MuHadBrems << 131 AddEmModel(1, mod, fm); << 132 } << 133 350 134 //....oooOO0OOooo........oooOO0OOooo........oo << 351 return CrossSection; >> 352 } 135 353 136 void G4MuPairProduction::StreamProcessInfo(std << 354 void G4MuPairProduction::MakeSamplingTables( >> 355 const G4ParticleDefinition* ParticleType) 137 { 356 { 138 G4ElementData* ed = EmModel()->GetElementDat << 357 G4int nbin; 139 if(ed) { << 358 G4double AtomicNumber,KineticEnergy ; 140 for(G4int Z=1; Z<93; ++Z) { << 359 G4double c,y,ymin,ymax,dy,yy,dx,x,ep ; 141 G4Physics2DVector* pv = ed->GetElement2D << 360 142 if(pv) { << 361 static const G4double sqrte = sqrt(exp(1.)) ; 143 out << " Sampling table " << pv-> << 362 144 << "x" << pv->GetLengthX() << "; from " << 363 for (G4int iz=0; iz<nzdat; iz++) 145 << std::exp(pv->GetY(0))/GeV << " GeV to << 364 { 146 << std::exp(pv->GetY(pv->GetLengthY()-1) << 365 AtomicNumber = zdat[iz]; 147 << " TeV " << G4endl; << 366 G4double z13 = exp(log(AtomicNumber)/3.) ; 148 break; << 367 >> 368 for (G4int it=0; it<ntdat; it++) >> 369 { >> 370 KineticEnergy = tdat[it]; >> 371 G4double MaxPairEnergy = KineticEnergy+ParticleMass*(1.-0.75*sqrte*z13) ; >> 372 >> 373 G4double CrossSection = 0.0 ; >> 374 >> 375 G4int NbofIntervals ; >> 376 c = log(MaxPairEnergy/MinPairEnergy) ; >> 377 >> 378 ymin = -5. ; >> 379 ymax = 0. ; >> 380 dy = (ymax-ymin)/NBIN ; >> 381 nbin=-1; >> 382 y = ymin - 0.5*dy ; >> 383 yy = ymin - dy ; >> 384 for (G4int i=0 ; i<NBIN; i++) >> 385 { >> 386 y += dy ; >> 387 x = exp(y) ; >> 388 yy += dy ; >> 389 dx = exp(yy+dy)-exp(yy) ; >> 390 ep = MinPairEnergy*exp(c*x) ; >> 391 CrossSection += ep*dx*ComputeDMicroscopicCrossSection(ParticleType, >> 392 KineticEnergy,AtomicNumber,ep); >> 393 if(nbin<NBIN) >> 394 { >> 395 nbin += 1 ; >> 396 ya[nbin]=y ; >> 397 proba[iz][it][nbin] = CrossSection ; >> 398 } >> 399 } >> 400 ya[NBIN]=0. ; >> 401 >> 402 if(CrossSection > 0.) >> 403 { >> 404 for(G4int ib=0; ib<=nbin; ib++) >> 405 { >> 406 proba[iz][it][ib] /= CrossSection ; >> 407 >> 408 } 149 } 409 } 150 } 410 } 151 } 411 } >> 412 } >> 413 >> 414 >> 415 G4double G4MuPairProduction::ComputeDDMicroscopicCrossSection( >> 416 const G4ParticleDefinition* ParticleType, >> 417 G4double KineticEnergy, G4double AtomicNumber, >> 418 G4double PairEnergy,G4double asymmetry) >> 419 // Calculates the double differential (DD) microscopic cross section >> 420 // using the cross section formula of R.P. Kokoulin (18/01/98) >> 421 { >> 422 static const G4double sqrte = sqrt(exp(1.)) ; >> 423 >> 424 G4double bbbtf= 183. ; >> 425 G4double bbbh = 202.4 ; >> 426 G4double g1tf = 1.95e-5 ; >> 427 G4double g2tf = 5.3e-5 ; >> 428 G4double g1h = 4.4e-5 ; >> 429 G4double g2h = 4.8e-5 ; >> 430 >> 431 G4double massratio = ParticleMass/electron_mass_c2 ; >> 432 G4double massratio2 = massratio*massratio ; >> 433 G4double TotalEnergy = KineticEnergy + ParticleMass ; >> 434 G4double z13 = exp(log(AtomicNumber)/3.) ; >> 435 G4double z23 = z13*z13 ; >> 436 G4double EnergyLoss = TotalEnergy - PairEnergy ; >> 437 >> 438 G4double c3 = 3.*sqrte*ParticleMass/4. ; >> 439 >> 440 G4double DDCrossSection = 0. ; >> 441 >> 442 if(EnergyLoss <= c3*z13) >> 443 return DDCrossSection ; >> 444 >> 445 G4double c7 = 4.*electron_mass_c2 ; >> 446 G4double c8 = 6.*ParticleMass*ParticleMass ; >> 447 G4double alf = c7/PairEnergy ; >> 448 G4double a3 = 1. - alf ; >> 449 >> 450 if(a3 <= 0.) >> 451 return DDCrossSection ; >> 452 >> 453 // zeta calculation >> 454 G4double bbb,g1,g2,zeta1,zeta2,zeta,z2 ; >> 455 if( AtomicNumber < 1.5 ) >> 456 { >> 457 bbb = bbbh ; >> 458 g1 = g1h ; >> 459 g2 = g2h ; >> 460 } >> 461 else >> 462 { >> 463 bbb = bbbtf ; >> 464 g1 = g1tf ; >> 465 g2 = g2tf ; >> 466 } >> 467 zeta1 = 0.073 * log(TotalEnergy/(ParticleMass+g1*z23*TotalEnergy))-0.26 ; >> 468 if( zeta1 > 0.) >> 469 { >> 470 zeta2 = 0.058*log(TotalEnergy/(ParticleMass+g2*z13*TotalEnergy))-0.14 ; >> 471 zeta = zeta1/zeta2 ; >> 472 } >> 473 else >> 474 { >> 475 zeta = 0. ; >> 476 } >> 477 >> 478 z2 = AtomicNumber*(AtomicNumber+zeta) ; >> 479 >> 480 G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*PairEnergy) ; >> 481 G4double a0 = TotalEnergy*EnergyLoss ; >> 482 G4double a1 = PairEnergy*PairEnergy/a0 ; >> 483 G4double bet = 0.5*a1 ; >> 484 G4double xi0 = 0.25*massratio2*a1 ; >> 485 G4double del = c8/a0 ; >> 486 >> 487 G4double romin = 0. ; >> 488 G4double romax = (1.-del)*sqrt(1.-c7/PairEnergy) ; >> 489 >> 490 if((asymmetry < romin) || (asymmetry > romax)) >> 491 return DDCrossSection ; >> 492 >> 493 G4double a4 = 1.-asymmetry ; >> 494 G4double a5 = a4*(2.-a4) ; >> 495 G4double a6 = 1.-a5 ; >> 496 G4double a7 = 1.+a6 ; >> 497 G4double a9 = 3.+a6 ; >> 498 G4double xi = xi0*a5 ; >> 499 G4double xii = 1./xi ; >> 500 G4double xi1 = 1.+xi ; >> 501 G4double screen = screen0*xi1/a5 ; >> 502 >> 503 G4double yeu = 5.-a6+4.*bet*a7 ; >> 504 G4double yed = 2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6) ; >> 505 G4double yel = 1.+yeu/yed ; >> 506 G4double ale=log(bbb/z13*sqrt(xi1*yel)/(1.+screen*yel)) ; >> 507 G4double cre = 0.5*log(1.+2.25/(massratio2*z23)*xi1*yel) ; >> 508 G4double be ; >> 509 if(xi <= 1.e3) >> 510 be = ((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9; >> 511 else >> 512 be = (3.-a6+a1*a7)/(2.+xi) ; >> 513 G4double fe = (ale-cre)*be ; >> 514 if( fe < 0.) >> 515 fe = 0. ; >> 516 >> 517 G4double ymu = 4.+a6 +3.*bet*a7 ; >> 518 G4double ymd = a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6 ; >> 519 G4double ym1 = 1.+ymu/ymd ; >> 520 G4double alm_crm = log(bbb*massratio/(1.5*z23*(1.+screen*ym1))) ; >> 521 G4double a10,bm ; >> 522 if( xi >= 1.e-3) >> 523 { >> 524 a10 = (1.+a1)*a5 ; >> 525 bm = (a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10 ; >> 526 } >> 527 else >> 528 bm = (5.-a6+bet*a9)*(xi/2.) ; >> 529 G4double fm = alm_crm*bm ; >> 530 if( fm < 0.) >> 531 fm = 0. ; >> 532 >> 533 DDCrossSection = (fe+fm/massratio2) ; >> 534 >> 535 DDCrossSection *= 4.*fine_structure_const*fine_structure_const >> 536 *classic_electr_radius*classic_electr_radius/(3.*pi) ; >> 537 >> 538 DDCrossSection *= z2*EnergyLoss/(TotalEnergy*PairEnergy) ; >> 539 >> 540 >> 541 return DDCrossSection ; >> 542 >> 543 } >> 544 >> 545 G4double G4MuPairProduction::GetDMicroscopicCrossSection( >> 546 const G4ParticleDefinition* ParticleType, >> 547 G4double KineticEnergy, G4double AtomicNumber, >> 548 G4double PairEnergy) >> 549 { >> 550 return ComputeDMicroscopicCrossSection(ParticleType,KineticEnergy, >> 551 AtomicNumber,PairEnergy) ; >> 552 } >> 553 >> 554 G4double G4MuPairProduction::ComputeDMicroscopicCrossSection( >> 555 const G4ParticleDefinition* ParticleType, >> 556 G4double KineticEnergy, G4double AtomicNumber, >> 557 G4double PairEnergy) >> 558 // Calculates the differential (D) microscopic cross section >> 559 // using the cross section formula of R.P. Kokoulin (18/01/98) >> 560 { >> 561 >> 562 static const G4double >> 563 xgi[] ={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 }; >> 564 >> 565 static const G4double >> 566 wgi[] ={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 }; >> 567 >> 568 G4double DCrossSection = 0. ; >> 569 >> 570 G4double TotalEnergy = KineticEnergy + ParticleMass ; >> 571 G4double EnergyLoss = TotalEnergy - PairEnergy ; >> 572 G4double a = 6.*ParticleMass*ParticleMass/(TotalEnergy*EnergyLoss) ; >> 573 G4double b = 4.*electron_mass_c2/PairEnergy ; >> 574 >> 575 if((b+2.*a*(1.-b))/(1.+(1.-a)*sqrt(1.-b)) <= 0.) return DCrossSection ; >> 576 >> 577 G4double tmn=log((b+2.*a*(1.-b))/(1.+(1.-a)*sqrt(1.-b))) ; >> 578 >> 579 // G4double DCrossSection = 0. ; >> 580 G4double ro ; >> 581 // Gaussian integration in ln(1-ro) ( with 8 points) >> 582 for (G4int i=0; i<7; i++) >> 583 { >> 584 ro = 1.-exp(tmn*xgi[i]) ; >> 585 >> 586 DCrossSection += (1.-ro)*ComputeDDMicroscopicCrossSection( >> 587 ParticleType,KineticEnergy, >> 588 AtomicNumber,PairEnergy,ro) >> 589 *wgi[i] ; >> 590 } >> 591 >> 592 DCrossSection *= -tmn ; >> 593 >> 594 return DCrossSection ; >> 595 152 } 596 } >> 597 >> 598 G4VParticleChange* G4MuPairProduction::PostStepDoIt(const G4Track& trackData, >> 599 const G4Step& stepData) >> 600 { >> 601 static const G4double esq = sqrt(exp(1.)); 153 602 154 //....oooOO0OOooo........oooOO0OOooo........oo << 603 aParticleChange.Initialize(trackData); >> 604 G4Material* aMaterial=trackData.GetMaterial() ; >> 605 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle(); >> 606 G4double KineticEnergy = aDynamicParticle->GetKineticEnergy(); >> 607 G4ParticleMomentum ParticleDirection = >> 608 aDynamicParticle->GetMomentumDirection(); >> 609 >> 610 // e-e+ cut in this material >> 611 G4double ElectronEnergyCut = electron_mass_c2+ >> 612 ((*G4Electron::Electron()).GetCutsInEnergy())[aMaterial->GetIndex()]; >> 613 G4double PositronEnergyCut = electron_mass_c2+ >> 614 ((*G4Positron::Positron()).GetCutsInEnergy())[aMaterial->GetIndex()]; >> 615 G4double CutInPairEnergy = ElectronEnergyCut + PositronEnergyCut ; >> 616 >> 617 if (CutInPairEnergy < MinPairEnergy) CutInPairEnergy = MinPairEnergy ; >> 618 >> 619 // check against insufficient energy >> 620 if(KineticEnergy < CutInPairEnergy ) >> 621 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData); >> 622 >> 623 // select randomly one element constituing the material >> 624 G4Element* anElement = SelectRandomAtom(aMaterial); >> 625 >> 626 // limits of the energy sampling >> 627 G4double TotalEnergy = KineticEnergy + ParticleMass ; >> 628 G4double TotalMomentum = sqrt(KineticEnergy*(TotalEnergy+ParticleMass)) ; >> 629 G4double Z3 = anElement->GetIonisation()->GetZ3() ; >> 630 G4double MaxPairEnergy = TotalEnergy-0.75*esq*ParticleMass*Z3 ; >> 631 >> 632 if(MinPairEnergy >= MaxPairEnergy) >> 633 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData); >> 634 >> 635 // sample e-e+ energy, pair energy first >> 636 G4double PairEnergy,xc,x,yc,y ; >> 637 G4int iZ,iT,iy ; >> 638 >> 639 // select sampling table ; >> 640 G4double lnZ = log(anElement->GetZ()) ; >> 641 G4double delmin = 1.e10 ; >> 642 G4double del ; >> 643 G4int izz,itt,NBINminus1 ; >> 644 NBINminus1 = NBIN-1 ; >> 645 for (G4int iz=0; iz<nzdat; iz++) >> 646 { >> 647 del = abs(lnZ-log(zdat[iz])) ; >> 648 if(del<delmin) >> 649 { >> 650 delmin=del ; >> 651 izz=iz ; >> 652 } >> 653 } >> 654 delmin = 1.e10 ; >> 655 for (G4int it=0; it<ntdat; it++) >> 656 { >> 657 del = abs(log(KineticEnergy)-log(tdat[it])) ; >> 658 if(del<delmin) >> 659 { >> 660 delmin=del; >> 661 itt=it ; >> 662 } >> 663 } >> 664 >> 665 if( CutInPairEnergy <= MinPairEnergy) >> 666 iy = 0 ; >> 667 else >> 668 { >> 669 xc = log(CutInPairEnergy/MinPairEnergy)/log(MaxPairEnergy/MinPairEnergy) ; >> 670 yc = log(xc) ; >> 671 >> 672 iy = -1 ; >> 673 do { >> 674 iy += 1 ; >> 675 } while ((ya[iy] < yc )&&(iy < NBINminus1)) ; >> 676 } >> 677 >> 678 G4double norm = proba[izz][itt][iy] ; >> 679 >> 680 G4double r = norm+G4UniformRand()*(1.-norm) ; >> 681 >> 682 iy -= 1 ; >> 683 do { >> 684 iy += 1 ; >> 685 } while ((proba[izz][itt][iy] < r)&&(iy < NBINminus1)) ; >> 686 >> 687 //sampling is uniformly in y in the bin >> 688 if( iy < NBIN ) >> 689 y = ya[iy] + G4UniformRand() * ( ya[iy+1] - ya[iy]) ; >> 690 else >> 691 y = ya[iy] ; >> 692 >> 693 x = exp(y) ; >> 694 >> 695 PairEnergy = MinPairEnergy*exp(x*log(MaxPairEnergy/MinPairEnergy)) ; >> 696 >> 697 // sample r=(E+-E-)/PairEnergy ( uniformly .....) >> 698 G4double rmax = (1.-6.*ParticleMass*ParticleMass/(TotalEnergy* >> 699 (TotalEnergy-PairEnergy))) >> 700 *sqrt(1.-MinPairEnergy/PairEnergy) ; >> 701 r = rmax * (-1.+2.*G4UniformRand()) ; >> 702 >> 703 // compute energies from PairEnergy,r >> 704 G4double ElectronEnergy=(1.-r)*PairEnergy/2. ; >> 705 G4double PositronEnergy=(1.+r)*PairEnergy/2. ; >> 706 >> 707 // angles of the emitted particles ( Z - axis along the parent particle) >> 708 // (mean theta for the moment) >> 709 G4double Teta = electron_mass_c2/TotalEnergy ; >> 710 >> 711 G4double Phi = twopi * G4UniformRand() ; >> 712 G4double dirx = sin(Teta)*cos(Phi) , diry = sin(Teta)*sin(Phi) , >> 713 dirz = cos(Teta) ; >> 714 >> 715 G4double LocalEnerDeposit = 0. ; >> 716 G4int numberofsecondaries = 1 ; >> 717 G4int flagelectron = 0 ; >> 718 G4int flagpositron = 1 ; >> 719 G4DynamicParticle *aParticle1,*aParticle2 ; >> 720 G4double ElectronMomentum , PositronMomentum ; >> 721 G4double finalPx,finalPy,finalPz ; >> 722 >> 723 G4double ElectKineEnergy = ElectronEnergy - electron_mass_c2 ; >> 724 >> 725 if((ElectKineEnergy > ElectronEnergyCut) || >> 726 (G4EnergyLossTables::GetRange( >> 727 G4Electron::Electron(),ElectKineEnergy,aMaterial) >= >> 728 stepData.GetPostStepPoint()->GetSafety())) >> 729 { >> 730 numberofsecondaries += 1 ; >> 731 flagelectron = 1 ; >> 732 ElectronMomentum = sqrt(ElectKineEnergy* >> 733 (ElectronEnergy+electron_mass_c2)); >> 734 G4ThreeVector ElectDirection ( dirx, diry, dirz ); >> 735 ElectDirection.rotateUz(ParticleDirection); >> 736 >> 737 // create G4DynamicParticle object for the particle1 >> 738 aParticle1= new G4DynamicParticle (G4Electron::Electron(), >> 739 ElectDirection, ElectKineEnergy); >> 740 } >> 741 else >> 742 { LocalEnerDeposit += ElectKineEnergy ; } >> 743 >> 744 // the e+ is always created (even with Ekine=0) for further annihilation. >> 745 >> 746 G4double PositKineEnergy = PositronEnergy - electron_mass_c2 ; >> 747 PositronMomentum = sqrt(PositKineEnergy*(PositronEnergy+electron_mass_c2)); >> 748 >> 749 if((PositKineEnergy < PositronEnergyCut) && >> 750 (G4EnergyLossTables::GetRange( >> 751 G4Positron::Positron(),PositKineEnergy,aMaterial) <= >> 752 stepData.GetPostStepPoint()->GetSafety())) >> 753 { >> 754 LocalEnerDeposit += PositKineEnergy ; >> 755 PositKineEnergy = 0. ; >> 756 } >> 757 G4ThreeVector PositDirection ( -dirx, -diry, dirz ); >> 758 PositDirection.rotateUz(ParticleDirection); >> 759 >> 760 // create G4DynamicParticle object for the particle2 >> 761 aParticle2= new G4DynamicParticle (G4Positron::Positron(), >> 762 PositDirection, PositKineEnergy); >> 763 >> 764 // fill particle change and update initial particle >> 765 aParticleChange.SetNumberOfSecondaries(numberofsecondaries) ; >> 766 if(flagelectron==1) >> 767 aParticleChange.AddSecondary( aParticle1 ) ; >> 768 if(flagpositron==1) >> 769 aParticleChange.AddSecondary( aParticle2 ) ; >> 770 >> 771 G4double NewKinEnergy = KineticEnergy - ElectronEnergy - PositronEnergy ; >> 772 G4double finalMomentum=sqrt(NewKinEnergy* >> 773 (NewKinEnergy+2.*ParticleMass)) ; >> 774 >> 775 aParticleChange.SetMomentumChange( ParticleDirection ); >> 776 >> 777 G4double KinEnergyCut = (aDynamicParticle->GetDefinition()-> >> 778 GetEnergyCuts())[aMaterial->GetIndex()]; >> 779 >> 780 if (NewKinEnergy > KinEnergyCut) >> 781 { >> 782 aParticleChange.SetEnergyChange( NewKinEnergy ); >> 783 } >> 784 else >> 785 { >> 786 aParticleChange.SetEnergyChange(0.); >> 787 LocalEnerDeposit += NewKinEnergy ; >> 788 aParticleChange.SetStatusChange(fStopButAlive); >> 789 } 155 790 156 void G4MuPairProduction::ProcessDescription(st << 791 aParticleChange.SetLocalEnergyDeposit( LocalEnerDeposit ) ; >> 792 >> 793 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData); >> 794 } >> 795 >> 796 G4Element* G4MuPairProduction::SelectRandomAtom(G4Material* aMaterial) const 157 { 797 { 158 out << " Electron-positron pair production << 798 // select randomly 1 element within the material 159 G4VEnergyLossProcess::ProcessDescription(out << 799 >> 800 const G4int Index = aMaterial->GetIndex(); >> 801 const G4int NumberOfElements = aMaterial->GetNumberOfElements(); >> 802 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); >> 803 >> 804 G4double rval = G4UniformRand()*((*PartialSumSigma(Index)) >> 805 (NumberOfElements-1)); >> 806 >> 807 for ( G4int i=0; i < NumberOfElements; i++ ) >> 808 { >> 809 if (rval <= (*PartialSumSigma(Index))(i)) return ((*theElementVector)(i)); >> 810 } >> 811 G4cout << " WARNING !!! - The Material '"<< aMaterial->GetName() >> 812 << "' has no elements, NULL pointer returned." << G4endl; >> 813 return NULL; 160 } 814 } >> 815 void G4MuPairProduction::PrintInfoDefinition() >> 816 { >> 817 G4String comments = "theoretical cross sections \n "; >> 818 comments += " Good description up to 1000 PeV."; >> 819 >> 820 G4cout << G4endl << GetProcessName() << ": " << comments >> 821 << "\n PhysicsTables from " << G4BestUnit(LowerBoundLambda, >> 822 "Energy") >> 823 << " to " << G4BestUnit(UpperBoundLambda,"Energy") >> 824 << " in " << NbinLambda << " bins. \n"; >> 825 } >> 826 161 827 162 //....oooOO0OOooo........oooOO0OOooo........oo << 163 828