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