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