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