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