| 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: G4MuBremsstrahlung.cc,v 1.27 2003/04/29 04:58:32 kurasige Exp $
>> 25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
28 // 26 //
29 // GEANT4 Class file <<
30 // 27 //
31 // << 28 //--------------- G4MuBremsstrahlung physics process ---------------------------
32 // File name: G4MuBremsstrahlung << 29 // by Laszlo Urban, September 1997
33 // <<
34 // Author: Laszlo Urban <<
35 // <<
36 // Creation date: 30.09.1997 <<
37 // <<
38 // Modifications: <<
39 // 30 //
40 // 08-04-98 remove 'tracking cut' of muon in o 31 // 08-04-98 remove 'tracking cut' of muon in oIt, MMa
41 // 26-10-98 new cross section of R.Kokoulin,cl << 32 // 26/10/98 new cross section of R.Kokoulin,cleanup , L.Urban
42 // 10-02-00 modifications , new e.m. structure << 33 // 10/02/00 modifications , new e.m. structure, L.Urban
43 // 29-05-01 V.Ivanchenko minor changes to prov << 34 // 29/05/01 V.Ivanchenko minor changes to provide ANSI -wall compilation
44 // 09-08-01 new methods Store/Retrieve Physics 35 // 09-08-01 new methods Store/Retrieve PhysicsTable (mma)
45 // 17-09-01 migration of Materials to pure STL 36 // 17-09-01 migration of Materials to pure STL (mma)
46 // 26-09-01 completion of store/retrieve Physi 37 // 26-09-01 completion of store/retrieve PhysicsTable (mma)
47 // 28-09-01 suppression of theMuonPlus ..etc.. 38 // 28-09-01 suppression of theMuonPlus ..etc..data members (mma)
48 // 29-10-01 all static functions no more inlin 39 // 29-10-01 all static functions no more inlined (mma)
49 // 08-11-01 particleMass becomes a local varia 40 // 08-11-01 particleMass becomes a local variable (mma)
50 // 19-08-02 V.Ivanchenko update to new design << 41 // 16-01-03 Migrade to cut per region (V.Ivanchenko)
51 // 23-12-02 Change interface in order to move << 42 // 26-04-03 fix problems of retrieve tables (V.Ivanchenko)
52 // 26-12-02 Secondary production moved to deri << 43 //------------------------------------------------------------------------------
53 // 08-08-03 STD substitute standard (V.Ivanch <<
54 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossPro <<
55 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko) <<
56 // 17-08-04 Utilise mu+ tables for mu- (V.Ivan <<
57 // 08-11-04 Migration to new interface of Stor <<
58 // 08-04-05 Major optimisation of internal int <<
59 // <<
60 // ------------------------------------------- <<
61 // <<
62 //....oooOO0OOooo........oooOO0OOooo........oo <<
63 //....oooOO0OOooo........oooOO0OOooo........oo <<
64 44
65 #include "G4MuBremsstrahlung.hh" 45 #include "G4MuBremsstrahlung.hh"
66 #include "G4SystemOfUnits.hh" << 46 #include "G4UnitsTable.hh"
67 #include "G4Gamma.hh" << 47 #include "G4ProductionCutsTable.hh"
68 #include "G4MuBremsstrahlungModel.hh" << 48
69 #include "G4EmParameters.hh" << 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 50
>> 51 // static members
>> 52 //
>> 53 G4int G4MuBremsstrahlung::nzdat = 5 ;
>> 54 G4double G4MuBremsstrahlung::zdat[]={1.,4.,13.,29.,92.};
>> 55 G4double G4MuBremsstrahlung::adat[]={1.01,9.01,26.98,63.55,238.03};
>> 56 G4int G4MuBremsstrahlung::ntdat = 8 ;
>> 57 G4double G4MuBremsstrahlung::tdat[]={1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
>> 58 G4int G4MuBremsstrahlung::NBIN = 1000; // 100 ;
>> 59 G4double G4MuBremsstrahlung::ya[1001];
>> 60 G4double G4MuBremsstrahlung::proba[5][8][1001];
>> 61 G4double G4MuBremsstrahlung::CutFixed=0.98*keV;
>> 62
>> 63 G4double G4MuBremsstrahlung::LowerBoundLambda = 1.*keV;
>> 64 G4double G4MuBremsstrahlung::UpperBoundLambda = 1000000.*TeV;
>> 65 G4int G4MuBremsstrahlung::NbinLambda = 150;
>> 66
>> 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 68
>> 69 // constructor
>> 70
>> 71 G4MuBremsstrahlung::G4MuBremsstrahlung(const G4String& processName)
>> 72 : G4VMuEnergyLoss(processName),
>> 73 theMeanFreePathTable(NULL)
>> 74 { }
>> 75
>> 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 77
>> 78 G4MuBremsstrahlung::~G4MuBremsstrahlung()
>> 79 {
>> 80 if (theMeanFreePathTable) {
>> 81 theMeanFreePathTable->clearAndDestroy();
>> 82
>> 83 delete theMeanFreePathTable;
>> 84 }
>> 85 PartialSumSigma.clearAndDestroy();
>> 86 }
>> 87
>> 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 89
>> 90 void G4MuBremsstrahlung::SetLowerBoundLambda(G4double val)
>> 91 {LowerBoundLambda = val;}
>> 92
>> 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 94
>> 95 void G4MuBremsstrahlung::SetUpperBoundLambda(G4double val)
>> 96 {UpperBoundLambda = val;}
>> 97
>> 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 99
>> 100 void G4MuBremsstrahlung::SetNbinLambda(G4int n)
>> 101 {NbinLambda = n;}
>> 102
>> 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 104
>> 105 G4double G4MuBremsstrahlung::GetLowerBoundLambda()
>> 106 { return LowerBoundLambda;}
>> 107
>> 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 109
>> 110 G4double G4MuBremsstrahlung::GetUpperBoundLambda()
>> 111 { return UpperBoundLambda;}
>> 112
>> 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 114
>> 115 G4int G4MuBremsstrahlung::GetNbinLambda()
>> 116 {return NbinLambda;}
>> 117
>> 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 119
>> 120 void G4MuBremsstrahlung::BuildPhysicsTable(
>> 121 const G4ParticleDefinition& aParticleType)
>> 122 {
>> 123 if( !CutsWhereModified() && theLossTable) return;
>> 124
>> 125 LowestKineticEnergy = GetLowerBoundEloss() ;
>> 126 HighestKineticEnergy = GetUpperBoundEloss() ;
>> 127 TotBin = GetNbinEloss() ;
>> 128
>> 129 BuildLossTable(aParticleType) ;
>> 130
>> 131 if(&aParticleType==G4MuonMinus::MuonMinus())
>> 132 {
>> 133 RecorderOfmuminusProcess[CounterOfmuminusProcess] = (*this).theLossTable;
>> 134 CounterOfmuminusProcess++;
>> 135 }
>> 136 else
>> 137 {
>> 138 RecorderOfmuplusProcess[CounterOfmuplusProcess] = (*this).theLossTable;
>> 139 CounterOfmuplusProcess++;
>> 140 }
>> 141
>> 142 if( !theMeanFreePathTable ) MakeSamplingTables(&aParticleType) ;
>> 143
>> 144 BuildLambdaTable(aParticleType) ;
>> 145
>> 146 G4VMuEnergyLoss::BuildDEDXTable(aParticleType);
>> 147
>> 148 if(&aParticleType == G4MuonPlus::MuonPlus()) PrintInfoDefinition();
>> 149 }
>> 150
>> 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 152
>> 153 void G4MuBremsstrahlung::BuildLossTable(
>> 154 const G4ParticleDefinition& aParticleType)
>> 155 {
>> 156 G4double KineticEnergy,TotalEnergy,bremloss,Z,
>> 157 loss,natom ;
>> 158
>> 159 const G4ProductionCutsTable* theCoupleTable=
>> 160 G4ProductionCutsTable::GetProductionCutsTable();
>> 161 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 162
>> 163 if (theLossTable) {theLossTable->clearAndDestroy(); delete theLossTable;}
>> 164 theLossTable = new G4PhysicsTable(numOfCouples);
>> 165
>> 166 secondaryEnergyCuts = theCoupleTable->GetEnergyCutsVector(0);
>> 167 G4double particleMass = aParticleType.GetPDGMass();
>> 168
>> 169 // loop for materials
>> 170 //
>> 171 for (size_t J=0; J<numOfCouples; J++)
>> 172 {
>> 173
>> 174 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
>> 175 LowestKineticEnergy,HighestKineticEnergy,TotBin);
>> 176
>> 177 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(J);
>> 178 const G4Material* material= couple->GetMaterial();
>> 179
>> 180 G4double Cut = SecondaryEnergyThreshold(J);
>> 181
>> 182 const G4ElementVector* theElementVector =
>> 183 material->GetElementVector() ;
>> 184 const G4double* theAtomicNumDensityVector =
>> 185 material->GetAtomicNumDensityVector() ;
>> 186 const G4int NumberOfElements = material->GetNumberOfElements() ;
>> 187
>> 188 for (G4int i=0; i<TotBin; i++)
>> 189 {
>> 190 KineticEnergy = aVector->GetLowEdgeEnergy(i) ;
>> 191 TotalEnergy = KineticEnergy+particleMass ;
>> 192
>> 193 if(Cut>KineticEnergy) Cut = KineticEnergy ;
>> 194 bremloss = 0.;
>> 195
>> 196 for (G4int iel=0; iel<NumberOfElements; iel++)
>> 197 {
>> 198 Z=(*theElementVector)[iel]->GetZ();
>> 199 natom = theAtomicNumDensityVector[iel] ;
>> 200 loss = ComputeBremLoss((&aParticleType),Z,
>> 201 (*theElementVector)[iel]->GetA(),
>> 202 KineticEnergy,Cut) ;
>> 203 bremloss += natom*loss ;
>> 204 }
>> 205 if(bremloss<0.) bremloss = 0. ;
>> 206
>> 207 aVector->PutValue(i,bremloss);
>> 208 }
>> 209 theLossTable->insert(aVector);
>> 210 }
>> 211 }
70 212
71 //....oooOO0OOooo........oooOO0OOooo........oo << 213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 214
73 using namespace std; << 215 G4double G4MuBremsstrahlung::ComputeBremLoss(
>> 216 const G4ParticleDefinition* aParticleType,
>> 217 G4double AtomicNumber,G4double AtomicMass,
>> 218 G4double KineticEnergy,G4double GammaEnergyCut)
>> 219 {
>> 220 G4double TotalEnergy,vcut,vmax,aaa,bbb,hhh,aa,x,ep ;
>> 221 G4int kkk ;
>> 222 G4double ak1=0.05 ;
>> 223 G4int k2=5 ;
>> 224 G4double xgi[]={0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
>> 225 G4double wgi[]={0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
>> 226 G4double loss = 0. ;
>> 227
>> 228 G4double particleMass = aParticleType->GetPDGMass();
>> 229 TotalEnergy=KineticEnergy+particleMass ;
>> 230 vcut = GammaEnergyCut/TotalEnergy ;
>> 231 vmax = KineticEnergy/TotalEnergy ;
>> 232
>> 233 aaa=0.;
>> 234 bbb=vcut ;
>> 235 if(vcut>vmax) bbb=vmax ;
>> 236 kkk=int((bbb-aaa)/ak1)+k2 ;
>> 237 hhh=(bbb-aaa)/float(kkk) ;
>> 238
>> 239 for(G4int l=0; l<kkk; l++)
>> 240 {
>> 241 aa=aaa+hhh*float(l) ;
>> 242 for(G4int i=0; i<6; i++)
>> 243 {
>> 244 x=aa+xgi[i]*hhh ;
>> 245 ep=x*TotalEnergy ;
>> 246 loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(
>> 247 aParticleType,KineticEnergy,
>> 248 AtomicNumber,AtomicMass,ep) ;
>> 249 }
>> 250 }
>> 251
>> 252 loss *=hhh*TotalEnergy ;
>> 253
>> 254 return loss ;
>> 255 }
>> 256
>> 257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 258
>> 259 void G4MuBremsstrahlung::BuildLambdaTable(
>> 260 const G4ParticleDefinition& ParticleType)
>> 261 {
>> 262
>> 263 G4double LowEdgeEnergy , Value;
>> 264 G4double FixedEnergy = (LowestKineticEnergy + HighestKineticEnergy)/2. ;
>> 265
>> 266 //create table
>> 267 //
>> 268 const G4ProductionCutsTable* theCoupleTable=
>> 269 G4ProductionCutsTable::GetProductionCutsTable();
>> 270 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 271
>> 272 //create table
>> 273 if (theMeanFreePathTable) {theMeanFreePathTable->clearAndDestroy();
>> 274 delete theMeanFreePathTable;
>> 275 }
>> 276 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
>> 277
>> 278 PartialSumSigma.clearAndDestroy();
>> 279 PartialSumSigma.resize(numOfCouples);
>> 280
>> 281 G4PhysicsLogVector* ptrVector;
>> 282 for ( size_t J=0; J<numOfCouples; J++ )
>> 283 {
>> 284 ptrVector = new G4PhysicsLogVector(
>> 285 LowerBoundLambda,UpperBoundLambda,NbinLambda);
>> 286
>> 287 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(J);
>> 288
>> 289 for ( G4int i = 0 ; i < NbinLambda ; i++ )
>> 290 {
>> 291 LowEdgeEnergy = ptrVector->GetLowEdgeEnergy( i ) ;
>> 292 Value = ComputeMeanFreePath( &ParticleType, LowEdgeEnergy, couple);
>> 293 ptrVector->PutValue( i , Value ) ;
>> 294 }
>> 295
>> 296 theMeanFreePathTable->insertAt( J , ptrVector );
>> 297
>> 298 // Compute the PartialSumSigma table at a given fixed energy
>> 299 ComputePartialSumSigma( &ParticleType, FixedEnergy, couple) ;
>> 300 }
>> 301 }
>> 302
>> 303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 304
>> 305 void G4MuBremsstrahlung::ComputePartialSumSigma(
>> 306 const G4ParticleDefinition* ParticleType,
>> 307 G4double KineticEnergy,
>> 308 const G4MaterialCutsCouple* couple)
>> 309 // Build the table of cross section per element.
>> 310 // The table is built for MATERIALS.
>> 311 // This table is used by DoIt to select randomly an element in the material.
>> 312 {
>> 313 const G4Material* aMaterial = couple->GetMaterial();
>> 314 size_t index = couple->GetIndex();
>> 315 G4int NbOfElements = aMaterial->GetNumberOfElements();
>> 316 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 317 const G4double* theAtomNumDensityVector =
>> 318 aMaterial->GetAtomicNumDensityVector();
>> 319 G4double GammaEnergyCut = SecondaryEnergyThreshold(index);
>> 320
>> 321 PartialSumSigma[index] = new G4DataVector();
>> 322
>> 323 G4double SIGMA = 0. ;
>> 324
>> 325 for ( G4int Ielem=0 ; Ielem < NbOfElements ; Ielem++ )
>> 326 {
>> 327 SIGMA += theAtomNumDensityVector[Ielem] *
>> 328 ComputeMicroscopicCrossSection( ParticleType, KineticEnergy,
>> 329 (*theElementVector)[Ielem]->GetZ(),
>> 330 (*theElementVector)[Ielem]->GetA(),
>> 331 GammaEnergyCut );
>> 332 PartialSumSigma[index]->push_back(SIGMA);
>> 333 }
>> 334 }
>> 335
>> 336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 337
>> 338 G4double G4MuBremsstrahlung::ComputeMicroscopicCrossSection(
>> 339 const G4ParticleDefinition* ParticleType,
>> 340 G4double KineticEnergy,
>> 341 G4double AtomicNumber,
>> 342 G4double AtomicMass,
>> 343 G4double GammaEnergyCut)
>> 344 // Cross section is calculated according to a formula of R.Kokoulin.
>> 345 {
>> 346 G4double TotalEnergy,vcut,vmax,aaa,bbb,hhh,aa,x,ep ;
>> 347 G4int kkk ;
>> 348 G4double ak1=2.3 ;
>> 349 G4int k2=4 ;
>> 350 G4double xgi[]={0.03377,0.16940,0.38069,0.61931,0.83060,0.96623};
>> 351 G4double wgi[]={0.08566,0.18038,0.23396,0.23396,0.18038,0.08566};
>> 352 G4double CrossSection = 0. ;
>> 353
>> 354 G4double particleMass = ParticleType->GetPDGMass();
>> 355 TotalEnergy=KineticEnergy+particleMass ;
>> 356 vcut = GammaEnergyCut/TotalEnergy ;
>> 357 vmax = KineticEnergy/TotalEnergy ;
>> 358 if(vmax <= vcut) return CrossSection;
>> 359
>> 360 // numerical integration
>> 361 aaa=log(vcut) ;
>> 362 bbb=log(vmax);
>> 363 kkk=int((bbb-aaa)/ak1)+k2 ;
>> 364 hhh=(bbb-aaa)/float(kkk) ;
>> 365
>> 366 for(G4int l=0; l<kkk; l++)
>> 367 {
>> 368 aa=aaa+hhh*float(l) ;
>> 369 for(G4int i=0; i<6; i++)
>> 370 {
>> 371 x=aa+xgi[i]*hhh ;
>> 372 ep=exp(x)*TotalEnergy ;
>> 373 CrossSection += ep*wgi[i]*ComputeDMicroscopicCrossSection(
>> 374 ParticleType,KineticEnergy,
>> 375 AtomicNumber,AtomicMass,ep) ;
>> 376 }
>> 377 }
>> 378
>> 379 CrossSection *= hhh ;
>> 380
>> 381 return CrossSection;
>> 382 }
>> 383
>> 384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 385
>> 386 G4double G4MuBremsstrahlung::GetDMicroscopicCrossSection(
>> 387 const G4ParticleDefinition* ParticleType,
>> 388 G4double KineticEnergy,
>> 389 G4double AtomicNumber,
>> 390 G4double AtomicMass,
>> 391 G4double GammaEnergy)
>> 392 // get differential cross section
>> 393 {
>> 394 return ComputeDMicroscopicCrossSection(ParticleType,KineticEnergy,
>> 395 AtomicNumber,AtomicMass,GammaEnergy);
>> 396 }
>> 397
>> 398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74 399
75 G4MuBremsstrahlung::G4MuBremsstrahlung(const G << 400 G4double G4MuBremsstrahlung::ComputeDMicroscopicCrossSection(
76 : G4VEnergyLossProcess(name), << 401 const G4ParticleDefinition* ParticleType,
77 lowestKinEnergy(0.1*CLHEP::GeV) << 402 G4double KineticEnergy,
>> 403 G4double AtomicNumber,
>> 404 G4double AtomicMass,
>> 405 G4double GammaEnergy)
>> 406 // differential cross section
78 { 407 {
79 SetProcessSubType(fBremsstrahlung); << 408 G4double particleMass = ParticleType->GetPDGMass();
80 SetSecondaryParticle(G4Gamma::Gamma()); << 409
81 SetIonisation(false); << 410 static const G4double sqrte=sqrt(exp(1.)) ;
>> 411 static const G4double bh=202.4,bh1=446.,btf=183.,btf1=1429. ;
>> 412 static const G4double rmass=particleMass/electron_mass_c2 ;
>> 413 static const G4double cc=classic_electr_radius/rmass ;
>> 414 static const G4double coeff= 16.*fine_structure_const*cc*cc/3. ;
>> 415
>> 416 G4double dxsection = 0.;
>> 417
>> 418 if( GammaEnergy > KineticEnergy) return dxsection ;
>> 419
>> 420 G4double A = AtomicMass/(g/mole) ; // !!!!!!!!!!!!!!!!!!!
>> 421 G4double E=KineticEnergy+particleMass ;
>> 422 G4double v=GammaEnergy/E ;
>> 423 G4double delta=0.5*particleMass*particleMass*v/(E-GammaEnergy) ;
>> 424 G4double rab0=delta*sqrte ;
>> 425
>> 426 G4double z13=exp(-log(AtomicNumber)/3.) ;
>> 427 G4double dn=1.54*exp(0.27*log(A)) ;
>> 428
>> 429 G4double b,b1,dnstar ;
>> 430
>> 431 if(AtomicNumber<1.5)
>> 432 {
>> 433 b=bh;
>> 434 b1=bh1;
>> 435 dnstar=dn ;
>> 436 }
>> 437 else
>> 438 {
>> 439 b=btf;
>> 440 b1=btf1;
>> 441 dnstar = exp((1.-1./AtomicNumber)*log(dn)) ;
>> 442 }
>> 443
>> 444 // nucleus contribution logarithm
>> 445 G4double rab1=b*z13;
>> 446 G4double fn=log(rab1/(dnstar*(electron_mass_c2+rab0*rab1))*
>> 447 (particleMass+delta*(dnstar*sqrte-2.))) ;
>> 448 if(fn <0.) fn = 0. ;
>> 449 // electron contribution logarithm
>> 450 G4double epmax1=E/(1.+0.5*particleMass*rmass/E) ;
>> 451 G4double fe=0.;
>> 452 if(GammaEnergy<epmax1)
>> 453 {
>> 454 G4double rab2=b1*z13*z13 ;
>> 455 fe=log(rab2*particleMass/((1.+delta*rmass/(electron_mass_c2*sqrte))*
>> 456 (electron_mass_c2+rab0*rab2))) ;
>> 457 if(fe<0.) fe=0. ;
>> 458 }
>> 459
>> 460 dxsection = coeff*(1.-v*(1.-0.75*v))*AtomicNumber*(fn*AtomicNumber+fe)/
>> 461 GammaEnergy ;
>> 462 return dxsection ;
82 } 463 }
83 464
84 //....oooOO0OOooo........oooOO0OOooo........oo <<
85 465
86 G4bool G4MuBremsstrahlung::IsApplicable(const << 466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 467
>> 468 void G4MuBremsstrahlung::MakeSamplingTables(
>> 469 const G4ParticleDefinition* ParticleType)
87 { 470 {
88 return (p.GetPDGCharge() != 0.0); << 471 G4int nbin;
>> 472 G4double AtomicNumber,AtomicWeight,KineticEnergy,
>> 473 TotalEnergy,Maxep ;
>> 474
>> 475 G4double particleMass = ParticleType->GetPDGMass() ;
>> 476
>> 477 for (G4int iz=0; iz<nzdat; iz++)
>> 478 {
>> 479 AtomicNumber = zdat[iz];
>> 480 AtomicWeight = adat[iz]*g/mole ;
>> 481
>> 482 for (G4int it=0; it<ntdat; it++)
>> 483 {
>> 484 KineticEnergy = tdat[it];
>> 485 TotalEnergy = KineticEnergy + particleMass;
>> 486 Maxep = KineticEnergy ;
>> 487
>> 488 G4double CrossSection = 0.0 ;
>> 489
>> 490 G4double c,y,ymin,ymax,dy,yy,dx,x,ep;
>> 491
>> 492 //G4int NbofIntervals ;
>> 493 // calculate the differential cross section
>> 494 // numerical integration in
>> 495 // log ...............
>> 496 c = log(Maxep/CutFixed) ;
>> 497 ymin = -5. ;
>> 498 ymax = 0. ;
>> 499 dy = (ymax-ymin)/NBIN ;
>> 500 nbin=-1;
>> 501
>> 502 y = ymin - 0.5*dy ;
>> 503 yy = ymin - dy ;
>> 504 for (G4int i=0 ; i<NBIN; i++)
>> 505 {
>> 506 y += dy ;
>> 507 x = exp(y) ;
>> 508 yy += dy ;
>> 509 dx = exp(yy+dy)-exp(yy) ;
>> 510
>> 511 ep = CutFixed*exp(c*x) ;
>> 512
>> 513 CrossSection += ep*dx*ComputeDMicroscopicCrossSection(ParticleType,
>> 514 KineticEnergy,AtomicNumber,
>> 515 AtomicWeight,ep) ;
>> 516 if(nbin<NBIN)
>> 517 {
>> 518 nbin += 1 ;
>> 519 ya[nbin]=y ;
>> 520 proba[iz][it][nbin] = CrossSection ;
>> 521 }
>> 522 }
>> 523
>> 524 ya[NBIN] = 0. ; // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> 525
>> 526 if(CrossSection > 0.)
>> 527 {
>> 528 for(G4int ib=0; ib<=nbin; ib++)
>> 529 {
>> 530 proba[iz][it][ib] /= CrossSection ;
>> 531 }
>> 532 }
>> 533 }
>> 534 }
89 } 535 }
90 536
91 //....oooOO0OOooo........oooOO0OOooo........oo << 537 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 538
92 539
93 G4double G4MuBremsstrahlung::MinPrimaryEnergy( << 540 G4VParticleChange* G4MuBremsstrahlung::PostStepDoIt(const G4Track& trackData,
94 const G4Material*, << 541 const G4Step& stepData)
95 G4double) <<
96 { 542 {
97 return lowestKinEnergy; << 543
>> 544 static G4double ysmall = -100. ;
>> 545 static G4double ytablelow = -5. ;
>> 546
>> 547 aParticleChange.Initialize(trackData);
>> 548 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
>> 549
>> 550 const G4DynamicParticle* aDynamicParticle=trackData.GetDynamicParticle();
>> 551
>> 552 G4double KineticEnergy = aDynamicParticle->GetKineticEnergy();
>> 553 G4ParticleMomentum ParticleDirection =
>> 554 aDynamicParticle->GetMomentumDirection();
>> 555
>> 556 // Gamma cut in this material
>> 557 G4double GammaEnergyCut = SecondaryEnergyThreshold(couple->GetIndex());
>> 558
>> 559 // check against insufficient energy
>> 560 if(KineticEnergy < GammaEnergyCut)
>> 561 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 562
>> 563 // select randomly one element constituing the material
>> 564 const G4Element* anElement = SelectRandomAtom(couple);
>> 565
>> 566 G4double TotalEnergy=KineticEnergy+aDynamicParticle->
>> 567 GetDefinition()->GetPDGMass() ;
>> 568
>> 569 G4double dy = 5./G4float(NBIN) ;
>> 570
>> 571 G4double ymin=log(log(GammaEnergyCut/CutFixed)/log(KineticEnergy/CutFixed)) ;
>> 572
>> 573 if(ymin < ysmall)
>> 574 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 575
>> 576 // sampling using tables
>> 577 //G4double v,xc,x,yc,y ;
>> 578 //G4int iZ,iT,iy ;
>> 579 G4double v,x,y ;
>> 580 G4int iy;
>> 581 // select sampling table ;
>> 582 G4double lnZ = log(anElement->GetZ()) ;
>> 583 G4double delmin = 1.e10 ;
>> 584 G4double del ;
>> 585 G4int izz = 0;
>> 586 G4int itt = 0;
>> 587 G4int NBINminus1;
>> 588 NBINminus1 = NBIN-1 ;
>> 589 for (G4int iz=0; iz<nzdat; iz++)
>> 590 {
>> 591 del = abs(lnZ-log(zdat[iz])) ;
>> 592 if(del<delmin)
>> 593 {
>> 594 delmin=del ;
>> 595 izz=iz ;
>> 596 }
>> 597 }
>> 598
>> 599 delmin = 1.e10 ;
>> 600 for (G4int it=0; it<ntdat; it++)
>> 601 {
>> 602 del = abs(log(KineticEnergy)-log(tdat[it])) ;
>> 603 if(del<delmin)
>> 604 {
>> 605 delmin=del;
>> 606 itt=it ;
>> 607 }
>> 608 }
>> 609 G4int iymin = G4int((ymin+5.)/dy+0.5) ;
>> 610
>> 611 if(ymin < ytablelow)
>> 612 {
>> 613 y = ymin + G4UniformRand()*(ytablelow-ymin) ;
>> 614 }
>> 615 else
>> 616 {
>> 617 G4double r = G4UniformRand() ;
>> 618
>> 619 iy = iymin-1 ;
>> 620 delmin = proba[izz][itt][NBINminus1]-proba[izz][itt][iymin] ;
>> 621 do {
>> 622 iy += 1 ;
>> 623 } while ((r > (proba[izz][itt][iy]-proba[izz][itt][iymin])/delmin)
>> 624 &&(iy < NBINminus1)) ;
>> 625
>> 626 //sampling is Done uniformly in y in the bin
>> 627 y = ya[iy] + G4UniformRand() * ( ya[iy+1] - ya[iy] ) ;
>> 628 }
>> 629
>> 630 x = exp(y) ;
>> 631
>> 632 v = CutFixed*exp(x*log(KineticEnergy/CutFixed)) ;
>> 633 if( v <= 0.)
>> 634 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
>> 635
>> 636 // create G4DynamicParticle object for the Gamma
>> 637 G4double GammaEnergy = v;
>> 638
>> 639 // angles of the emitted gamma. ( Z - axis along the parent particle)
>> 640 // Teta = electron_mass_c2/TotalEnergy for the moment .....
>> 641
>> 642 G4double Teta = electron_mass_c2/TotalEnergy ;
>> 643 G4double Phi = twopi * G4UniformRand() ;
>> 644 G4double dirx = sin(Teta)*cos(Phi) , diry = sin(Teta)*sin(Phi) ,
>> 645 dirz = cos(Teta) ;
>> 646
>> 647 G4ThreeVector GammaDirection ( dirx, diry, dirz);
>> 648 GammaDirection.rotateUz(ParticleDirection);
>> 649
>> 650 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
>> 651 GammaDirection, GammaEnergy);
>> 652
>> 653 aParticleChange.SetNumberOfSecondaries(1);
>> 654 aParticleChange.AddSecondary(aGamma);
>> 655
>> 656 // Update the incident particle
>> 657 G4double NewKinEnergy = KineticEnergy - GammaEnergy;
>> 658 if (NewKinEnergy > 0.)
>> 659 {
>> 660 aParticleChange.SetMomentumChange(ParticleDirection);
>> 661 aParticleChange.SetEnergyChange(NewKinEnergy);
>> 662 aParticleChange.SetLocalEnergyDeposit (0.);
>> 663 }
>> 664 else
>> 665 {
>> 666 aParticleChange.SetEnergyChange(0.);
>> 667 aParticleChange.SetLocalEnergyDeposit (0.);
>> 668 aParticleChange.SetStatusChange(fStopButAlive);
>> 669 }
>> 670
>> 671 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
98 } 672 }
99 673
100 //....oooOO0OOooo........oooOO0OOooo........oo << 674 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 675
102 void G4MuBremsstrahlung::InitialiseEnergyLossP << 676 const G4Element* G4MuBremsstrahlung::SelectRandomAtom(
103 const G4ParticleDefinition*, << 677 const G4MaterialCutsCouple* couple) const
104 const G4ParticleDefinition*) <<
105 { 678 {
106 if(isInitialised) { return; } << 679 // select randomly 1 element within the material
107 isInitialised = true; << 680 size_t index = couple->GetIndex();
>> 681 const G4Material* aMaterial = couple->GetMaterial();
>> 682 const G4int NumberOfElements = aMaterial->GetNumberOfElements();
>> 683 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
>> 684
>> 685 G4double rval = G4UniformRand()
>> 686 *((*PartialSumSigma[index])[NumberOfElements-1]);
>> 687 for ( G4int i=0; i < NumberOfElements; i++ )
>> 688 if (rval <= (*PartialSumSigma[index])[i]) return ((*theElementVector)[i]);
>> 689 G4cout << " WARNING !!! - The Material " << aMaterial->GetName()
>> 690 << " has no elements, NULL pointer returned." << G4endl;
>> 691 return NULL;
>> 692 }
108 693
109 if (nullptr == EmModel(0)) { SetEmModel(new << 694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 695
111 G4VEmFluctuationModel* fm = nullptr; << 696 G4bool G4MuBremsstrahlung::StorePhysicsTable(G4ParticleDefinition* particle,
112 G4EmParameters* param = G4EmParameters::Inst << 697 const G4String& directory,
113 EmModel(0)->SetLowEnergyLimit(param->MinKinE << 698 G4bool ascii)
114 EmModel(0)->SetHighEnergyLimit(param->MaxKin << 699 {
115 EmModel(0)->SetSecondaryThreshold(param->MuH << 700 G4String filename;
116 AddEmModel(1, EmModel(0), fm); << 701
>> 702 // store stopping power table
>> 703 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 704 if ( !theLossTable->StorePhysicsTable(filename, ascii) ){
>> 705 G4cout << " FAIL theLossTable->StorePhysicsTable in " << filename
>> 706 << G4endl;
>> 707 return false;
>> 708 }
>> 709 // store mean free path table
>> 710 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 711 if ( !theMeanFreePathTable->StorePhysicsTable(filename, ascii) ){
>> 712 G4cout << " FAIL theMeanFreePathTable->StorePhysicsTable in " << filename
>> 713 << G4endl;
>> 714 return false;
>> 715 }
>> 716
>> 717 // store PartialSumSigma table (G4OrderedTable)
>> 718 filename = GetPhysicsTableFileName(particle,directory,"PartSumSigma",ascii);
>> 719 if ( !PartialSumSigma.Store(filename, ascii) ){
>> 720 G4cout << " FAIL PartialSumSigma.store in " << filename
>> 721 << G4endl;
>> 722 return false;
>> 723 }
>> 724
>> 725 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 726 << ": Success to store the PhysicsTables in "
>> 727 << directory << G4endl;
>> 728
>> 729 return true;
117 } 730 }
118 731
119 //....oooOO0OOooo........oooOO0OOooo........oo << 732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 733
121 void G4MuBremsstrahlung::ProcessDescription(st << 734 G4bool G4MuBremsstrahlung::RetrievePhysicsTable(G4ParticleDefinition* particle,
>> 735 const G4String& directory,
>> 736 G4bool ascii)
122 { 737 {
123 out << " Muon bremsstrahlung"; << 738 // delete theLossTable and theMeanFreePathTable
124 G4VEnergyLossProcess::ProcessDescription(out << 739 if (theLossTable != 0) {
>> 740 theLossTable->clearAndDestroy();
>> 741 delete theLossTable;
>> 742 }
>> 743 if (theMeanFreePathTable != 0) {
>> 744 theMeanFreePathTable->clearAndDestroy();
>> 745 delete theMeanFreePathTable;
>> 746 }
>> 747
>> 748 // get bining from EnergyLoss
>> 749 LowestKineticEnergy = GetLowerBoundEloss();
>> 750 HighestKineticEnergy = GetUpperBoundEloss();
>> 751 TotBin = GetNbinEloss();
>> 752
>> 753 G4String filename;
>> 754 const G4ProductionCutsTable* theCoupleTable=
>> 755 G4ProductionCutsTable::GetProductionCutsTable();
>> 756 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 757
>> 758 secondaryEnergyCuts = theCoupleTable->GetEnergyCutsVector(0);
>> 759 PartialSumSigma.clearAndDestroy();
>> 760 PartialSumSigma.reserve(numOfCouples);
>> 761
>> 762 // retreive stopping power table
>> 763 filename = GetPhysicsTableFileName(particle,directory,"StoppingPower",ascii);
>> 764 theLossTable = new G4PhysicsTable(numOfCouples);
>> 765 if ( !theLossTable->RetrievePhysicsTable(filename, ascii) ){
>> 766 G4cout << " FAIL theLossTable->RetrievePhysicsTable in " << filename
>> 767 << G4endl;
>> 768 return false;
>> 769 }
>> 770
>> 771 // retreive mean free path table
>> 772 filename = GetPhysicsTableFileName(particle,directory,"MeanFreePath",ascii);
>> 773 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
>> 774 if ( !theMeanFreePathTable->RetrievePhysicsTable(filename, ascii) ){
>> 775 G4cout << " FAIL theMeanFreePathTable->RetrievePhysicsTable in " << filename
>> 776 << G4endl;
>> 777 return false;
>> 778 }
>> 779
>> 780 // retrieve PartialSumSigma table (G4OrderedTable)
>> 781 filename = GetPhysicsTableFileName(particle,directory,"PartSumSigma",ascii);
>> 782 if ( !PartialSumSigma.Retrieve(filename, ascii) ){
>> 783 G4cout << " FAIL PartialSumSigma.retrieve in " << filename
>> 784 << G4endl;
>> 785 return false;
>> 786 }
>> 787
>> 788 G4cout << GetProcessName() << " for " << particle->GetParticleName()
>> 789 << ": Success to retrieve the PhysicsTables from "
>> 790 << directory << G4endl;
>> 791
>> 792 if (particle->GetPDGCharge() < 0.)
>> 793 {
>> 794 RecorderOfmuminusProcess[CounterOfmuminusProcess] = (*this).theLossTable;
>> 795 CounterOfmuminusProcess++;
>> 796 }
>> 797 else
>> 798 {
>> 799 RecorderOfmuplusProcess[CounterOfmuplusProcess] = (*this).theLossTable;
>> 800 CounterOfmuplusProcess++;
>> 801 }
>> 802
>> 803 secondaryEnergyCuts = theCoupleTable->GetEnergyCutsVector(0);
>> 804 MakeSamplingTables(particle);
>> 805
>> 806 G4VMuEnergyLoss::BuildDEDXTable(*particle);
>> 807 if(particle==G4MuonPlus::MuonPlus()) PrintInfoDefinition();
>> 808
>> 809 return true;
125 } 810 }
126 811
127 //....oooOO0OOooo........oooOO0OOooo........oo << 812 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 813
>> 814 void G4MuBremsstrahlung::PrintInfoDefinition()
>> 815 {
>> 816 G4String comments = "theoretical cross section \n ";
>> 817 comments += " Good description up to 1000 PeV.";
>> 818
>> 819 G4cout << G4endl << GetProcessName() << ": " << comments
>> 820 << "\n PhysicsTables from " << G4BestUnit(LowerBoundLambda,
>> 821 "Energy")
>> 822 << " to " << G4BestUnit(UpperBoundLambda,"Energy")
>> 823 << " in " << NbinLambda << " bins. \n";
>> 824 }
>> 825
>> 826 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 827
128 828