Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/muons/src/G4MuPairProductionModel.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/G4MuPairProductionModel.cc (Version 11.3.0) and /processes/electromagnetic/muons/src/G4MuPairProductionModel.cc (Version 6.2.p2)


  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 //
                                                   >>  23 // $Id: G4MuPairProductionModel.cc,v 1.19 2004/05/05 18:45:42 vnivanch Exp $
                                                   >>  24 // GEANT4 tag $Name: geant4-06-02 $
 26 //                                                 25 //
 27 // -------------------------------------------     26 // -------------------------------------------------------------------
 28 //                                                 27 //
 29 // GEANT4 Class file                               28 // GEANT4 Class file
 30 //                                                 29 //
 31 //                                                 30 //
 32 // File name:     G4MuPairProductionModel          31 // File name:     G4MuPairProductionModel
 33 //                                                 32 //
 34 // Author:        Vladimir Ivanchenko on base      33 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 35 //                                                 34 //
 36 // Creation date: 24.06.2002                       35 // Creation date: 24.06.2002
 37 //                                                 36 //
 38 // Modifications:                                  37 // Modifications:
 39 //                                                 38 //
 40 // 04-12-02 Change G4DynamicParticle construct     39 // 04-12-02 Change G4DynamicParticle constructor in PostStep (V.Ivanchenko)
 41 // 23-12-02 Change interface in order to move      40 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
 42 // 24-01-03 Fix for compounds (V.Ivanchenko)       41 // 24-01-03 Fix for compounds (V.Ivanchenko)
 43 // 27-01-03 Make models region aware (V.Ivanch     42 // 27-01-03 Make models region aware (V.Ivanchenko)
 44 // 13-02-03 Add model (V.Ivanchenko)               43 // 13-02-03 Add model (V.Ivanchenko)
 45 // 06-06-03 Fix in cross section calculation f     44 // 06-06-03 Fix in cross section calculation for high energy (V.Ivanchenko)
 46 // 20-10-03 2*xi in ComputeDDMicroscopicCrossS     45 // 20-10-03 2*xi in ComputeDDMicroscopicCrossSection   (R.Kokoulin)
 47 //          8 integration points in ComputeDMi     46 //          8 integration points in ComputeDMicroscopicCrossSection
 48 // 12-01-04 Take min cut of e- and e+ not its      47 // 12-01-04 Take min cut of e- and e+ not its sum (V.Ivanchenko)
 49 // 10-02-04 Update parameterisation using R.Ko     48 // 10-02-04 Update parameterisation using R.Kokoulin model (V.Ivanchenko)
 50 // 28-04-04 For complex materials repeat calcu     49 // 28-04-04 For complex materials repeat calculation of max energy for each
 51 //          material (V.Ivanchenko)                50 //          material (V.Ivanchenko)
 52 // 01-11-04 Fix bug inside ComputeDMicroscopic << 
 53 // 08-04-05 Major optimisation of internal int << 
 54 // 03-08-05 Add SetParticle method (V.Ivantche << 
 55 // 23-10-05 Add protection in sampling of e+e- << 
 56 //          low cuts (V.Ivantchenko)           << 
 57 // 13-02-06 Add ComputeCrossSectionPerAtom (mm << 
 58 // 24-04-07 Add protection in SelectRandomAtom << 
 59 // 12-05-06 Updated sampling (use cut) in Sele << 
 60 // 11-10-07 Add ignoreCut flag (V.Ivanchenko)  << 
 61                                                    51 
 62 //                                                 52 //
 63 // Class Description:                              53 // Class Description:
 64 //                                                 54 //
 65 //                                                 55 //
 66 // -------------------------------------------     56 // -------------------------------------------------------------------
 67 //                                                 57 //
 68 //....oooOO0OOooo........oooOO0OOooo........oo     58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 69 //....oooOO0OOooo........oooOO0OOooo........oo     59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 70                                                    60 
 71 #include "G4MuPairProductionModel.hh"              61 #include "G4MuPairProductionModel.hh"
 72 #include "G4PhysicalConstants.hh"              << 
 73 #include "G4SystemOfUnits.hh"                  << 
 74 #include "G4EmParameters.hh"                   << 
 75 #include "G4Electron.hh"                           62 #include "G4Electron.hh"
 76 #include "G4Positron.hh"                           63 #include "G4Positron.hh"
 77 #include "G4MuonMinus.hh"                          64 #include "G4MuonMinus.hh"
 78 #include "G4MuonPlus.hh"                           65 #include "G4MuonPlus.hh"
 79 #include "Randomize.hh"                            66 #include "Randomize.hh"
 80 #include "G4Material.hh"                           67 #include "G4Material.hh"
 81 #include "G4Element.hh"                            68 #include "G4Element.hh"
 82 #include "G4ElementVector.hh"                      69 #include "G4ElementVector.hh"
 83 #include "G4ElementDataRegistry.hh"            << 
 84 #include "G4ProductionCutsTable.hh"                70 #include "G4ProductionCutsTable.hh"
 85 #include "G4ParticleChangeForLoss.hh"          << 
 86 #include "G4ModifiedMephi.hh"                  << 
 87 #include "G4Log.hh"                            << 
 88 #include "G4Exp.hh"                            << 
 89 #include "G4AutoLock.hh"                       << 
 90                                                << 
 91 #include <iostream>                            << 
 92 #include <fstream>                             << 
 93                                                    71 
 94 //....oooOO0OOooo........oooOO0OOooo........oo     72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 95                                                    73 
 96 const G4int G4MuPairProductionModel::ZDATPAIR[ <<  74 // static members
 97                                                <<  75 //
 98 const G4double G4MuPairProductionModel::xgi[]  <<  76 G4double G4MuPairProductionModel::zdat[]={1.,4.,13.,29.,92.};
 99     0.0198550717512320, 0.1016667612931865, 0. <<  77 G4double G4MuPairProductionModel::adat[]={1.01,9.01,26.98,63.55,238.03};
100     0.5917173212478250, 0.7627662049581645, 0. <<  78 G4double G4MuPairProductionModel::tdat[]={1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
101   };                                           <<  79 G4double G4MuPairProductionModel::xgi[]={ 0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801 };
102                                                <<  80 G4double G4MuPairProductionModel::wgi[]={ 0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506 };
103 const G4double G4MuPairProductionModel::wgi[]  << 
104     0.0506142681451880, 0.1111905172266870, 0. << 
105     0.1813418916891810, 0.1568533229389435, 0. << 
106   };                                           << 
107                                                << 
108 namespace                                      << 
109 {                                              << 
110   G4Mutex theMuPairMutex = G4MUTEX_INITIALIZER << 
111                                                << 
112   const G4double ak1 = 6.9;                    << 
113   const G4double ak2 = 1.0;                    << 
114 }                                              << 
115                                                    81 
                                                   >>  82  
