Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc (Version 11.3.0) and /processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc (Version 5.2)


  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 $
 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