Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4PAIPhotModel.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/standard/src/G4PAIPhotModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4PAIPhotModel.cc (Version 10.3.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4PAIPhotModel.cc 73607 2013-09-02 10:04:03Z gcosmo $
 26 //                                                 27 //
 27 // -------------------------------------------     28 // -------------------------------------------------------------------
 28 //                                                 29 //
 29 // GEANT4 Class                                    30 // GEANT4 Class
 30 // File name:     G4PAIPhotModel.cc                31 // File name:     G4PAIPhotModel.cc
 31 //                                                 32 //
 32 // Author: Vladimir.Grichine@cern.ch on base o     33 // Author: Vladimir.Grichine@cern.ch on base of G4PAIModel MT interface
 33 //                                                 34 //
 34 // Creation date: 07.10.2013                       35 // Creation date: 07.10.2013
 35 //                                                 36 //
 36 // Modifications:                                  37 // Modifications:
 37 //                                                 38 //
 38 //                                                 39 //
 39                                                    40 
 40 #include "G4PAIPhotModel.hh"                       41 #include "G4PAIPhotModel.hh"
 41                                                    42 
 42 #include "G4SystemOfUnits.hh"                      43 #include "G4SystemOfUnits.hh"
 43 #include "G4PhysicalConstants.hh"                  44 #include "G4PhysicalConstants.hh"
 44 #include "G4Region.hh"                             45 #include "G4Region.hh"
                                                   >>  46 #include "G4PhysicsLogVector.hh"
                                                   >>  47 #include "G4PhysicsFreeVector.hh"
                                                   >>  48 #include "G4PhysicsTable.hh"
 45 #include "G4ProductionCutsTable.hh"                49 #include "G4ProductionCutsTable.hh"
 46 #include "G4MaterialCutsCouple.hh"                 50 #include "G4MaterialCutsCouple.hh"
 47 #include "G4MaterialTable.hh"                      51 #include "G4MaterialTable.hh"
                                                   >>  52 #include "G4SandiaTable.hh"
                                                   >>  53 #include "G4OrderedTable.hh"
 48 #include "G4RegionStore.hh"                        54 #include "G4RegionStore.hh"
 49                                                    55 
 50 #include "Randomize.hh"                            56 #include "Randomize.hh"
 51 #include "G4Electron.hh"                           57 #include "G4Electron.hh"
 52 #include "G4Positron.hh"                           58 #include "G4Positron.hh"
 53 #include "G4Gamma.hh"                              59 #include "G4Gamma.hh"
 54 #include "G4Poisson.hh"                            60 #include "G4Poisson.hh"
 55 #include "G4Step.hh"                               61 #include "G4Step.hh"
 56 #include "G4Material.hh"                           62 #include "G4Material.hh"
 57 #include "G4DynamicParticle.hh"                    63 #include "G4DynamicParticle.hh"
 58 #include "G4ParticleDefinition.hh"                 64 #include "G4ParticleDefinition.hh"
 59 #include "G4ParticleChangeForLoss.hh"              65 #include "G4ParticleChangeForLoss.hh"
 60 #include "G4PAIPhotData.hh"                        66 #include "G4PAIPhotData.hh"
 61 #include "G4DeltaAngle.hh"                         67 #include "G4DeltaAngle.hh"
 62                                                    68 
 63 //////////////////////////////////////////////     69 ////////////////////////////////////////////////////////////////////////
 64                                                    70 
 65 using namespace std;                               71 using namespace std;
 66                                                    72 
 67 G4PAIPhotModel::G4PAIPhotModel(const G4Particl     73 G4PAIPhotModel::G4PAIPhotModel(const G4ParticleDefinition* p, const G4String& nam)
 68   : G4VEmModel(nam),G4VEmFluctuationModel(nam)     74   : G4VEmModel(nam),G4VEmFluctuationModel(nam),
 69     fVerbose(0),                                   75     fVerbose(0),
 70     fModelData(nullptr),                           76     fModelData(nullptr),
 71     fParticle(nullptr)                             77     fParticle(nullptr)
 72 {                                                  78 {  
 73   fElectron = G4Electron::Electron();              79   fElectron = G4Electron::Electron();
 74   fPositron = G4Positron::Positron();              80   fPositron = G4Positron::Positron();
 75                                                    81 
 76   fParticleChange = nullptr;                       82   fParticleChange = nullptr;
 77                                                    83 
 78   if(p) { SetParticle(p); }                        84   if(p) { SetParticle(p); }
 79   else  { SetParticle(fElectron); }                85   else  { SetParticle(fElectron); }
 80                                                    86 
 81   // default generator                             87   // default generator
 82   SetAngularDistribution(new G4DeltaAngle());      88   SetAngularDistribution(new G4DeltaAngle());
 83   fLowestTcut = 12.5*CLHEP::eV;                << 
 84 }                                                  89 }
 85                                                    90 
 86 //////////////////////////////////////////////     91 ////////////////////////////////////////////////////////////////////////////
 87                                                    92 
 88 G4PAIPhotModel::~G4PAIPhotModel()                  93 G4PAIPhotModel::~G4PAIPhotModel()
 89 {                                                  94 {
                                                   >>  95   //G4cout << "G4PAIPhotModel::~G4PAIPhotModel() " << this << G4endl;
 90   if(IsMaster()) { delete fModelData; fModelDa     96   if(IsMaster()) { delete fModelData; fModelData = nullptr; }
 91 }                                                  97 }
 92                                                    98 
 93 //////////////////////////////////////////////     99 ////////////////////////////////////////////////////////////////////////////
 94                                                   100 
 95 void G4PAIPhotModel::Initialise(const G4Partic    101 void G4PAIPhotModel::Initialise(const G4ParticleDefinition* p,
 96               const G4DataVector& cuts)        << 102           const G4DataVector& cuts)
 97 {                                                 103 {
 98   if(fVerbose > 1)                             << 104   if(fVerbose > 0) 
 99   {                                               105   {
100     G4cout<<"G4PAIPhotModel::Initialise for "<    106     G4cout<<"G4PAIPhotModel::Initialise for "<<p->GetParticleName()<<G4endl;
101   }                                               107   }
102   SetParticle(p);                                 108   SetParticle(p);
103   fParticleChange = GetParticleChangeForLoss()    109   fParticleChange = GetParticleChangeForLoss();
104                                                   110 
105   if( IsMaster() )                                111   if( IsMaster() ) 
106   {                                               112   { 
                                                   >> 113 
                                                   >> 114     InitialiseElementSelectors(p, cuts);
                                                   >> 115 
107     delete fModelData;                            116     delete fModelData;
108     fMaterialCutsCoupleVector.clear();            117     fMaterialCutsCoupleVector.clear(); 
109                                                   118 
110     G4double tmin = LowEnergyLimit()*fRatio;      119     G4double tmin = LowEnergyLimit()*fRatio;
111     G4double tmax = HighEnergyLimit()*fRatio;     120     G4double tmax = HighEnergyLimit()*fRatio;
112     fModelData = new G4PAIPhotData(tmin, tmax,    121     fModelData = new G4PAIPhotData(tmin, tmax, fVerbose);
113                                                   122     
114     // Prepare initialization                     123     // Prepare initialization
115     const G4MaterialTable* theMaterialTable =     124     const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
116     size_t numOfMat   = G4Material::GetNumberO    125     size_t numOfMat   = G4Material::GetNumberOfMaterials();
117     size_t numRegions = fPAIRegionVector.size(    126     size_t numRegions = fPAIRegionVector.size();
118                                                   127 
119     // protect for unit tests                     128     // protect for unit tests
120     if(0 == numRegions) {                         129     if(0 == numRegions) {
121       G4Exception("G4PAIModel::Initialise()","    130       G4Exception("G4PAIModel::Initialise()","em0106",JustWarning,
122                   "no G4Regions are registered    131                   "no G4Regions are registered for the PAI model - World is used");
123       fPAIRegionVector.push_back(G4RegionStore    132       fPAIRegionVector.push_back(G4RegionStore::GetInstance()
124          ->GetRegion("DefaultRegionForTheWorld    133          ->GetRegion("DefaultRegionForTheWorld", false));
125       numRegions = 1;                             134       numRegions = 1;
126     }                                             135     }
127                                                   136 
128     for( size_t iReg = 0; iReg < numRegions; +    137     for( size_t iReg = 0; iReg < numRegions; ++iReg ) 
129     {                                             138     {
130       const G4Region* curReg = fPAIRegionVecto    139       const G4Region* curReg = fPAIRegionVector[iReg];
131       G4Region* reg = const_cast<G4Region*>(cu    140       G4Region* reg = const_cast<G4Region*>(curReg);
132                                                   141 
133       for(size_t jMat = 0; jMat < numOfMat; ++    142       for(size_t jMat = 0; jMat < numOfMat; ++jMat) 
134       {                                           143       {
135   G4Material* mat = (*theMaterialTable)[jMat];    144   G4Material* mat = (*theMaterialTable)[jMat];
136   const G4MaterialCutsCouple* cutCouple = reg-    145   const G4MaterialCutsCouple* cutCouple = reg->FindCouple(mat);
137   if(nullptr != cutCouple)                     << 146   //G4cout << "Couple <" << fCutCouple << G4endl;
                                                   >> 147   if(cutCouple) 
138         {                                         148         {
139     if(fVerbose > 1)                           << 149     if(fVerbose>0)
140     {                                             150     {
141       G4cout << "Reg <" <<curReg->GetName() <<    151       G4cout << "Reg <" <<curReg->GetName() << ">  mat <" 
142       << mat->GetName() << ">  fCouple= "         152       << mat->GetName() << ">  fCouple= " 
143       << cutCouple << ", idx= " << cutCouple->    153       << cutCouple << ", idx= " << cutCouple->GetIndex()
144       <<"  " << p->GetParticleName()              154       <<"  " << p->GetParticleName() 
145       <<", cuts.size() = " << cuts.size() << G    155       <<", cuts.size() = " << cuts.size() << G4endl;
146     }                                             156     }
147     // check if this couple is not already ini    157     // check if this couple is not already initialized
148                                                   158 
149     size_t n = fMaterialCutsCoupleVector.size(    159     size_t n = fMaterialCutsCoupleVector.size();
150                                                   160 
151     G4bool isnew = true;                          161     G4bool isnew = true;
152     if( 0 < n )                                   162     if( 0 < n ) 
153           {                                       163           {
154       for(size_t i=0; i<fMaterialCutsCoupleVec    164       for(size_t i=0; i<fMaterialCutsCoupleVector.size(); ++i) 
155             {                                     165             {
156         if(cutCouple == fMaterialCutsCoupleVec    166         if(cutCouple == fMaterialCutsCoupleVector[i]) {
157     isnew = false;                                167     isnew = false;
158     break;                                        168     break;
159         }                                         169         }
160       }                                           170       }
161     }                                             171     }
162     // initialise data banks                      172     // initialise data banks
163     if(isnew) {                                   173     if(isnew) {
164       fMaterialCutsCoupleVector.push_back(cutC    174       fMaterialCutsCoupleVector.push_back(cutCouple);
165       G4double deltaCutInKinEnergy = cuts[cutC    175       G4double deltaCutInKinEnergy = cuts[cutCouple->GetIndex()];
166       fModelData->Initialise(cutCouple, deltaC    176       fModelData->Initialise(cutCouple, deltaCutInKinEnergy, this);
167     }                                             177     }
168   }                                               178   }
169       }                                           179       }
170     }                                             180     }
171     InitialiseElementSelectors(p, cuts);       << 
172   }                                               181   }
173 }                                                 182 }
174                                                   183 
175 //////////////////////////////////////////////    184 /////////////////////////////////////////////////////////////////////////
176                                                   185 
177 void G4PAIPhotModel::InitialiseLocal(const G4P    186 void G4PAIPhotModel::InitialiseLocal(const G4ParticleDefinition*, 
178          G4VEmModel* masterModel)                 187          G4VEmModel* masterModel)
179 {                                                 188 {
180   fModelData = static_cast<G4PAIPhotModel*>(ma    189   fModelData = static_cast<G4PAIPhotModel*>(masterModel)->GetPAIPhotData();
181   fMaterialCutsCoupleVector = static_cast<G4PA    190   fMaterialCutsCoupleVector = static_cast<G4PAIPhotModel*>(masterModel)->GetVectorOfCouples();
182   SetElementSelectors( masterModel->GetElement    191   SetElementSelectors( masterModel->GetElementSelectors() );
183 }                                                 192 }
184                                                   193 
185 //////////////////////////////////////////////    194 //////////////////////////////////////////////////////////////////////////////
186                                                   195 
187 G4double G4PAIPhotModel::MinEnergyCut(const G4 << 
188               const G4MaterialCutsCouple*)     << 
189 {                                              << 
190   return fLowestTcut;                          << 
191 }                                              << 
192                                                << 
193 ////////////////////////////////////////////// << 
194                                                << 
195 G4double G4PAIPhotModel::ComputeDEDXPerVolume(    196 G4double G4PAIPhotModel::ComputeDEDXPerVolume(const G4Material*,
196             const G4ParticleDefinition* p,        197             const G4ParticleDefinition* p,
197             G4double kineticEnergy,               198             G4double kineticEnergy,
198             G4double cutEnergy)                   199             G4double cutEnergy)
199 {                                                 200 {
200   G4int coupleIndex = FindCoupleIndex(CurrentC    201   G4int coupleIndex = FindCoupleIndex(CurrentCouple());
201   if(0 > coupleIndex) { return 0.0; }             202   if(0 > coupleIndex) { return 0.0; }
202                                                   203 
203   G4double cut = std::min(MaxSecondaryEnergy(p    204   G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
                                                   >> 205 
204   G4double scaledTkin = kineticEnergy*fRatio;     206   G4double scaledTkin = kineticEnergy*fRatio;
205   G4double dedx = fChargeSquare*fModelData->DE << 207  
206   return dedx;                                 << 208   return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin, cut);
207 }                                                 209 }
208                                                   210 
209 //////////////////////////////////////////////    211 /////////////////////////////////////////////////////////////////////////
210                                                   212 
211 G4double G4PAIPhotModel::CrossSectionPerVolume    213 G4double G4PAIPhotModel::CrossSectionPerVolume( const G4Material*,
212               const G4ParticleDefinition* p,      214               const G4ParticleDefinition* p,
213               G4double kineticEnergy,             215               G4double kineticEnergy,
214               G4double cutEnergy,                 216               G4double cutEnergy,
215               G4double maxEnergy  )               217               G4double maxEnergy  ) 
216 {                                                 218 {
217   G4int coupleIndex = FindCoupleIndex(CurrentC    219   G4int coupleIndex = FindCoupleIndex(CurrentCouple());
218   if(0 > coupleIndex) { return 0.0; }          << 220 
                                                   >> 221   if(0 > coupleIndex) return 0.0; 
219                                                   222 
220   G4double tmax = std::min(MaxSecondaryEnergy(    223   G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
221   if(tmax <= cutEnergy) { return 0.0; }        << 224 
                                                   >> 225   if(tmax <= cutEnergy)  return 0.0; 
222                                                   226 
223   G4double scaledTkin = kineticEnergy*fRatio;     227   G4double scaledTkin = kineticEnergy*fRatio;
224   G4double xs = fChargeSquare*fModelData->Cros << 228   G4double xsc = fChargeSquare*fModelData->CrossSectionPerVolume(coupleIndex, 
225                                           scal << 229                scaledTkin, 
226   return xs;                                   << 230                cutEnergy, tmax);
                                                   >> 231   return xsc;
227 }                                                 232 }
228                                                   233 
229 //////////////////////////////////////////////    234 ///////////////////////////////////////////////////////////////////////////
230 //                                                235 //
231 // It is analog of PostStepDoIt in terms of se    236 // It is analog of PostStepDoIt in terms of secondary electron.
232 //                                                237 //
233                                                   238  
234 void G4PAIPhotModel::SampleSecondaries(std::ve    239 void G4PAIPhotModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
235            const G4MaterialCutsCouple* matCC,     240            const G4MaterialCutsCouple* matCC,
236            const G4DynamicParticle* dp,           241            const G4DynamicParticle* dp,
237            G4double tmin,                         242            G4double tmin,
238            G4double maxEnergy)                    243            G4double maxEnergy)
239 {                                                 244 {
240   G4int coupleIndex = FindCoupleIndex(matCC);     245   G4int coupleIndex = FindCoupleIndex(matCC);
241   if(0 > coupleIndex) { return; }                 246   if(0 > coupleIndex) { return; }
242                                                   247 
243   SetParticle(dp->GetDefinition());               248   SetParticle(dp->GetDefinition());
244                                                   249 
245   G4double kineticEnergy = dp->GetKineticEnerg    250   G4double kineticEnergy = dp->GetKineticEnergy();
246                                                   251 
247   G4double tmax = MaxSecondaryEnergy(fParticle    252   G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy);
248                                                   253 
249   if( maxEnergy <  tmax) tmax = maxEnergy;        254   if( maxEnergy <  tmax) tmax = maxEnergy; 
250   if( tmin      >= tmax) return;                  255   if( tmin      >= tmax) return; 
251                                                   256 
252   G4ThreeVector direction = dp->GetMomentumDir    257   G4ThreeVector direction = dp->GetMomentumDirection();
253   G4double scaledTkin     = kineticEnergy*fRat    258   G4double scaledTkin     = kineticEnergy*fRatio;
254   G4double totalEnergy    = kineticEnergy + fM    259   G4double totalEnergy    = kineticEnergy + fMass;
255   G4double totalMomentum  = sqrt(kineticEnergy    260   G4double totalMomentum  = sqrt(kineticEnergy*(totalEnergy + fMass));
256   G4double plRatio        = fModelData->GetPla    261   G4double plRatio        = fModelData->GetPlasmonRatio(coupleIndex, scaledTkin);
257                                                   262 
258   if( G4UniformRand() <= plRatio )                263   if( G4UniformRand() <= plRatio )
259   {                                               264   {
260     G4double deltaTkin = fModelData->SamplePos    265     G4double deltaTkin = fModelData->SamplePostStepPlasmonTransfer(coupleIndex, scaledTkin);
261                                                   266 
262     // G4cout<<"G4PAIPhotModel::SampleSecondar    267     // G4cout<<"G4PAIPhotModel::SampleSecondaries; dp "<<dp->GetParticleDefinition()->GetParticleName()
263     // <<"; Tkin = "<<kineticEnergy/keV<<" keV    268     // <<"; Tkin = "<<kineticEnergy/keV<<" keV; transfer = "<<deltaTkin/keV<<" keV "<<G4endl;
264                                                   269 
265     if( deltaTkin <= 0. && fVerbose > 0)          270     if( deltaTkin <= 0. && fVerbose > 0) 
266     {                                             271     {
267       G4cout<<"G4PAIPhotModel::SampleSecondary << 272     G4cout<<"G4PAIPhotModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
268     }                                             273     }
269     if( deltaTkin <= 0.) { return; }           << 274     if( deltaTkin <= 0.) return; 
270                                                   275 
271     if( deltaTkin > tmax) { deltaTkin = tmax;  << 276     if( deltaTkin > tmax) deltaTkin = tmax;
272                                                   277 
273     const G4Element* anElement = SelectTargetA << 278     const G4Element* anElement = SelectRandomAtom(matCC,fParticle,kineticEnergy);
274                                                << 279     G4int Z = G4lrint(anElement->GetZ());
275     G4int Z = anElement->GetZasInt();          << 
276                                                   280  
277     auto deltaRay = new G4DynamicParticle(fEle << 281     G4DynamicParticle* deltaRay = new G4DynamicParticle(fElectron, 
278       GetAngularDistribution()->SampleDirectio    282       GetAngularDistribution()->SampleDirection(dp, deltaTkin,
279                   Z, matCC->GetMaterial()),       283                   Z, matCC->GetMaterial()),
280                   deltaTkin);                     284                   deltaTkin);
281                                                   285 
282     // primary change                             286     // primary change
283                                                   287 
284     kineticEnergy -= deltaTkin;                   288     kineticEnergy -= deltaTkin;
285                                                   289 
286     if( kineticEnergy <= 0. ) // kill primary     290     if( kineticEnergy <= 0. ) // kill primary as local? energy deposition
287     {                                             291     {
288       fParticleChange->SetProposedKineticEnerg    292       fParticleChange->SetProposedKineticEnergy(0.0);
289       fParticleChange->ProposeLocalEnergyDepos    293       fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy+deltaTkin);
                                                   >> 294       // fParticleChange->ProposeTrackStatus(fStopAndKill);
290       return;                                     295       return; 
291     }                                             296     }
292     else                                          297     else
293     {                                             298     {
294       G4ThreeVector dir = totalMomentum*direct    299       G4ThreeVector dir = totalMomentum*direction - deltaRay->GetMomentum();
295       direction = dir.unit();                     300       direction = dir.unit();
296       fParticleChange->SetProposedKineticEnerg    301       fParticleChange->SetProposedKineticEnergy(kineticEnergy);
297       fParticleChange->SetProposedMomentumDire    302       fParticleChange->SetProposedMomentumDirection(direction);
298       vdp->push_back(deltaRay);                   303       vdp->push_back(deltaRay);
299     }                                             304     }
300   }                                               305   }
301   else // secondary X-ray CR photon               306   else // secondary X-ray CR photon
302   {                                               307   {
303     G4double deltaTkin     = fModelData->Sampl    308     G4double deltaTkin     = fModelData->SamplePostStepPhotonTransfer(coupleIndex, scaledTkin);
304                                                   309 
305     //  G4cout<<"PAIPhotonModel PhotonPostStep    310     //  G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl; 
306                                                   311 
307     if( deltaTkin <= 0. )                         312     if( deltaTkin <= 0. )
308     {                                             313     {
309       G4cout<<"G4PAIPhotonModel::SampleSeconda    314       G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
310     }                                             315     }
311     if( deltaTkin <= 0.) return;                  316     if( deltaTkin <= 0.) return;
312                                                   317 
313     if( deltaTkin >= kineticEnergy ) // stop p    318     if( deltaTkin >= kineticEnergy ) // stop primary
314     {                                             319     {
315       deltaTkin = kineticEnergy;                  320       deltaTkin = kineticEnergy;
316       kineticEnergy = 0.0;                        321       kineticEnergy = 0.0;
317     }                                             322     }
318     G4double costheta = 0.; // G4UniformRand()    323     G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
319     G4double sintheta = sqrt((1.+costheta)*(1.    324     G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
320                                                   325 
321     //  direction of the 'Cherenkov' photon       326     //  direction of the 'Cherenkov' photon  
322     G4double phi = twopi*G4UniformRand();         327     G4double phi = twopi*G4UniformRand(); 
323     G4double dirx = sintheta*cos(phi), diry =     328     G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
324                                                   329 
325     G4ThreeVector deltaDirection(dirx,diry,dir    330     G4ThreeVector deltaDirection(dirx,diry,dirz);
326     deltaDirection.rotateUz(direction);           331     deltaDirection.rotateUz(direction);
327                                                   332 
328     if( kineticEnergy > 0.) // primary change     333     if( kineticEnergy > 0.) // primary change
329     {                                             334     {
330       kineticEnergy -= deltaTkin;                 335       kineticEnergy -= deltaTkin;
331       fParticleChange->SetProposedKineticEnerg    336       fParticleChange->SetProposedKineticEnergy(kineticEnergy);
332     }                                             337     }
333     else // stop primary, but pass X-ray CR       338     else // stop primary, but pass X-ray CR
334     {                                             339     {
335       // fParticleChange->ProposeLocalEnergyDe    340       // fParticleChange->ProposeLocalEnergyDeposit(deltaTkin);
336       fParticleChange->SetProposedKineticEnerg    341       fParticleChange->SetProposedKineticEnergy(0.0);
337     }                                             342     }
338     // create G4DynamicParticle object for pho    343     // create G4DynamicParticle object for photon ray
339                                                   344  
340     auto photonRay = new G4DynamicParticle;    << 345     G4DynamicParticle* photonRay = new G4DynamicParticle;
341     photonRay->SetDefinition( G4Gamma::Gamma()    346     photonRay->SetDefinition( G4Gamma::Gamma() );
342     photonRay->SetKineticEnergy( deltaTkin );     347     photonRay->SetKineticEnergy( deltaTkin );
343     photonRay->SetMomentumDirection(deltaDirec    348     photonRay->SetMomentumDirection(deltaDirection); 
344                                                   349 
345     vdp->push_back(photonRay);                    350     vdp->push_back(photonRay);
346   }                                               351   }
347   return;                                         352   return;
348 }                                                 353 }
349                                                   354 
350 //////////////////////////////////////////////    355 ///////////////////////////////////////////////////////////////////////
351                                                   356 
352 G4double G4PAIPhotModel::SampleFluctuations(   << 357 G4double G4PAIPhotModel::SampleFluctuations( const G4MaterialCutsCouple* matCC,
353                          const G4MaterialCutsC << 358                                          const G4DynamicParticle* aParticle,
354                          const G4DynamicPartic << 359            G4double, G4double step,
355                          const G4double, const << 360            G4double eloss)
356                          const G4double step,  << 
357 {                                                 361 {
358   // return 0.;                                   362   // return 0.;
359   G4int coupleIndex = FindCoupleIndex(matCC);     363   G4int coupleIndex = FindCoupleIndex(matCC);
360   if(0 > coupleIndex) { return eloss; }           364   if(0 > coupleIndex) { return eloss; }
361                                                   365 
362   SetParticle(aParticle->GetDefinition());        366   SetParticle(aParticle->GetDefinition());
363                                                   367 
                                                   >> 368   
364   // G4cout << "G4PAIPhotModel::SampleFluctuat    369   // G4cout << "G4PAIPhotModel::SampleFluctuations step(mm)= "<< step/mm
365   // << "  Eloss(keV)= " << eloss/keV  << " in    370   // << "  Eloss(keV)= " << eloss/keV  << " in " 
366   // << matCC->GetMaterial()->GetName() << G4e    371   // << matCC->GetMaterial()->GetName() << G4endl;
                                                   >> 372   
367                                                   373 
368   G4double Tkin       = aParticle->GetKineticE    374   G4double Tkin       = aParticle->GetKineticEnergy();
369   G4double scaledTkin = Tkin*fRatio;              375   G4double scaledTkin = Tkin*fRatio;
370                                                   376 
371   G4double loss = fModelData->SampleAlongStepP    377   G4double loss = fModelData->SampleAlongStepPhotonTransfer(coupleIndex, Tkin,
372                         scaledTkin,            << 378                   scaledTkin,
373                         step*fChargeSquare);   << 379                   step*fChargeSquare);
374   loss += fModelData->SampleAlongStepPlasmonTr << 380            loss += fModelData->SampleAlongStepPlasmonTransfer(coupleIndex, Tkin,
375                                                << 381                   scaledTkin,
                                                   >> 382                     step*fChargeSquare);
                                                   >> 383 
376                                                   384   
377   // G4cout<<"  PAIPhotModel::SampleFluctuatio    385   // G4cout<<"  PAIPhotModel::SampleFluctuations loss = "<<loss/keV<<" keV, on step = "
378   // <<step/mm<<" mm"<<G4endl;                    386   // <<step/mm<<" mm"<<G4endl; 
379   return loss;                                    387   return loss;
380                                                   388 
381 }                                                 389 }
382                                                   390 
383 //////////////////////////////////////////////    391 //////////////////////////////////////////////////////////////////////
384 //                                                392 //
385 // Returns the statistical estimation of the e    393 // Returns the statistical estimation of the energy loss distribution variance
386 //                                                394 //
387                                                   395 
388                                                   396 
389 G4double G4PAIPhotModel::Dispersion(const G4Ma << 397 G4double G4PAIPhotModel::Dispersion( const G4Material* material, 
390                                     const G4Dy << 398                                  const G4DynamicParticle* aParticle,
391             const G4double tcut,               << 399                G4double tmax, 
392             const G4double tmax,               << 400                      G4double step       )
393                   const G4double step)         << 
394 {                                                 401 {
395   G4double particleMass  = aParticle->GetMass(    402   G4double particleMass  = aParticle->GetMass();
396   G4double electronDensity = material->GetElec    403   G4double electronDensity = material->GetElectronDensity();
397   G4double kineticEnergy = aParticle->GetKinet    404   G4double kineticEnergy = aParticle->GetKineticEnergy();
398   G4double q = aParticle->GetCharge()/eplus;      405   G4double q = aParticle->GetCharge()/eplus;
399   G4double etot = kineticEnergy + particleMass    406   G4double etot = kineticEnergy + particleMass;
400   G4double beta2 = kineticEnergy*(kineticEnerg    407   G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
401   G4double siga  = (tmax/beta2 - 0.5*tcut) * t << 408   G4double siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
402                  * electronDensity * q * q;       409                  * electronDensity * q * q;
403                                                   410 
404   return siga;                                    411   return siga;
                                                   >> 412   /*
                                                   >> 413   G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
                                                   >> 414   for(G4int i = 0; i < fMeanNumber; i++)
                                                   >> 415   {
                                                   >> 416     loss      = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
                                                   >> 417     sumLoss  += loss;
                                                   >> 418     sumLoss2 += loss*loss;
                                                   >> 419   }
                                                   >> 420   meanLoss = sumLoss/fMeanNumber;
                                                   >> 421   sigma2   = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
                                                   >> 422   return sigma2;
                                                   >> 423   */
405 }                                                 424 }
406                                                   425 
407 //////////////////////////////////////////////    426 /////////////////////////////////////////////////////////////////////
408                                                   427 
409 G4double G4PAIPhotModel::MaxSecondaryEnergy( c    428 G4double G4PAIPhotModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
410            G4double kinEnergy)                    429            G4double kinEnergy) 
411 {                                                 430 {
412   SetParticle(p);                                 431   SetParticle(p);
413   G4double tmax = kinEnergy;                      432   G4double tmax = kinEnergy;
414   if(p == fElectron) { tmax *= 0.5; }             433   if(p == fElectron) { tmax *= 0.5; }
415   else if(p != fPositron) {                       434   else if(p != fPositron) { 
416     G4double ratio= electron_mass_c2/fMass;       435     G4double ratio= electron_mass_c2/fMass;
417     G4double gamma= kinEnergy/fMass + 1.0;        436     G4double gamma= kinEnergy/fMass + 1.0;
418     tmax = 2.0*electron_mass_c2*(gamma*gamma -    437     tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
419                   (1. + 2.0*gamma*ratio + rati    438                   (1. + 2.0*gamma*ratio + ratio*ratio);
420   }                                               439   }
421   return tmax;                                    440   return tmax;
422 }                                                 441 }
423                                                   442 
424 //////////////////////////////////////////////    443 ///////////////////////////////////////////////////////////////
425                                                   444 
426 void G4PAIPhotModel::DefineForRegion(const G4R    445 void G4PAIPhotModel::DefineForRegion(const G4Region* r) 
427 {                                                 446 {
428   fPAIRegionVector.push_back(r);                  447   fPAIRegionVector.push_back(r);
429 }                                                 448 }
430                                                   449 
431 //                                                450 //
432 //                                                451 //
433 //////////////////////////////////////////////    452 /////////////////////////////////////////////////
434                                                   453