116 //....oooOO0OOooo........oooOO0OOooo........oo     83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117                                                    84 
118 G4MuPairProductionModel::G4MuPairProductionMod <<  85 G4MuPairProductionModel::G4MuPairProductionModel(const G4ParticleDefinition*,
119                                                    86                                                  const G4String& nam)
120   : G4VEmModel(nam),                               87   : G4VEmModel(nam),
121   factorForCross(CLHEP::fine_structure_const*C <<  88   minPairEnergy(4.*electron_mass_c2),
122      CLHEP::classic_electr_radius*CLHEP::class <<  89   highKinEnergy(1000000.*TeV),
123      4./(3.*CLHEP::pi)),                       <<  90   lowKinEnergy(minPairEnergy),
124   sqrte(std::sqrt(G4Exp(1.))),                 <<  91   lowestKinEnergy(1.*GeV),
125   minPairEnergy(4.*CLHEP::electron_mass_c2),   <<  92   factorForCross(4.*fine_structure_const*fine_structure_const
126   lowestKinEnergy(0.85*CLHEP::GeV)             <<  93                    *classic_electr_radius*classic_electr_radius/(3.*pi)),
127 {                                              <<  94   sqrte(sqrt(exp(1.))),
128   nist = G4NistManager::Instance();            <<  95   particleMass(G4MuonPlus::MuonPlus()->GetPDGMass()),
                                                   >>  96   currentZ(0),
                                                   >>  97   particle(G4MuonPlus::MuonPlus()),
                                                   >>  98   nzdat(5),
                                                   >>  99   ntdat(8),
                                                   >> 100   nbiny(1000),
                                                   >> 101   nmaxElements(0),
                                                   >> 102   ymin(-5.),
                                                   >> 103   ymax(0.),
                                                   >> 104   dy((ymax-ymin)/nbiny),
                                                   >> 105   samplingTablesAreFilled(false)
                                                   >> 106 {}
129                                                   107 
130   theElectron = G4Electron::Electron();        << 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131   thePositron = G4Positron::Positron();        << 
132                                                   109 
133   if(nullptr != p) {                           << 110 G4MuPairProductionModel::~G4MuPairProductionModel()
134     SetParticle(p);                            << 111 {}
135     lowestKinEnergy = std::max(lowestKinEnergy << 
136   }                                            << 
137   emin = lowestKinEnergy;                      << 
138   emax = emin*10000.;                          << 
139   SetAngularDistribution(new G4ModifiedMephi() << 
140 }                                              << 
141                                                   112 
142 //....oooOO0OOooo........oooOO0OOooo........oo    113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143                                                   114 
144 G4double G4MuPairProductionModel::MinPrimaryEn << 115 G4double G4MuPairProductionModel::HighEnergyLimit(const G4ParticleDefinition*)
145                                                << 
146                                                << 
147 {                                                 116 {
148   return std::max(lowestKinEnergy, cut);       << 117   return highKinEnergy;
149 }                                                 118 }
150                                                   119 
151 //....oooOO0OOooo........oooOO0OOooo........oo    120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152                                                   121 
153 void G4MuPairProductionModel::Initialise(const << 122 G4double G4MuPairProductionModel::LowEnergyLimit(const G4ParticleDefinition*)
154                                          const << 123 {
155 {                                              << 124   return lowKinEnergy;
156   SetParticle(p);                              << 125 }
157                                                << 
158   if (nullptr == fParticleChange) {            << 
159     fParticleChange = GetParticleChangeForLoss << 
160                                                   126 
161     // define scale of internal table for each << 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162     if (0 == nbine) {                          << 
163       emin = std::max(lowestKinEnergy, LowEner << 
164       emax = std::max(HighEnergyLimit(), emin* << 
165       nbine = std::size_t(nYBinPerDecade*std:: << 
166       if(nbine < 3) { nbine = 3; }             << 
167                                                   128 
168       ymin = G4Log(minPairEnergy/emin);        << 129 G4double G4MuPairProductionModel::MinEnergyCut(const G4ParticleDefinition*,
169       dy = -ymin/G4double(nbiny);              << 130                                                const G4MaterialCutsCouple* )
170     }                                          << 131 {
171     if (p == particle) {                       << 132   return minPairEnergy;
172       G4int pdg = std::abs(p->GetPDGEncoding() << 133 }
173       if (pdg == 2212) {                       << 
174   dataName = "pEEPairProd";                    << 
175       } else if (pdg == 321) {                 << 
176   dataName = "kaonEEPairProd";                 << 
177       } else if (pdg == 211) {                 << 
178   dataName = "pionEEPairProd";                 << 
179       } else if (pdg == 11) {                  << 
180   dataName = "eEEPairProd";                    << 
181       } else if (pdg == 13) {                  << 
182         if (GetName() == "muToMuonPairProd") { << 
183           dataName = "muMuMuPairProd";         << 
184   } else {                                     << 
185     dataName = "muEEPairProd";                 << 
186   }                                            << 
187       }                                        << 
188     }                                          << 
189   }                                            << 
190                                                   134 
191   // for low-energy application this process s << 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192   if(lowestKinEnergy >= HighEnergyLimit()) { r << 
193                                                   136 
194   if (p == particle) {                         << 137 G4bool G4MuPairProductionModel::IsInCharge(const G4ParticleDefinition* p)
195     fElementData =                             << 138 {
196       G4ElementDataRegistry::Instance()->GetEl << 139   return (p == G4MuonMinus::MuonMinus() || p == G4MuonPlus::MuonPlus());
197     if (nullptr == fElementData) {             << 
198       G4AutoLock l(&theMuPairMutex);           << 
199       fElementData =                           << 
200   G4ElementDataRegistry::Instance()->GetElemen << 
201       if (nullptr == fElementData) {           << 
202   fElementData = new G4ElementData(NZDATPAIR); << 
203   fElementData->SetName(dataName);             << 
204       }                                        << 
205       G4bool useDataFile = G4EmParameters::Ins << 
206       if (useDataFile)  { useDataFile = Retrie << 
207       if (!useDataFile) { MakeSamplingTables() << 
208       if (fTableToFile) { StoreTables(); }     << 
209       l.unlock();                              << 
210     }                                          << 
211     if (IsMaster()) {                          << 
212       InitialiseElementSelectors(p, cuts);     << 
213     }                                          << 
214   }                                            << 
215 }                                                 140 }
216                                                   141 
217 //....oooOO0OOooo........oooOO0OOooo........oo << 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218                                                   143 
219 void G4MuPairProductionModel::InitialiseLocal( << 144 void G4MuPairProductionModel::Initialise(const G4ParticleDefinition*,
220                                                << 145                                          const G4DataVector&)
221 {                                              << 146 { 
222   if(p == particle && lowestKinEnergy < HighEn << 147   if (!samplingTablesAreFilled) MakeSamplingTables();
223     SetElementSelectors(masterModel->GetElemen << 
224   }                                            << 
225 }                                                 148 }
226                                                   149 
227 //....oooOO0OOooo........oooOO0OOooo........oo    150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
228                                                   151 
229 G4double G4MuPairProductionModel::ComputeDEDXP << 152 G4double G4MuPairProductionModel::ComputeDEDX(const G4MaterialCutsCouple* couple,
230                                                << 
231                                                   153                                               const G4ParticleDefinition*,
232                                                   154                                                     G4double kineticEnergy,
233                                                   155                                                     G4double cutEnergy)
234 {                                                 156 {
235   G4double dedx = 0.0;                            157   G4double dedx = 0.0;
236   if (cutEnergy <= minPairEnergy || kineticEne << 158   if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy) return dedx;
237     { return dedx; }                           << 
238                                                   159 
                                                   >> 160   const G4Material* material = couple->GetMaterial();
239   const G4ElementVector* theElementVector = ma    161   const G4ElementVector* theElementVector = material->GetElementVector();
240   const G4double* theAtomicNumDensityVector =     162   const G4double* theAtomicNumDensityVector =
241                                    material->G    163                                    material->GetAtomicNumDensityVector();
242                                                   164 
243   //  loop for elements in the material           165   //  loop for elements in the material
244   for (std::size_t i=0; i<material->GetNumberO << 166   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
245      G4double Z = (*theElementVector)[i]->GetZ    167      G4double Z = (*theElementVector)[i]->GetZ();
246      G4double tmax = MaxSecondaryEnergyForElem << 168      SetCurrentElement(Z);
247      G4double loss = ComputMuPairLoss(Z, kinet << 169      G4double tmax = MaxSecondaryEnergy(particle, kineticEnergy);
                                                   >> 170      G4double cut  = std::min(cutEnergy,tmax);
                                                   >> 171      G4double loss = ComputMuPairLoss(Z, kineticEnergy, cut, tmax);
248      dedx += loss*theAtomicNumDensityVector[i]    172      dedx += loss*theAtomicNumDensityVector[i];
249   }                                               173   }
250   dedx = std::max(dedx, 0.0);                  << 174   if (dedx < 0.) dedx = 0.;
251   return dedx;                                    175   return dedx;
252 }                                                 176 }
253                                                   177 
254 //....oooOO0OOooo........oooOO0OOooo........oo    178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255                                                   179 
256 G4double G4MuPairProductionModel::ComputMuPair << 180 G4double G4MuPairProductionModel::ComputMuPairLoss(G4double Z,
257                                                << 181                                                    G4double tkin, G4double cutEnergy,
258                                                << 182                G4double tmax)
259                                                << 
260 {                                                 183 {
                                                   >> 184   SetCurrentElement(Z);
261   G4double loss = 0.0;                            185   G4double loss = 0.0;
262                                                   186 
263   G4double cut = std::min(cutEnergy, tmax);    << 187   G4double cut = cutEnergy;
264   if(cut <= minPairEnergy) { return loss; }    << 188   if(tmax <= cutEnergy || cut <= minPairEnergy) return loss;
265                                                   189 
266   // calculate the rectricted loss                190   // calculate the rectricted loss
267   // numerical integration in log(PairEnergy)     191   // numerical integration in log(PairEnergy)
268   G4double aaa = G4Log(minPairEnergy);         << 192   G4double ak1=6.9;
269   G4double bbb = G4Log(cut);                   << 193   G4double ak2=1.0;
270                                                << 194   G4double aaa = log(minPairEnergy);
271   G4int kkk = std::min(std::max(G4lrint((bbb-a << 195   G4double bbb = log(cut);
272   G4double hhh = (bbb-aaa)/kkk;                << 196   G4int    kkk = (G4int)((bbb-aaa)/ak1+ak2);
                                                   >> 197   if (kkk > 8) kkk = 8;
                                                   >> 198   G4double hhh = (bbb-aaa)/(G4double)kkk;
273   G4double x = aaa;                               199   G4double x = aaa;
274                                                   200 
275   for (G4int l=0 ; l<kkk; ++l) {               << 201   for (G4int l=0 ; l<kkk; l++)
276     for (G4int ll=0; ll<NINTPAIR; ++ll) {      << 202   {
277       G4double ep = G4Exp(x+xgi[ll]*hhh);      << 203 
                                                   >> 204     for (G4int ll=0; ll<8; ll++)
                                                   >> 205     {
                                                   >> 206       G4double ep = exp(x+xgi[ll]*hhh);
278       loss += wgi[ll]*ep*ep*ComputeDMicroscopi    207       loss += wgi[ll]*ep*ep*ComputeDMicroscopicCrossSection(tkin, Z, ep);
279     }                                             208     }
280     x += hhh;                                     209     x += hhh;
281   }                                               210   }
282   loss *= hhh;                                    211   loss *= hhh;
283   loss = std::max(loss, 0.0);                  << 212   if (loss < 0.) loss = 0.;
284   return loss;                                    213   return loss;
285 }                                                 214 }
286                                                   215 
287 //....oooOO0OOooo........oooOO0OOooo........oo    216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288                                                   217 
289 G4double G4MuPairProductionModel::ComputeMicro    218 G4double G4MuPairProductionModel::ComputeMicroscopicCrossSection(
290                                            G4d    219                                            G4double tkin,
291                                            G4d    220                                            G4double Z,
292                                            G4d << 221                                            G4double cut)
                                                   >> 222 
293 {                                                 223 {
294   G4double cross = 0.;                         << 224   G4double cross = 0. ;
295   G4double tmax = MaxSecondaryEnergyForElement << 225 
296   G4double cut  = std::max(cutEnergy, minPairE << 226   SetCurrentElement(Z);
297   if (tmax <= cut) { return cross; }           << 227   G4double tmax = MaxSecondaryEnergy(particle, tkin);
298                                                << 
299   G4double aaa = G4Log(cut);                   << 
300   G4double bbb = G4Log(tmax);                  << 
301   G4int kkk = std::min(std::max(G4lrint((bbb-a << 
302                                                   228 
303   G4double hhh = (bbb-aaa)/(kkk);              << 229   if (tmax <= cut) return cross;
                                                   >> 230 
                                                   >> 231   G4double ak1=6.9 ;
                                                   >> 232   G4double ak2=1.0 ;
                                                   >> 233   G4double aaa = log(cut);
                                                   >> 234   G4double bbb = log(tmax);
                                                   >> 235   G4int kkk = (G4int)((bbb-aaa)/ak1 + ak2);
                                                   >> 236   if(kkk > 8) kkk = 8;
                                                   >> 237   G4double hhh = (bbb-aaa)/float(kkk);
304   G4double x = aaa;                               238   G4double x = aaa;
305                                                   239 
306   for (G4int l=0; l<kkk; ++l) {                << 240   for(G4int l=0; l<kkk; l++)
307     for (G4int i=0; i<NINTPAIR; ++i) {         << 241   {
308       G4double ep = G4Exp(x + xgi[i]*hhh);     << 242     for(G4int i=0; i<8; i++)
                                                   >> 243     {
                                                   >> 244       G4double ep = exp(x + xgi[i]*hhh);
309       cross += ep*wgi[i]*ComputeDMicroscopicCr    245       cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
310     }                                             246     }
311     x += hhh;                                     247     x += hhh;
312   }                                               248   }
313                                                   249 
314   cross *= hhh;                                << 250   cross *=hhh;
315   cross = std::max(cross, 0.0);                << 251   if(cross < 0.0) cross = 0.0;
316   return cross;                                   252   return cross;
317 }                                                 253 }
318                                                   254 
319 //....oooOO0OOooo........oooOO0OOooo........oo << 
320                                                << 
321 G4double G4MuPairProductionModel::ComputeDMicr    255 G4double G4MuPairProductionModel::ComputeDMicroscopicCrossSection(
322                                            G4d    256                                            G4double tkin,
323                                            G4d    257                                            G4double Z,
324                                            G4d    258                                            G4double pairEnergy)
325 // Calculates the  differential (D) microscopi << 259  // Calculates the  differential (D) microscopic cross section
326 // using the cross section formula of R.P. Kok << 260  // using the cross section formula of R.P. Kokoulin (18/01/98)
327 // Code modified by R.P. Kokoulin, V.N. Ivanch << 261  // Code modified by R.P. Kokoulin, V.N. Ivanchenko (27/01/04)
328 {                                              << 262 {
329   static const G4double bbbtf= 183. ;          << 263   G4double bbbtf= 183. ;
330   static const G4double bbbh = 202.4 ;         << 264   G4double bbbh = 202.4 ;
331   static const G4double g1tf = 1.95e-5 ;       << 265   G4double g1tf = 1.95e-5 ;
332   static const G4double g2tf = 5.3e-5 ;        << 266   G4double g2tf = 5.3e-5 ;
333   static const G4double g1h  = 4.4e-5 ;        << 267   G4double g1h  = 4.4e-5 ;
334   static const G4double g2h  = 4.8e-5 ;        << 268   G4double g2h  = 4.8e-5 ;
335                                                << 
336   if (pairEnergy <= minPairEnergy)             << 
337     return 0.0;                                << 
338                                                   269 
339   G4double totalEnergy  = tkin + particleMass;    270   G4double totalEnergy  = tkin + particleMass;
340   G4double residEnergy  = totalEnergy - pairEn    271   G4double residEnergy  = totalEnergy - pairEnergy;
                                                   >> 272   G4double massratio    = particleMass/electron_mass_c2 ;
                                                   >> 273   G4double massratio2   = massratio*massratio ;
                                                   >> 274   G4double cross = 0.;
                                                   >> 275 
                                                   >> 276   SetCurrentElement(Z);
341                                                   277 
342   if (residEnergy <= 0.75*sqrte*z13*particleMa << 278   G4double c3 = 0.75*sqrte*particleMass;
343     return 0.0;                                << 279   if (residEnergy <= c3*z13) return cross;
344                                                   280 
345   G4double a0 = 1.0 / (totalEnergy * residEner << 281   G4double c7 = 4.*electron_mass_c2;
346   G4double alf = 4.0 * electron_mass_c2 / pair << 282   G4double c8 = 6.*particleMass*particleMass;
347   G4double rt = std::sqrt(1.0 - alf);          << 283   G4double alf = c7/pairEnergy;
348   G4double delta = 6.0 * particleMass * partic << 284   G4double a3 = 1. - alf;
349   G4double tmnexp = alf/(1.0 + rt) + delta*rt; << 285   if (a3 <= 0.) return cross;
350                                                << 
351   if(tmnexp >= 1.0) { return 0.0; }            << 
352                                                << 
353   G4double tmn = G4Log(tmnexp);                << 
354                                                << 
355   G4double massratio = particleMass/CLHEP::ele << 
356   G4double massratio2 = massratio*massratio;   << 
357   G4double inv_massratio2 = 1.0 / massratio2;  << 
358                                                   286 
359   // zeta calculation                             287   // zeta calculation
360   G4double bbb,g1,g2;                             288   G4double bbb,g1,g2;
361   if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 =    289   if( Z < 1.5 ) { bbb = bbbh ; g1 = g1h ; g2 = g2h ; }
362   else          { bbb = bbbtf; g1 = g1tf; g2 =    290   else          { bbb = bbbtf; g1 = g1tf; g2 = g2tf; }
363                                                   291 
364   G4double zeta = 0.0;                         << 292   G4double zeta = 0;
365   G4double z1exp = totalEnergy / (particleMass << 293   G4double zeta1 = 0.073 * log(totalEnergy/(particleMass+g1*z23*totalEnergy))-0.26 ;
366                                                << 294   if ( zeta1 > 0.)
367   // 35.221047195922 is the root of zeta1(x) = << 
368   // condition below is the same as zeta1 > 0. << 
369   if (z1exp > 35.221047195922)                 << 
370   {                                               295   {
371     G4double z2exp = totalEnergy / (particleMa << 296     G4double zeta2 = 0.058*log(totalEnergy/(particleMass+g2*z13*totalEnergy))-0.14 ;
372     zeta = (0.073 * G4Log(z1exp) - 0.26) / (0. << 297     zeta  = zeta1/zeta2 ;
373   }                                               298   }
374                                                   299 
375   G4double z2 = Z*(Z+zeta);                       300   G4double z2 = Z*(Z+zeta);
376   G4double screen0 = 2.*electron_mass_c2*sqrte    301   G4double screen0 = 2.*electron_mass_c2*sqrte*bbb/(z13*pairEnergy);
377   G4double beta = 0.5*pairEnergy*pairEnergy*a0 << 302   G4double a0 = totalEnergy*residEnergy;
378   G4double xi0 = 0.5*massratio2*beta;          << 303   G4double a1 = pairEnergy*pairEnergy/a0;
                                                   >> 304   G4double bet = 0.5*a1;
                                                   >> 305   G4double xi0 = 0.25*massratio2*a1;
                                                   >> 306   G4double del = c8/a0;
                                                   >> 307 
                                                   >> 308   G4double tmnexp = (alf+2.*del*a3)/(1.+(1.-del)*sqrt(a3));
                                                   >> 309   if(tmnexp <= 0.) return cross;
                                                   >> 310   G4double tmn = log(tmnexp);
                                                   >> 311   G4double sum = 0.;
379                                                   312 
380   // Gaussian integration in ln(1-ro) ( with 8    313   // Gaussian integration in ln(1-ro) ( with 8 points)
381   G4double rho[NINTPAIR];                      << 314   for (G4int i=0; i<8; i++)
382   G4double rho2[NINTPAIR];                     << 
383   G4double xi[NINTPAIR];                       << 
384   G4double xi1[NINTPAIR];                      << 
385   G4double xii[NINTPAIR];                      << 
386                                                << 
387   for (G4int i = 0; i < NINTPAIR; ++i)         << 
388   {                                               315   {
389     rho[i] = G4Exp(tmn*xgi[i]) - 1.0; // rho = << 316     G4double a4 = exp(tmn*xgi[i]);     // a4 = (1.-asymmetry)
390     rho2[i] = rho[i] * rho[i];                 << 317     G4double a5 = a4*(2.-a4) ;
391     xi[i] = xi0*(1.0-rho2[i]);                 << 318     G4double a6 = 1.-a5 ;
392     xi1[i] = 1.0 + xi[i];                      << 319     G4double a7 = 1.+a6 ;
393     xii[i] = 1.0 / xi[i];                      << 320     G4double a9 = 3.+a6 ;
                                                   >> 321     G4double xi = xi0*a5 ;
                                                   >> 322     G4double xii = 1./xi ;
                                                   >> 323     G4double xi1 = 1.+xi ;
                                                   >> 324     G4double screen = screen0*xi1/a5 ;
                                                   >> 325     G4double yeu = 5.-a6+4.*bet*a7 ;
                                                   >> 326     G4double yed = 2.*(1.+3.*bet)*log(3.+xii)-a6-a1*(2.-a6) ;
                                                   >> 327     G4double ye1 = 1.+yeu/yed ;
                                                   >> 328     G4double ale=log(bbb/z13*sqrt(xi1*ye1)/(1.+screen*ye1)) ;
                                                   >> 329     G4double cre = 0.5*log(1.+2.25*z23*xi1*ye1/massratio2) ;
                                                   >> 330     G4double be;
                                                   >> 331 
                                                   >> 332     if (xi <= 1.e3) be = ((2.+a6)*(1.+bet)+xi*a9)*log(1.+xii)+(a5-bet)/xi1-a9;
                                                   >> 333     else            be = (3.-a6+a1*a7)/(2.*xi);
                                                   >> 334 
                                                   >> 335     G4double fe = (ale-cre)*be;
                                                   >> 336     if ( fe < 0.) fe = 0. ;
                                                   >> 337 
                                                   >> 338     G4double ymu = 4.+a6 +3.*bet*a7 ;
                                                   >> 339     G4double ymd = a7*(1.5+a1)*log(3.+xi)+1.-1.5*a6 ;
                                                   >> 340     G4double ym1 = 1.+ymu/ymd ;
                                                   >> 341     G4double alm_crm = log(bbb*massratio/(1.5*z23*(1.+screen*ym1)));
                                                   >> 342     G4double a10,bm;
                                                   >> 343     if ( xi >= 1.e-3)
                                                   >> 344     {
                                                   >> 345       a10 = (1.+a1)*a5 ;
                                                   >> 346       bm  = (a7*(1.+1.5*bet)-a10*xii)*log(xi1)+xi*(a5-bet)/xi1+a10;
                                                   >> 347     } else {
                                                   >> 348       bm = (5.-a6+bet*a9)*(xi/2.);
                                                   >> 349     }
                                                   >> 350 
                                                   >> 351     G4double fm = alm_crm*bm;
                                                   >> 352     if ( fm < 0.) fm = 0. ;
                                                   >> 353 
                                                   >> 354     sum += wgi[i]*a4*(fe+fm/massratio2);
394   }                                               355   }
395                                                   356 
396   G4double ye1[NINTPAIR];                      << 357   cross = -tmn*sum*factorForCross*z2*residEnergy/(totalEnergy*pairEnergy);
397   G4double ym1[NINTPAIR];                      << 
398                                                   358 
399   G4double b40 = 4.0 * beta;                   << 359   return cross;
400   G4double b62 = 6.0 * beta + 2.0;             << 360 }
401                                                   361 
402   for (G4int i = 0; i < NINTPAIR; ++i)         << 362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
403   {                                            << 
404     G4double yeu = (b40 + 5.0) + (b40 - 1.0) * << 
405     G4double yed = b62*G4Log(3.0 + xii[i]) + ( << 
406                                                   363 
407     G4double ymu = b62 * (1.0 + rho2[i]) + 6.0 << 364 G4double G4MuPairProductionModel::CrossSection(const G4MaterialCutsCouple* couple,
408     G4double ymd = (b40 + 3.0)*(1.0 + rho2[i]) << 365                                                const G4ParticleDefinition*,
409       + 2.0 - 3.0 * rho2[i];                   << 366                                                      G4double kineticEnergy,
                                                   >> 367                                                      G4double cutEnergy,
                                                   >> 368                  G4double maxEnergy)
                                                   >> 369 {
                                                   >> 370   G4double cross = 0.0;
                                                   >> 371   if (kineticEnergy <= lowestKinEnergy) return cross;
410                                                   372 
411     ye1[i] = 1.0 + yeu / yed;                  << 373   maxEnergy += particleMass;
412     ym1[i] = 1.0 + ymu / ymd;                  << 
413   }                                            << 
414                                                   374 
415   G4double be[NINTPAIR];                       << 375   const G4Material* material = couple->GetMaterial();
416   G4double bm[NINTPAIR];                       << 376   const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 377   const G4double* theAtomNumDensityVector=material->GetAtomicNumDensityVector();
417                                                   378 
418   for(G4int i = 0; i < NINTPAIR; ++i) {        << 379   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
419     if(xi[i] <= 1000.0) {                      << 380     G4double Z = (*theElementVector)[i]->GetZ();
420       be[i] = ((2.0 + rho2[i])*(1.0 + beta) +  << 381     SetCurrentElement(Z);
421          xi[i]*(3.0 + rho2[i]))*G4Log(1.0 + xi << 382     G4double tmax = std::min(maxEnergy,MaxSecondaryEnergy(particle, kineticEnergy));
422   (1.0 - rho2[i] - beta)/xi1[i] - (3.0 + rho2[ << 383     G4double cut  = std::max(minPairEnergy,cutEnergy);
423     } else {                                   << 384     if(cut < tmax) {
424       be[i] = 0.5*(3.0 - rho2[i] + 2.0*beta*(1 << 385       G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut)
425     }                                          << 386                   - ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
426                                                   387 
427     if(xi[i] >= 0.001) {                       << 388       cross += theAtomNumDensityVector[i] * cr;
428       G4double a10 = (1.0 + 2.0 * beta) * (1.0 << 
429       bm[i] = ((1.0 + rho2[i])*(1.0 + 1.5 * be << 
430                 xi[i] * (1.0 - rho2[i] - beta) << 
431     } else {                                   << 
432       bm[i] = 0.5*(5.0 - rho2[i] + beta * (3.0 << 
433     }                                             389     }
434   }                                               390   }
                                                   >> 391   return cross;
                                                   >> 392 }
435                                                   393 
436   G4double sum = 0.0;                          << 394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 395 
                                                   >> 396 void G4MuPairProductionModel::MakeSamplingTables()
                                                   >> 397 {
                                                   >> 398   for (G4int iz=0; iz<nzdat; iz++)
                                                   >> 399   {
                                                   >> 400     G4double Z = zdat[iz];
                                                   >> 401     SetCurrentElement(Z);
437                                                   402 
438   for (G4int i = 0; i < NINTPAIR; ++i) {       << 403     for (G4int it=0; it<ntdat; it++)
439     G4double screen = screen0*xi1[i]/(1.0 - rh << 404     {
440     G4double ale = G4Log(bbb/z13*std::sqrt(xi1 << 405       G4double kineticEnergy = tdat[it];
441     G4double cre = 0.5*G4Log(1. + 2.25*z23*xi1 << 406       G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
                                                   >> 407 
                                                   >> 408       G4double CrossSection = 0.0 ;
                                                   >> 409 
                                                   >> 410       G4double y = ymin - 0.5*dy ;
                                                   >> 411       G4double yy = ymin - dy ;
                                                   >> 412       G4double x = exp(y);
                                                   >> 413       G4double fac = exp(dy);
                                                   >> 414       G4double dx = exp(yy)*(fac - 1.0);
                                                   >> 415 
                                                   >> 416       G4double c = log(maxPairEnergy/minPairEnergy);
                                                   >> 417 
                                                   >> 418       for (G4int i=0 ; i<nbiny; i++)
                                                   >> 419       {
                                                   >> 420         y += dy ;
                                                   >> 421         if(c > 0.0) {
                                                   >> 422     x *= fac;
                                                   >> 423           dx*= fac;
                                                   >> 424           G4double ep = minPairEnergy*exp(c*x) ;
                                                   >> 425           CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
                                                   >> 426                                       kineticEnergy, Z, ep);
                                                   >> 427         }
                                                   >> 428         ya[i] = y;
                                                   >> 429         proba[iz][it][i] = CrossSection;
                                                   >> 430       }
442                                                   431 
443     G4double fe = (ale-cre)*be[i];             << 432       ya[nbiny]=ymax;
444     fe = std::max(fe, 0.0);                    << 
445                                                   433 
446     G4double alm_crm = G4Log(bbb*massratio/(1. << 434       proba[iz][it][nbiny] = CrossSection;
447     G4double fm = std::max(alm_crm*bm[i], 0.0) << 
448                                                   435 
449     sum += wgi[i]*(1.0 + rho[i])*(fe + fm);    << 436     }
450   }                                               437   }
451                                                << 438   samplingTablesAreFilled = true;
452   return -tmn*sum*factorForCross*z2*residEnerg << 
453 }                                                 439 }
454                                                   440 
455 //....oooOO0OOooo........oooOO0OOooo........oo    441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
456                                                   442 
457 G4double G4MuPairProductionModel::ComputeCross << 443 G4DynamicParticle* G4MuPairProductionModel::SampleSecondary(
458                                            con << 444                              const G4MaterialCutsCouple*,
459                                                << 445                              const G4DynamicParticle*,
460                                                << 446                                    G4double,
461                                                << 447                                    G4double)
462                                                << 
463 {                                                 448 {
464   G4double cross = 0.0;                        << 449   return 0;
465   if (kineticEnergy <= lowestKinEnergy) { retu << 
466                                                << 
467   G4double maxPairEnergy = MaxSecondaryEnergyF << 
468   G4double tmax = std::min(maxEnergy, maxPairE << 
469   G4double cut  = std::max(cutEnergy, minPairE << 
470   if (cut >= tmax) { return cross; }           << 
471                                                << 
472   cross = ComputeMicroscopicCrossSection(kinet << 
473   if(tmax < kineticEnergy) {                   << 
474     cross -= ComputeMicroscopicCrossSection(ki << 
475   }                                            << 
476   return cross;                                << 
477 }                                                 450 }
478                                                   451 
479 //....oooOO0OOooo........oooOO0OOooo........oo    452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
480                                                   453 
481 void G4MuPairProductionModel::MakeSamplingTabl << 454 std::vector<G4DynamicParticle*>* G4MuPairProductionModel::SampleSecondaries(
                                                   >> 455                              const G4MaterialCutsCouple* couple,
                                                   >> 456                              const G4DynamicParticle* aDynamicParticle,
                                                   >> 457                                    G4double cut,
                                                   >> 458                                    G4double tmax)
482 {                                                 459 {
483   G4double factore = G4Exp(G4Log(emax/emin)/G4 << 460   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
                                                   >> 461   G4double totalEnergy   = kineticEnergy + particleMass ;
                                                   >> 462   G4ParticleMomentum ParticleDirection = aDynamicParticle->GetMomentumDirection();
484                                                   463 
485   for (G4int iz=0; iz<NZDATPAIR; ++iz) {       << 464   G4int it;
                                                   >> 465   for(it=1; it<ntdat; it++) {if(kineticEnergy <= tdat[it]) break;}
                                                   >> 466   if(it == ntdat) it--;
                                                   >> 467   G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
486                                                   468 
487     G4double Z = ZDATPAIR[iz];                 << 469   // select randomly one element constituing the material
488     G4Physics2DVector* pv = new G4Physics2DVec << 470   G4int iymin = 0;
489     G4double kinEnergy = emin;                 << 471   G4int iymax = nbiny-1;
490                                                << 472   const G4Element* anElement = SelectRandomAtom(kineticEnergy, dt, it, couple);
491     for (std::size_t it=0; it<=nbine; ++it) {  << 473   SetCurrentElement(anElement->GetZ());
492                                                << 474 
493       pv->PutY(it, G4Log(kinEnergy/CLHEP::MeV) << 475   G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
494       G4double maxPairEnergy = MaxSecondaryEne << 476   G4double maxEnergy     = std::min(tmax, maxPairEnergy);
495       /*                                       << 477   G4double minEnergy     = std::min(maxEnergy, cut);
496       G4cout << "it= " << it << " E= " << kinE << 
497              << "  " << particle->GetParticleN << 
498              << " maxE= " << maxPairEnergy <<  << 
499              << " ymin= " << ymin << G4endl;   << 
500       */                                       << 
501       G4double coef = G4Log(minPairEnergy/kinE << 
502       G4double ymax = G4Log(maxPairEnergy/kinE << 
503       G4double fac  = (ymax - ymin)/dy;        << 
504       std::size_t imax   = (std::size_t)fac;   << 
505       fac -= (G4double)imax;                   << 
506                                                << 
507       G4double xSec = 0.0;                     << 
508       G4double x = ymin;                       << 
509       /*                                       << 
510       G4cout << "Z= " << currentZ << " z13= "  << 
511              << " mE= " << maxPairEnergy << "  << 
512              << " dy= " << dy << "  c= " << co << 
513       */                                       << 
514       // start from zero                       << 
515       pv->PutValue(0, it, 0.0);                << 
516       if(0 == it) { pv->PutX(nbiny, 0.0); }    << 
517                                                << 
518       for (std::size_t i=0; i<nbiny; ++i) {    << 
519                                                << 
520         if(0 == it) { pv->PutX(i, x); }        << 
521                                                << 
522         if(i < imax) {                         << 
523           G4double ep = kinEnergy*G4Exp(coef*( << 
524                                                << 
525           // not multiplied by interval, becau << 
526           // will be used only for sampling    << 
527           //G4cout << "i= " << i << " x= " <<  << 
528           //         << " Egamma= " << ep << G << 
529           xSec += ep*ComputeDMicroscopicCrossS << 
530                                                << 
531           // last bin before the kinematic lim << 
532         } else if(i == imax) {                 << 
533           G4double ep = kinEnergy*G4Exp(coef*( << 
534           xSec += ep*fac*ComputeDMicroscopicCr << 
535         }                                      << 
536         pv->PutValue(i + 1, it, xSec);         << 
537         x += dy;                               << 
538       }                                        << 
539       kinEnergy *= factore;                    << 
540                                                   478 
541       // to avoid precision lost               << 479   if( minEnergy > minPairEnergy)
542       if(it+1 == nbine) { kinEnergy = emax; }  << 480   {
543     }                                          << 481     G4double xc = log(minEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
544     fElementData->InitialiseForElement(iz, pv) << 482     iymin = (G4int)((log(xc) - ymin)/dy);
                                                   >> 483     if(iymin >= nbiny) iymin = nbiny-1;
                                                   >> 484     xc = log(maxEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
                                                   >> 485     iymax = (G4int)((log(xc) - ymin)/dy) + 1;
                                                   >> 486     if(iymax >= nbiny) iymax = nbiny-1;
545   }                                               487   }
546 }                                              << 
547                                                   488 
548 //....oooOO0OOooo........oooOO0OOooo........oo << 489   // sample e-e+ energy, pair energy first
                                                   >> 490   G4int iz, iy;
549                                                   491 
550 void G4MuPairProductionModel::SampleSecondarie << 492   for(iz=1; iz<nzdat; iz++) {if(currentZ <= zdat[iz]) break;}
551                               std::vector<G4Dy << 493   if(iz == nzdat) iz--;
552                               const G4Material << 
553                               const G4DynamicP << 
554                               G4double tmin,   << 
555                               G4double tmax)   << 
556 {                                              << 
557   G4double kinEnergy = aDynamicParticle->GetKi << 
558   //G4cout << "------- G4MuPairProductionModel << 
559   //         << kinEnergy << "  "              << 
560   //         << aDynamicParticle->GetDefinitio << 
561   G4double totalEnergy   = kinEnergy + particl << 
562   G4double totalMomentum =                     << 
563     std::sqrt(kinEnergy*(kinEnergy + 2.0*parti << 
564                                                   494 
565   G4ThreeVector partDirection = aDynamicPartic << 495   G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
566                                                   496 
567   // select randomly one element constituing t << 497   G4double pmin = InterpolatedIntegralCrossSection(dt, dz, iz, it, iymin, currentZ);
568   const G4Element* anElement = SelectRandomAto << 498   G4double pmax = InterpolatedIntegralCrossSection(dt, dz, iz, it, iymax, currentZ);
569                                                   499 
570   // define interval of energy transfer        << 500   G4double p = pmin+G4UniformRand()*(pmax - pmin);
571   G4double maxPairEnergy = MaxSecondaryEnergyF << 
572                                                << 
573   G4double maxEnergy = std::min(tmax, maxPairE << 
574   G4double minEnergy = std::max(tmin, minPairE << 
575                                                << 
576   if (minEnergy >= maxEnergy) { return; }      << 
577   //G4cout << "emin= " << minEnergy << " emax= << 
578   // << " minPair= " << minPairEnergy << " max << 
579   //    << " ymin= " << ymin << " dy= " << dy  << 
580                                                << 
581   G4double coeff = G4Log(minPairEnergy/kinEner << 
582                                                << 
583   // compute limits                            << 
584   G4double yymin = G4Log(minEnergy/kinEnergy)/ << 
585   G4double yymax = G4Log(maxEnergy/kinEnergy)/ << 
586                                                << 
587   //G4cout << "yymin= " << yymin << "  yymax=  << 
588                                                   501 
589   // units should not be used, bacause table w << 502   // interpolate sampling vector;
590   G4double logTkin = G4Log(kinEnergy/CLHEP::Me << 503   G4double p1 = pmin;
                                                   >> 504   G4double p2 = pmin;
                                                   >> 505   for(iy=iymin+1; iy<=iymax; iy++) {
                                                   >> 506     p1 = p2;
                                                   >> 507     p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
                                                   >> 508     if(p <= p2) break;
                                                   >> 509   }
                                                   >> 510   G4double y = ya[iy-1] + dy*(p - p1)/(p2 - p1);
591                                                   511 
592   // sample e-e+ energy, pair energy first     << 512   G4double PairEnergy = minPairEnergy*exp(exp(y)*log(maxPairEnergy/minPairEnergy));
                                                   >> 513   if(PairEnergy < minEnergy) PairEnergy = minEnergy;
                                                   >> 514   if(PairEnergy > maxEnergy) PairEnergy = maxEnergy;
593                                                   515 
594   // select sample table via Z                 << 516   // sample r=(E+-E-)/PairEnergy  ( uniformly .....)
595   G4int iz1(0), iz2(0);                        << 517   G4double rmax = 
596   for (G4int iz=0; iz<NZDATPAIR; ++iz) {       << 518     (1.-6.*particleMass*particleMass/(totalEnergy*(totalEnergy-PairEnergy)))
597     if(currentZ == ZDATPAIR[iz]) {             << 519                                        *sqrt(1.-minPairEnergy/PairEnergy);
598       iz1 = iz2 = iz;                          << 
599       break;                                   << 
600     } else if(currentZ < ZDATPAIR[iz]) {       << 
601       iz2 = iz;                                << 
602       if(iz > 0) { iz1 = iz-1; }               << 
603       else { iz1 = iz2; }                      << 
604       break;                                   << 
605     }                                          << 
606   }                                            << 
607   if (0 == iz1) { iz1 = iz2 = NZDATPAIR-1; }   << 
608                                                << 
609   G4double pairEnergy = 0.0;                   << 
610   G4int count = 0;                             << 
611   //G4cout << "start loop Z1= " << iz1 << " Z2 << 
612   do {                                         << 
613     ++count;                                   << 
614     // sampling using only one random number   << 
615     G4double rand = G4UniformRand();           << 
616                                                << 
617     G4double x = FindScaledEnergy(iz1, rand, l << 
618     if(iz1 != iz2) {                           << 
619       G4double x2 = FindScaledEnergy(iz2, rand << 
620       G4double lz1= nist->GetLOGZ(ZDATPAIR[iz1 << 
621       G4double lz2= nist->GetLOGZ(ZDATPAIR[iz2 << 
622       //G4cout << count << ".  x= " << x << "  << 
623       //             << " Z1= " << iz1 << " Z2 << 
624       x += (x2 - x)*(lnZ - lz1)/(lz2 - lz1);   << 
625     }                                          << 
626     //G4cout << "x= " << x << "  coeff= " << c << 
627     pairEnergy = kinEnergy*G4Exp(x*coeff);     << 
628                                                << 
629     // Loop checking, 03-Aug-2015, Vladimir Iv << 
630   } while((pairEnergy < minEnergy || pairEnerg << 
631                                                << 
632   //G4cout << "## pairEnergy(GeV)= " << pairEn << 
633   //         << " Etot(GeV)= " << totalEnergy/ << 
634                                                << 
635   // sample r=(E+-E-)/pairEnergy  ( uniformly  << 
636   G4double rmax =                              << 
637     (1.-6.*particleMass*particleMass/(totalEne << 
638     *std::sqrt(1.-minPairEnergy/pairEnergy);   << 
639   G4double r = rmax * (-1.+2.*G4UniformRand())    520   G4double r = rmax * (-1.+2.*G4UniformRand()) ;
640                                                   521 
641   // compute energies from pairEnergy,r        << 522   // compute energies from PairEnergy,r
642   G4double eEnergy = (1.-r)*pairEnergy*0.5;    << 523   G4double ElectronEnergy = (1.-r)*PairEnergy*0.5;
643   G4double pEnergy = pairEnergy - eEnergy;     << 524   G4double PositronEnergy = PairEnergy - ElectronEnergy;
644                                                << 525 
645   // Sample angles                             << 526   //  angles of the emitted particles ( Z - axis along the parent particle)
646   G4ThreeVector eDirection, pDirection;        << 527   //      (mean theta for the moment)
647   //                                           << 528   G4double Teta = electron_mass_c2/totalEnergy ;
648   GetAngularDistribution()->SamplePairDirectio << 529 
649                                                << 530   G4double Phi  = twopi * G4UniformRand() ;
650                                                << 531   G4double dirx = sin(Teta)*cos(Phi);
651   // create G4DynamicParticle object for e+e-  << 532   G4double diry = sin(Teta)*sin(Phi);
652   eEnergy = std::max(eEnergy - CLHEP::electron << 533   G4double dirz = cos(Teta) ;
653   pEnergy = std::max(pEnergy - CLHEP::electron << 534 
654   auto aParticle1 = new G4DynamicParticle(theE << 535   //G4double finalPx,finalPy,finalPz ;
655   auto aParticle2 = new G4DynamicParticle(theP << 536   G4double ElectKineEnergy = ElectronEnergy - electron_mass_c2 ;
656   // Fill output vector                        << 537 
                                                   >> 538   G4ThreeVector ElectDirection ( dirx, diry, dirz );
                                                   >> 539   ElectDirection.rotateUz(ParticleDirection);
                                                   >> 540 
                                                   >> 541   // create G4DynamicParticle object for the particle1
                                                   >> 542   G4DynamicParticle* aParticle1= new G4DynamicParticle(G4Electron::Electron(),
                                                   >> 543                                                        ElectDirection,
                                                   >> 544                    ElectKineEnergy);
                                                   >> 545 
                                                   >> 546   G4double PositKineEnergy = PositronEnergy - electron_mass_c2 ;
                                                   >> 547 
                                                   >> 548   G4ThreeVector PositDirection ( -dirx, -diry, dirz );
                                                   >> 549   PositDirection.rotateUz(ParticleDirection);
                                                   >> 550 
                                                   >> 551   // create G4DynamicParticle object for the particle2
                                                   >> 552   G4DynamicParticle* aParticle2= new G4DynamicParticle(G4Positron::Positron(),
                                                   >> 553                                                        PositDirection,
                                                   >> 554                    PositKineEnergy);
                                                   >> 555 
                                                   >> 556 
                                                   >> 557   std::vector<G4DynamicParticle*>* vdp = new std::vector<G4DynamicParticle*>;
657   vdp->push_back(aParticle1);                     558   vdp->push_back(aParticle1);
658   vdp->push_back(aParticle2);                     559   vdp->push_back(aParticle2);
659                                                   560 
660   // primary change                            << 561   return vdp;
661   kinEnergy -= pairEnergy;                     << 
662   partDirection *= totalMomentum;              << 
663   partDirection -= (aParticle1->GetMomentum()  << 
664   partDirection = partDirection.unit();        << 
665                                                << 
666   // if energy transfer is higher than thresho << 
667   // then stop tracking the primary particle a << 
668   if (pairEnergy > SecondaryThreshold()) {     << 
669     fParticleChange->ProposeTrackStatus(fStopA << 
670     fParticleChange->SetProposedKineticEnergy( << 
671     auto newdp = new G4DynamicParticle(particl << 
672     vdp->push_back(newdp);                     << 
673   } else { // continue tracking the primary e- << 
674     fParticleChange->SetProposedMomentumDirect << 
675     fParticleChange->SetProposedKineticEnergy( << 
676   }                                            << 
677   //G4cout << "-- G4MuPairProductionModel::Sam << 
678 }                                                 562 }
679                                                   563 
680 //....oooOO0OOooo........oooOO0OOooo........oo    564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
681                                                   565 
682 G4double                                       << 566 const G4Element* G4MuPairProductionModel::SelectRandomAtom(
683 G4MuPairProductionModel::FindScaledEnergy(G4in << 567                  G4double kinEnergy, G4double dt, G4int it,
684             G4double logTkin,                  << 568            const G4MaterialCutsCouple* couple)
685             G4double yymin, G4double yymax)    << 
686 {                                                 569 {
687   G4double res = yymin;                        << 570   // select randomly 1 element within the material
688   G4Physics2DVector* pv = fElementData->GetEle << 571 
689   if (nullptr != pv) {                         << 572   const G4Material* material = couple->GetMaterial();
690     G4double pmin = pv->Value(yymin, logTkin); << 573   size_t nElements = material->GetNumberOfElements();
691     G4double pmax = pv->Value(yymax, logTkin); << 574   const G4ElementVector* theElementVector = material->GetElementVector();
692     G4double p0   = pv->Value(0.0, logTkin);   << 575   if (nElements == 1) return (*theElementVector)[0];
693     if(p0 <= 0.0) { DataCorrupted(ZDATPAIR[iz] << 576 
694     else { res = pv->FindLinearX((pmin + rand* << 577   if(nElements > nmaxElements) {
695   } else {                                     << 578     nmaxElements = nElements;
696     DataCorrupted(ZDATPAIR[iz], logTkin);      << 579     partialSum.resize(nmaxElements);
697   }                                               580   }
698   return res;                                  << 
699 }                                              << 
700                                                   581 
701 //....oooOO0OOooo........oooOO0OOooo........oo << 582   const G4double* theAtomNumDensityVector=material->GetAtomicNumDensityVector();
702                                                   583 
703 void G4MuPairProductionModel::DataCorrupted(G4 << 584   G4double sum = 0.0;
704 {                                              << 
705   G4ExceptionDescription ed;                   << 
706   ed << "G4ElementData is not properly initial << 
707      << " Ekin(MeV)= " << G4Exp(logTkin)       << 
708      << " IsMasterThread= " << IsMaster()      << 
709      << " Model " << GetName();                << 
710   G4Exception("G4MuPairProductionModel::()", " << 
711 }                                              << 
712                                                   585 
713 //....oooOO0OOooo........oooOO0OOooo........oo << 586   size_t i;
                                                   >> 587   for (i=0; i<nElements; i++) {
                                                   >> 588     G4double Z = ((*theElementVector)[i])->GetZ();
                                                   >> 589     SetCurrentElement(Z);
                                                   >> 590     G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
                                                   >> 591 
                                                   >> 592     G4int iz;
                                                   >> 593     for(iz=1; iz<nzdat; iz++) {if(Z <= zdat[iz]) break;}
                                                   >> 594     if(iz == nzdat) iz--;
                                                   >> 595     G4double dz = log(Z/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
                                                   >> 596 
                                                   >> 597     G4double xc = log(kinEnergy/minPairEnergy)/log(maxPairEnergy/minPairEnergy);
                                                   >> 598     G4int iy = (G4int)((log(xc) - ymin)/dy);
                                                   >> 599     if(iy >= nbiny) iy = nbiny-1;
                                                   >> 600 
                                                   >> 601     G4double sigtot = InterpolatedIntegralCrossSection(dt, dz, iz, it, nbiny, Z);
                                                   >> 602     G4double sigcut = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, Z);
                                                   >> 603     sum += (sigtot - sigcut)*theAtomNumDensityVector[i];
                                                   >> 604     partialSum[i] = sum;
                                                   >> 605   }
                                                   >> 606 
                                                   >> 607   G4double rval = G4UniformRand()*sum;
                                                   >> 608   for (i=0; i<nElements; i++) {if(rval<=partialSum[i]) break;}
                                                   >> 609   return (*theElementVector)[i];
714                                                   610 
715 void G4MuPairProductionModel::StoreTables() co << 
716 {                                              << 
717   for (G4int iz=0; iz<NZDATPAIR; ++iz) {       << 
718     G4int Z = ZDATPAIR[iz];                    << 
719     G4Physics2DVector* pv = fElementData->GetE << 
720     if(nullptr == pv) {                        << 
721       DataCorrupted(Z, 1.0);                   << 
722       return;                                  << 
723     }                                          << 
724     std::ostringstream ss;                     << 
725     ss << "mupair/" << particle->GetParticleNa << 
726     std::ofstream outfile(ss.str());           << 
727     pv->Store(outfile);                        << 
728   }                                            << 
729 }                                                 611 }
730                                                   612 
731 //....oooOO0OOooo........oooOO0OOooo........oo    613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
732                                                   614 
733 G4bool G4MuPairProductionModel::RetrieveTables << 
734 {                                              << 
735   for (G4int iz=0; iz<NZDATPAIR; ++iz) {       << 
736     G4double Z = ZDATPAIR[iz];                 << 
737     G4Physics2DVector* pv = new G4Physics2DVec << 
738     std::ostringstream ss;                     << 
739     ss << G4EmParameters::Instance()->GetDirLE << 
740        << particle->GetParticleName() << Z <<  << 
741     std::ifstream infile(ss.str(), std::ios::i << 
742     if(!pv->Retrieve(infile)) {                << 
743       delete pv;                               << 
744       return false;                            << 
745     }                                          << 
746     fElementData->InitialiseForElement(iz, pv) << 
747   }                                            << 
748   return true;                                 << 
749 }                                              << 
750                                                   615 
751 //....oooOO0OOooo........oooOO0OOooo........oo << 
752                                                   616