| Geant4 Cross Reference |
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: G4PAIModel.cc,v 1.55 2010-11-21 10:55:44 grichine Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class 31 // GEANT4 Class
30 // File name: G4PAIModel.cc 32 // File name: G4PAIModel.cc
31 // 33 //
32 // Author: Vladimir.Grichine@cern.ch on base o << 34 // Author: Vladimir.Grichine@cern.ch on base of Vladimir Ivanchenko model interface
33 // 35 //
34 // Creation date: 05.10.2003 36 // Creation date: 05.10.2003
35 // 37 //
36 // Modifications: 38 // Modifications:
37 // 39 //
38 // 17.08.04 V.Grichine, bug fixed for Tkin<=0 40 // 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary
39 // 16.08.04 V.Grichine, bug fixed in massRatio << 41 // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary
40 // SampleSecondary <<
41 // 08.04.05 Major optimisation of internal int 42 // 08.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
42 // 26.07.09 Fixed logic to work with several m 43 // 26.07.09 Fixed logic to work with several materials (V.Ivantchenko)
43 // 21.11.10 V. Grichine verbose flag for proto << 44 // 21.11.10 V. Grichine verbose flag for protons and G4PAYySection to check sandia table
44 // check sandia table <<
45 // 12.06.13 V. Grichine Bug fixed in SampleSec <<
46 // (fMass -> proton_mass_c2) <<
47 // 19.08.13 V.Ivanchenko extract data handling <<
48 // added sharing of internal data bet <<
49 // 45 //
50 46
51 #include "G4PAIModel.hh" <<
52 <<
53 #include "G4SystemOfUnits.hh" <<
54 #include "G4PhysicalConstants.hh" <<
55 #include "G4Region.hh" 47 #include "G4Region.hh"
>> 48 #include "G4PhysicsLogVector.hh"
>> 49 #include "G4PhysicsFreeVector.hh"
>> 50 #include "G4PhysicsTable.hh"
>> 51 #include "G4ProductionCutsTable.hh"
56 #include "G4MaterialCutsCouple.hh" 52 #include "G4MaterialCutsCouple.hh"
57 #include "G4MaterialTable.hh" 53 #include "G4MaterialTable.hh"
58 #include "G4RegionStore.hh" << 54 #include "G4SandiaTable.hh"
>> 55 #include "G4OrderedTable.hh"
59 56
>> 57 #include "G4PAIModel.hh"
60 #include "Randomize.hh" 58 #include "Randomize.hh"
61 #include "G4Electron.hh" 59 #include "G4Electron.hh"
62 #include "G4Positron.hh" 60 #include "G4Positron.hh"
63 #include "G4Poisson.hh" 61 #include "G4Poisson.hh"
64 #include "G4Step.hh" 62 #include "G4Step.hh"
65 #include "G4Material.hh" 63 #include "G4Material.hh"
66 #include "G4DynamicParticle.hh" 64 #include "G4DynamicParticle.hh"
67 #include "G4ParticleDefinition.hh" 65 #include "G4ParticleDefinition.hh"
68 #include "G4ParticleChangeForLoss.hh" 66 #include "G4ParticleChangeForLoss.hh"
69 #include "G4PAIModelData.hh" <<
70 #include "G4DeltaAngle.hh" <<
71 67
72 ////////////////////////////////////////////// 68 ////////////////////////////////////////////////////////////////////////
73 69
74 using namespace std; 70 using namespace std;
75 71
76 G4PAIModel::G4PAIModel(const G4ParticleDefinit 72 G4PAIModel::G4PAIModel(const G4ParticleDefinition* p, const G4String& nam)
77 : G4VEmModel(nam),G4VEmFluctuationModel(nam) 73 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
78 fVerbose(0), << 74 fVerbose(0),
79 fModelData(nullptr), << 75 fLowestGamma(1.005),
80 fParticle(nullptr) << 76 fHighestGamma(10000.),
>> 77 fTotBin(200),
>> 78 fMeanNumber(20),
>> 79 fParticle(0),
>> 80 fHighKinEnergy(100.*TeV),
>> 81 fTwoln10(2.0*log(10.0)),
>> 82 fBg2lim(0.0169),
>> 83 fTaulim(8.4146e-3)
81 { 84 {
82 fElectron = G4Electron::Electron(); 85 fElectron = G4Electron::Electron();
83 fPositron = G4Positron::Positron(); 86 fPositron = G4Positron::Positron();
84 87
85 fParticleChange = nullptr; << 88 fPAItransferTable = 0;
>> 89 fPAIdEdxTable = 0;
>> 90 fSandiaPhotoAbsCof = 0;
>> 91 fdEdxVector = 0;
>> 92 fLambdaVector = 0;
>> 93 fdNdxCutVector = 0;
>> 94 fParticleEnergyVector = 0;
>> 95 fSandiaIntervalNumber = 0;
>> 96 fMatIndex = 0;
>> 97 fDeltaCutInKinEnergy = 0.0;
>> 98 fParticleChange = 0;
>> 99 fMaterial = 0;
>> 100 fCutCouple = 0;
86 101
87 if(p) { SetParticle(p); } 102 if(p) { SetParticle(p); }
88 else { SetParticle(fElectron); } 103 else { SetParticle(fElectron); }
89 104
90 // default generator << 105 isInitialised = false;
91 SetAngularDistribution(new G4DeltaAngle()); <<
92 fLowestTcut = 12.5*CLHEP::eV; <<
93 } 106 }
94 107
95 ////////////////////////////////////////////// 108 ////////////////////////////////////////////////////////////////////////////
96 109
97 G4PAIModel::~G4PAIModel() 110 G4PAIModel::~G4PAIModel()
98 { 111 {
99 if(IsMaster()) { delete fModelData; } << 112 // G4cout << "PAI: start destruction" << G4endl;
>> 113 if(fParticleEnergyVector) delete fParticleEnergyVector;
>> 114 if(fdEdxVector) delete fdEdxVector ;
>> 115 if(fLambdaVector) delete fLambdaVector;
>> 116 if(fdNdxCutVector) delete fdNdxCutVector;
>> 117
>> 118 if( fPAItransferTable )
>> 119 {
>> 120 fPAItransferTable->clearAndDestroy();
>> 121 delete fPAItransferTable ;
>> 122 }
>> 123 if( fPAIdEdxTable )
>> 124 {
>> 125 fPAIdEdxTable->clearAndDestroy();
>> 126 delete fPAIdEdxTable ;
>> 127 }
>> 128 if(fSandiaPhotoAbsCof)
>> 129 {
>> 130 for(G4int i=0;i<fSandiaIntervalNumber;i++)
>> 131 {
>> 132 delete[] fSandiaPhotoAbsCof[i];
>> 133 }
>> 134 delete[] fSandiaPhotoAbsCof;
>> 135 }
>> 136 //G4cout << "PAI: end destruction" << G4endl;
>> 137 }
>> 138
>> 139 ///////////////////////////////////////////////////////////////////////////////
>> 140
>> 141 void G4PAIModel::SetParticle(const G4ParticleDefinition* p)
>> 142 {
>> 143 if(fParticle == p) { return; }
>> 144 fParticle = p;
>> 145 fMass = fParticle->GetPDGMass();
>> 146 fSpin = fParticle->GetPDGSpin();
>> 147 G4double q = fParticle->GetPDGCharge()/eplus;
>> 148 fChargeSquare = q*q;
>> 149 fLowKinEnergy = 0.2*MeV*fMass/proton_mass_c2;
>> 150 fRatio = electron_mass_c2/fMass;
>> 151 fQc = fMass/fRatio;
>> 152 fLowestKineticEnergy = fMass*(fLowestGamma - 1.0);
>> 153 fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
100 } 154 }
101 155
102 ////////////////////////////////////////////// 156 ////////////////////////////////////////////////////////////////////////////
103 157
104 void G4PAIModel::Initialise(const G4ParticleDe 158 void G4PAIModel::Initialise(const G4ParticleDefinition* p,
105 const G4DataVector& cuts) << 159 const G4DataVector&)
106 { 160 {
107 if(fVerbose > 1) { << 161 if( fVerbose > 0 && p->GetParticleName()=="proton")
>> 162 {
108 G4cout<<"G4PAIModel::Initialise for "<<p-> 163 G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
>> 164 fPAIySection.SetVerbose(1);
109 } 165 }
>> 166 else fPAIySection.SetVerbose(0);
>> 167
>> 168 if(isInitialised) { return; }
>> 169 isInitialised = true;
>> 170
110 SetParticle(p); 171 SetParticle(p);
111 fParticleChange = GetParticleChangeForLoss() <<
112 172
113 if(IsMaster()) { << 173 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
>> 174 fHighestKineticEnergy,
>> 175 fTotBin);
114 176
115 delete fModelData; << 177 fParticleChange = GetParticleChangeForLoss();
116 fMaterialCutsCoupleVector.clear(); << 178
117 if(fVerbose > 1) { << 179 // Prepare initialization
118 G4cout << "G4PAIModel instantiates data << 180 const G4ProductionCutsTable* theCoupleTable =
119 << G4endl; << 181 G4ProductionCutsTable::GetProductionCutsTable();
120 } << 182 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
121 G4double tmin = LowEnergyLimit()*fRatio; << 183 size_t numOfMat = G4Material::GetNumberOfMaterials();
122 G4double tmax = HighEnergyLimit()*fRatio; << 184 size_t numRegions = fPAIRegionVector.size();
123 fModelData = new G4PAIModelData(tmin, tmax << 185
124 << 186 for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
125 // Prepare initialization << 187 {
126 const G4MaterialTable* theMaterialTable = << 188 const G4Region* curReg = fPAIRegionVector[iReg];
127 std::size_t numOfMat = G4Material::GetNu << 189
128 std::size_t numRegions = fPAIRegionVector. << 190 for(size_t jMat = 0; jMat < numOfMat; ++jMat) // region material loop
129 << 191 {
130 // protect for unit tests << 192 fMaterial = (*theMaterialTable)[jMat];
131 if(0 == numRegions) { << 193 fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial,
132 G4Exception("G4PAIModel::Initialise()"," << 194 curReg->GetProductionCuts() );
133 "no G4Regions are registered << 195 //G4cout << "Reg <" <<curReg->GetName() << "> mat <"
134 fPAIRegionVector.push_back(G4RegionStore << 196 // << fMaterial->GetName() << "> fCouple= "
135 ->GetRegion("DefaultRegionForTheWorld << 197 // << fCutCouple<<" " << p->GetParticleName() <<G4endl;
136 numRegions = 1; << 198 if( fCutCouple )
137 } << 199 {
138 << 200 fMaterialCutsCoupleVector.push_back(fCutCouple);
139 if(fVerbose > 1) { << 201
140 G4cout << "G4PAIModel is defined for " < << 202 fPAItransferTable = new G4PhysicsTable(fTotBin+1);
141 << "; number of materials " << numOfMat << 203 fPAIdEdxTable = new G4PhysicsTable(fTotBin+1);
142 } << 204
143 for(std::size_t iReg = 0; iReg<numRegions; << 205 fDeltaCutInKinEnergy =
144 const G4Region* curReg = fPAIRegionVecto << 206 (*theCoupleTable->GetEnergyCutsVector(1))[fCutCouple->GetIndex()];
145 G4Region* reg = const_cast<G4Region*>(cu << 207
146 << 208 //ComputeSandiaPhotoAbsCof();
147 for(std::size_t jMat = 0; jMat<numOfMat; << 209 BuildPAIonisationTable();
148 G4Material* mat = (*theMaterialTable)[jMat]; << 210
149 const G4MaterialCutsCouple* cutCouple = reg- << 211 fPAIxscBank.push_back(fPAItransferTable);
150 std::size_t n = fMaterialCutsCoupleVector.si << 212 fPAIdEdxBank.push_back(fPAIdEdxTable);
151 if(nullptr != cutCouple) { << 213 fdEdxTable.push_back(fdEdxVector);
152 if(fVerbose > 1) { << 214
153 G4cout << "Region <" << curReg->GetName( << 215 BuildLambdaVector();
154 << mat->GetName() << "> CoupleIndex= " << 216 fdNdxCutTable.push_back(fdNdxCutVector);
155 << cutCouple->GetIndex() << 217 fLambdaTable.push_back(fLambdaVector);
156 << " " << p->GetParticleName() <<
157 << " cutsize= " << cuts.size() << G4end <<
158 } <<
159 // check if this couple is not already ini <<
160 G4bool isnew = true; <<
161 if(0 < n) { <<
162 for(std::size_t i=0; i<n; ++i) { <<
163 G4cout << i << G4endl; <<
164 if(cutCouple == fMaterialCutsCoupleVec <<
165 isnew = false; <<
166 break; <<
167 } <<
168 } <<
169 } <<
170 // initialise data banks <<
171 // G4cout << " isNew: " << isnew << " " <<
172 if(isnew) { <<
173 fMaterialCutsCoupleVector.push_back(cutC <<
174 fModelData->Initialise(cutCouple, this); <<
175 } <<
176 } <<
177 } 218 }
178 } 219 }
179 InitialiseElementSelectors(p, cuts); <<
180 } 220 }
181 } 221 }
182 222
183 ////////////////////////////////////////////// << 223 //////////////////////////////////////////////////////////////////
184 224
185 void G4PAIModel::InitialiseLocal(const G4Parti << 225 void G4PAIModel::InitialiseMe(const G4ParticleDefinition*)
186 G4VEmModel* masterModel) << 226 {}
>> 227
>> 228 //////////////////////////////////////////////////////////////////
>> 229
>> 230 void G4PAIModel::ComputeSandiaPhotoAbsCof()
>> 231 {
>> 232 G4int i, j;
>> 233 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 234
>> 235 G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
>> 236 G4int numberOfElements =
>> 237 (*theMaterialTable)[fMatIndex]->GetNumberOfElements();
>> 238
>> 239 G4int* thisMaterialZ = new G4int[numberOfElements] ;
>> 240
>> 241 for(i=0;i<numberOfElements;i++)
>> 242 {
>> 243 thisMaterialZ[i] =
>> 244 (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
>> 245 }
>> 246 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
>> 247 (thisMaterialZ,numberOfElements) ;
>> 248
>> 249 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
>> 250 ( thisMaterialZ ,
>> 251 (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
>> 252 numberOfElements,fSandiaIntervalNumber) ;
>> 253
>> 254 fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
>> 255
>> 256 for(i=0; i<fSandiaIntervalNumber; i++)
>> 257 {
>> 258 fSandiaPhotoAbsCof[i] = new G4double[5];
>> 259 }
>> 260
>> 261 for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
>> 262 {
>> 263 fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0);
>> 264
>> 265 for( j = 1; j < 5 ; j++ )
>> 266 {
>> 267 fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,j)*
>> 268 (*theMaterialTable)[fMatIndex]->GetDensity() ;
>> 269 }
>> 270 }
>> 271 delete[] thisMaterialZ;
>> 272 }
>> 273
>> 274 ////////////////////////////////////////////////////////////////////////////
>> 275 //
>> 276 // Build tables for the ionization energy loss
>> 277 // the tables are built for MATERIALS
>> 278 // *********
>> 279
>> 280 void G4PAIModel::BuildPAIonisationTable()
187 { 281 {
188 fModelData = static_cast<G4PAIModel*>(master << 282 G4double LowEdgeEnergy , ionloss ;
189 fMaterialCutsCoupleVector = << 283 G4double tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2 ;
190 static_cast<G4PAIModel*>(masterModel)->Get << 284
191 SetElementSelectors(masterModel->GetElementS << 285 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
>> 286 fHighestKineticEnergy,
>> 287 fTotBin);
>> 288 G4SandiaTable* sandia = fMaterial->GetSandiaTable();
>> 289
>> 290 Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV;
>> 291 deltaLow = 100.*eV; // 0.5*eV ;
>> 292
>> 293 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy
>> 294 {
>> 295 LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
>> 296 tau = LowEdgeEnergy/fMass ;
>> 297 //gamma = tau +1. ;
>> 298 // G4cout<<"gamma = "<<gamma<<endl ;
>> 299 bg2 = tau*( tau + 2. );
>> 300 Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
>> 301 // Tmax = std::min(fDeltaCutInKinEnergy, Tmax);
>> 302 Tkin = Tmax ;
>> 303
>> 304 if ( Tmax < Tmin + deltaLow ) // low energy safety
>> 305 Tkin = Tmin + deltaLow ;
>> 306
>> 307 fPAIySection.Initialize(fMaterial, Tkin, bg2);
>> 308
>> 309 // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
>> 310 // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
>> 311 // G4cout<<"protonPAI.GetSplineSize() = "<<
>> 312 // protonPAI.GetSplineSize()<<G4endl<<G4endl ;
>> 313
>> 314 G4int n = fPAIySection.GetSplineSize();
>> 315 G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n) ;
>> 316 G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
>> 317
>> 318 for( G4int k = 0 ; k < n; k++ )
>> 319 {
>> 320 transferVector->PutValue( k ,
>> 321 fPAIySection.GetSplineEnergy(k+1),
>> 322 fPAIySection.GetIntegralPAIySection(k+1) ) ;
>> 323 dEdxVector->PutValue( k ,
>> 324 fPAIySection.GetSplineEnergy(k+1),
>> 325 fPAIySection.GetIntegralPAIdEdx(k+1) ) ;
>> 326 }
>> 327 ionloss = fPAIySection.GetMeanEnergyLoss() ; // total <dE/dx>
>> 328
>> 329 if ( ionloss < DBL_MIN) { ionloss = 0.0; }
>> 330 fdEdxVector->PutValue(i,ionloss) ;
>> 331
>> 332 fPAItransferTable->insertAt(i,transferVector) ;
>> 333 fPAIdEdxTable->insertAt(i,dEdxVector) ;
>> 334
>> 335 } // end of Tkin loop
192 } 336 }
193 337
194 ////////////////////////////////////////////// << 338 ///////////////////////////////////////////////////////////////////////
>> 339 //
>> 340 // Build mean free path tables for the delta ray production process
>> 341 // tables are built for MATERIALS
>> 342 //
195 343
196 G4double G4PAIModel::MinEnergyCut(const G4Part << 344 void G4PAIModel::BuildLambdaVector()
197 const G4MaterialCutsCouple*) <<
198 { 345 {
199 return fLowestTcut; << 346 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
>> 347 fHighestKineticEnergy,
>> 348 fTotBin ) ;
>> 349 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
>> 350 fHighestKineticEnergy,
>> 351 fTotBin ) ;
>> 352 if(fVerbose > 1)
>> 353 {
>> 354 G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
>> 355 <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
>> 356 }
>> 357 for (G4int i = 0 ; i <= fTotBin ; i++ )
>> 358 {
>> 359 G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
>> 360 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
>> 361 if(dNdxCut < 0.0) dNdxCut = 0.0;
>> 362 // G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
>> 363 G4double lambda = DBL_MAX;
>> 364 if(dNdxCut > 0.0) lambda = 1.0/dNdxCut;
>> 365
>> 366 fLambdaVector->PutValue(i, lambda) ;
>> 367 fdNdxCutVector->PutValue(i, dNdxCut) ;
>> 368 }
>> 369 }
>> 370
>> 371 ///////////////////////////////////////////////////////////////////////
>> 372 //
>> 373 // Returns integral PAI cross section for energy transfers >= transferCut
>> 374
>> 375 G4double
>> 376 G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
>> 377 {
>> 378 G4int iTransfer;
>> 379 G4double x1, x2, y1, y2, dNdxCut;
>> 380 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
>> 381 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
>> 382 // <<G4endl;
>> 383 for( iTransfer = 0 ;
>> 384 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
>> 385 iTransfer++)
>> 386 {
>> 387 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
>> 388 {
>> 389 break ;
>> 390 }
>> 391 }
>> 392 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
>> 393 {
>> 394 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
>> 395 }
>> 396 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
>> 397 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
>> 398 if(y1 < 0.0) y1 = 0.0;
>> 399 if(y2 < 0.0) y2 = 0.0;
>> 400 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
>> 401 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
>> 402 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
>> 403 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
>> 404
>> 405 if ( y1 == y2 ) dNdxCut = y2 ;
>> 406 else
>> 407 {
>> 408 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 409 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 410 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
>> 411 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
>> 412 }
>> 413 // G4cout<<""<<dNdxCut<<G4endl;
>> 414 if(dNdxCut < 0.0) dNdxCut = 0.0;
>> 415 return dNdxCut ;
>> 416 }
>> 417 ///////////////////////////////////////////////////////////////////////
>> 418 //
>> 419 // Returns integral dEdx for energy transfers >= transferCut
>> 420
>> 421 G4double
>> 422 G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
>> 423 {
>> 424 G4int iTransfer;
>> 425 G4double x1, x2, y1, y2, dEdxCut;
>> 426 //G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
>> 427 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
>> 428 // <<G4endl;
>> 429 for( iTransfer = 0 ;
>> 430 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
>> 431 iTransfer++)
>> 432 {
>> 433 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
>> 434 {
>> 435 break ;
>> 436 }
>> 437 }
>> 438 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
>> 439 {
>> 440 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
>> 441 }
>> 442 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
>> 443 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
>> 444 if(y1 < 0.0) y1 = 0.0;
>> 445 if(y2 < 0.0) y2 = 0.0;
>> 446 //G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
>> 447 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
>> 448 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
>> 449 //G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
>> 450
>> 451 if ( y1 == y2 ) dEdxCut = y2 ;
>> 452 else
>> 453 {
>> 454 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 455 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 456 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
>> 457 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
>> 458 }
>> 459 //G4cout<<""<<dEdxCut<<G4endl;
>> 460 if(dEdxCut < 0.0) dEdxCut = 0.0;
>> 461 return dEdxCut ;
200 } 462 }
201 463
202 ////////////////////////////////////////////// 464 //////////////////////////////////////////////////////////////////////////////
203 465
204 G4double G4PAIModel::ComputeDEDXPerVolume(cons 466 G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
205 const G4ParticleDefinition* p, 467 const G4ParticleDefinition* p,
206 G4double kineticEnergy, 468 G4double kineticEnergy,
207 G4double cutEnergy) 469 G4double cutEnergy)
208 { << 470 {
209 G4int coupleIndex = FindCoupleIndex(CurrentC << 471 G4int iTkin,iPlace;
210 if(0 > coupleIndex) { return 0.0; } << 472
>> 473 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
>> 474 G4double cut = cutEnergy;
211 475
212 G4double cut = std::min(MaxSecondaryEnergy(p << 476 G4double massRatio = fMass/p->GetPDGMass();
213 G4double scaledTkin = kineticEnergy*fRatio; << 477 G4double scaledTkin = kineticEnergy*massRatio;
214 G4double dedx = fChargeSquare*fModelData->DE << 478 G4double charge = p->GetPDGCharge();
215 return dedx; << 479 G4double charge2 = charge*charge;
>> 480 const G4MaterialCutsCouple* matCC = CurrentCouple();
>> 481
>> 482 size_t jMat = 0;
>> 483 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
>> 484 {
>> 485 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
>> 486 }
>> 487 //G4cout << "jMat= " << jMat
>> 488 // << " jMax= " << fMaterialCutsCoupleVector.size()
>> 489 // << " matCC: " << matCC;
>> 490 //if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
>> 491 // G4cout << G4endl;
>> 492 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
>> 493
>> 494 fPAIdEdxTable = fPAIdEdxBank[jMat];
>> 495 fdEdxVector = fdEdxTable[jMat];
>> 496 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
>> 497 {
>> 498 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 499 }
>> 500 //G4cout << "E= " << scaledTkin << " iTkin= " << iTkin
>> 501 // << " " << fdEdxVector->GetVectorLength()<<G4endl;
>> 502 iPlace = iTkin - 1;
>> 503 if(iPlace < 0) iPlace = 0;
>> 504 else if(iPlace >= fTotBin) iPlace = fTotBin-1;
>> 505 G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
>> 506 if( dEdx < 0.) dEdx = 0.;
>> 507 return dEdx;
216 } 508 }
217 509
218 ////////////////////////////////////////////// 510 /////////////////////////////////////////////////////////////////////////
219 511
220 G4double G4PAIModel::CrossSectionPerVolume( co 512 G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
221 const G4ParticleDefinition* p, 513 const G4ParticleDefinition* p,
222 G4double kineticEnergy, 514 G4double kineticEnergy,
223 G4double cutEnergy, 515 G4double cutEnergy,
224 G4double maxEnergy ) 516 G4double maxEnergy )
225 { 517 {
226 G4int coupleIndex = FindCoupleIndex(CurrentC << 518 G4int iTkin,iPlace;
227 if(0 > coupleIndex) { return 0.0; } <<
228 <<
229 G4double tmax = std::min(MaxSecondaryEnergy( 519 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
230 if(tmax <= cutEnergy) { return 0.0; } << 520 if(tmax <= cutEnergy) return 0.0;
>> 521 G4double massRatio = fMass/p->GetPDGMass();
>> 522 G4double scaledTkin = kineticEnergy*massRatio;
>> 523 G4double charge = p->GetPDGCharge();
>> 524 G4double charge2 = charge*charge, cross, cross1, cross2;
>> 525 const G4MaterialCutsCouple* matCC = CurrentCouple();
>> 526
>> 527 size_t jMat = 0;
>> 528 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
>> 529 {
>> 530 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
>> 531 }
>> 532 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
>> 533
>> 534 fPAItransferTable = fPAIxscBank[jMat];
>> 535
>> 536 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
>> 537 {
>> 538 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 539 }
>> 540 iPlace = iTkin - 1;
>> 541 if(iPlace < 0) iPlace = 0;
>> 542 else if(iPlace >= fTotBin) iPlace = fTotBin-1;
>> 543
>> 544 //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
>> 545 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
>> 546 cross1 = GetdNdxCut(iPlace,tmax) ;
>> 547 //G4cout<<"cross1 = "<<cross1<<G4endl;
>> 548 cross2 = GetdNdxCut(iPlace,cutEnergy) ;
>> 549 //G4cout<<"cross2 = "<<cross2<<G4endl;
>> 550 cross = (cross2-cross1)*charge2;
>> 551 //G4cout<<"cross = "<<cross<<G4endl;
>> 552 if( cross < DBL_MIN) cross = 0.0;
>> 553 // if( cross2 < DBL_MIN) cross2 = DBL_MIN;
231 554
232 G4double scaledTkin = kineticEnergy*fRatio; << 555 // return cross2;
233 G4double xs = fChargeSquare*fModelData->Cros << 556 return cross;
234 scal <<
235 return xs; <<
236 } 557 }
237 558
238 ////////////////////////////////////////////// 559 ///////////////////////////////////////////////////////////////////////////
239 // 560 //
240 // It is analog of PostStepDoIt in terms of se 561 // It is analog of PostStepDoIt in terms of secondary electron.
241 // 562 //
242 563
243 void G4PAIModel::SampleSecondaries(std::vector 564 void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
244 const G4MaterialCutsCouple* matCC, 565 const G4MaterialCutsCouple* matCC,
245 const G4DynamicParticle* dp, 566 const G4DynamicParticle* dp,
246 G4double tmin, 567 G4double tmin,
247 G4double maxEnergy) 568 G4double maxEnergy)
248 { 569 {
249 G4int coupleIndex = FindCoupleIndex(matCC); << 570 size_t jMat = 0;
250 << 571 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
251 //G4cout << "G4PAIModel::SampleSecondaries: << 572 {
252 if(0 > coupleIndex) { return; } << 573 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
253 << 574 }
254 SetParticle(dp->GetDefinition()); << 575 if(jMat == fMaterialCutsCoupleVector.size()) return;
255 G4double kineticEnergy = dp->GetKineticEnerg <<
256 576
257 G4double tmax = MaxSecondaryEnergy(fParticle << 577 fPAItransferTable = fPAIxscBank[jMat];
258 if(maxEnergy < tmax) { tmax = maxEnergy; } << 578 fdNdxCutVector = fdNdxCutTable[jMat];
259 if(tmin >= tmax) { return; } <<
260 579
>> 580 G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy);
>> 581 if( tmin >= tmax && fVerbose > 0)
>> 582 {
>> 583 G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
>> 584 }
261 G4ThreeVector direction= dp->GetMomentumDire 585 G4ThreeVector direction= dp->GetMomentumDirection();
262 G4double scaledTkin = kineticEnergy*fRati << 586 G4double particleMass = dp->GetMass();
263 G4double totalEnergy = kineticEnergy + fMa << 587 G4double kineticEnergy = dp->GetKineticEnergy();
264 G4double totalMomentum = sqrt(kineticEnergy* <<
265 588
266 G4double deltaTkin = << 589 G4double massRatio = fMass/particleMass;
267 fModelData->SamplePostStepTransfer(coupleI << 590 G4double scaledTkin = kineticEnergy*massRatio;
>> 591 G4double totalEnergy = kineticEnergy + particleMass;
>> 592 G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
>> 593
>> 594 G4double deltaTkin = GetPostStepTransfer(scaledTkin);
268 595
269 //G4cout<<"G4PAIModel::SampleSecondaries; de << 596 // G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl;
270 // <<" keV "<< " Escaled(MeV)= " << scaledT <<
271 597
272 if( !(deltaTkin <= 0.) && !(deltaTkin > 0)) << 598 if( deltaTkin <= 0. && fVerbose > 0)
273 G4cout<<"G4PAIModel::SampleSecondaries; de << 599 {
274 <<" keV "<< " Escaled(MeV)= " << scaledTki << 600 G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
275 return; <<
276 } 601 }
277 if( deltaTkin <= 0.) { return; } << 602 if( deltaTkin <= 0.) return;
278 603
279 if( deltaTkin > tmax) { deltaTkin = tmax; } << 604 if( deltaTkin > tmax) deltaTkin = tmax;
280 605
281 const G4Element* anElement = SelectTargetAto << 606 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
282 << 607 G4double totalMomentum = sqrt(pSquare);
283 << 608 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
284 G4int Z = anElement->GetZasInt(); << 609 /(deltaTotalMomentum * totalMomentum);
285 << 610
286 auto deltaRay = new G4DynamicParticle(fElect << 611 if( costheta > 0.99999 ) costheta = 0.99999;
287 GetAngularDistribution()->SampleDirectio << 612 G4double sintheta = 0.0;
288 Z, matCC->GetMaterial()), << 613 G4double sin2 = 1. - costheta*costheta;
289 deltaTkin); << 614 if( sin2 > 0.) sintheta = sqrt(sin2);
>> 615
>> 616 // direction of the delta electron
>> 617 G4double phi = twopi*G4UniformRand();
>> 618 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
>> 619
>> 620 G4ThreeVector deltaDirection(dirx,diry,dirz);
>> 621 deltaDirection.rotateUz(direction);
>> 622 deltaDirection.unit();
290 623
291 // primary change 624 // primary change
292 kineticEnergy -= deltaTkin; 625 kineticEnergy -= deltaTkin;
293 G4ThreeVector dir = totalMomentum*direction << 626 G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
294 direction = dir.unit(); 627 direction = dir.unit();
295 fParticleChange->SetProposedKineticEnergy(ki 628 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
296 fParticleChange->SetProposedMomentumDirectio 629 fParticleChange->SetProposedMomentumDirection(direction);
297 630
>> 631 // create G4DynamicParticle object for e- delta ray
>> 632 G4DynamicParticle* deltaRay = new G4DynamicParticle;
>> 633 deltaRay->SetDefinition(G4Electron::Electron());
>> 634 deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ???
>> 635 deltaRay->SetMomentumDirection(deltaDirection);
>> 636
298 vdp->push_back(deltaRay); 637 vdp->push_back(deltaRay);
299 } 638 }
300 639
>> 640
301 ////////////////////////////////////////////// 641 ///////////////////////////////////////////////////////////////////////
>> 642 //
>> 643 // Returns post step PAI energy transfer > cut electron energy according to passed
>> 644 // scaled kinetic energy of particle
302 645
303 G4double G4PAIModel::SampleFluctuations(const << 646 G4double
304 const << 647 G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
305 const << 648 {
306 const G4double, << 649 // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
307 const << 650
308 const G4double eloss) << 651 G4int iTkin, iTransfer, iPlace ;
309 { << 652 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
310 G4int coupleIndex = FindCoupleIndex(matCC); << 653
311 if(0 > coupleIndex) { return eloss; } << 654 for(iTkin=0;iTkin<=fTotBin;iTkin++)
312 << 655 {
313 SetParticle(aParticle->GetDefinition()); << 656 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
314 << 657 }
315 /* << 658 iPlace = iTkin - 1 ;
316 G4cout << "G4PAIModel::SampleFluctuations st << 659 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
317 << " Eloss(keV)= " << eloss/keV << " in " << 660 if(iPlace < 0) iPlace = 0;
318 << matCC->Getmaterial()->GetName() << G4end << 661 else if(iPlace >= fTotBin) iPlace = fTotBin-1;
319 */ << 662 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
320 << 663 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
321 G4double Tkin = aParticle->GetKineticE << 664
322 G4double scaledTkin = Tkin*fRatio; << 665
323 << 666 if(iTkin == fTotBin) // Fermi plato, try from left
324 G4double loss = fModelData->SampleAlongStepT << 667 {
325 scaledTkin, tcut, << 668 position = dNdxCut1*G4UniformRand() ;
326 step*fChargeSquare); << 669
327 << 670 for( iTransfer = 0;
328 // G4cout<<"PAIModel AlongStepLoss = "<<loss << 671 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
329 //<<step/mm<<" mm"<<G4endl; << 672 {
330 return loss; << 673 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
>> 674 }
>> 675 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
>> 676 }
>> 677 else
>> 678 {
>> 679 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
>> 680 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
>> 681 if(iTkin == 0) // Tkin is too small, trying from right only
>> 682 {
>> 683 position = dNdxCut2*G4UniformRand() ;
>> 684
>> 685 for( iTransfer = 0;
>> 686 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
>> 687 {
>> 688 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
>> 689 }
>> 690 transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
>> 691 }
>> 692 else // general case: Tkin between two vectors of the material
>> 693 {
>> 694 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
>> 695 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
>> 696 W = 1.0/(E2 - E1) ;
>> 697 W1 = (E2 - scaledTkin)*W ;
>> 698 W2 = (scaledTkin - E1)*W ;
>> 699
>> 700 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
>> 701
>> 702 // G4cout<<position<<"\t" ;
>> 703
>> 704 G4int iTrMax1, iTrMax2, iTrMax;
>> 705
>> 706 iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
>> 707 iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength());
>> 708
>> 709 if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
>> 710 else iTrMax = iTrMax1;
>> 711
>> 712
>> 713 for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
>> 714 {
>> 715 if( position >=
>> 716 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
>> 717 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
>> 718 }
>> 719 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
>> 720 }
>> 721 }
>> 722 // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
>> 723 if(transfer < 0.0 ) transfer = 0.0 ;
>> 724 // if(transfer < DBL_MIN ) transfer = DBL_MIN;
>> 725
>> 726 return transfer ;
>> 727 }
>> 728
>> 729 ///////////////////////////////////////////////////////////////////////
>> 730 //
>> 731 // Returns random PAI energy transfer according to passed
>> 732 // indexes of particle kinetic
>> 733
>> 734 G4double
>> 735 G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
>> 736 {
>> 737 G4int iTransferMax;
>> 738 G4double x1, x2, y1, y2, energyTransfer;
>> 739
>> 740 if(iTransfer == 0)
>> 741 {
>> 742 energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
>> 743 }
>> 744 else
>> 745 {
>> 746 iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
>> 747
>> 748 if ( iTransfer >= iTransferMax ) iTransfer = iTransferMax - 1;
>> 749
>> 750 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
>> 751 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
>> 752
>> 753 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
>> 754 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
>> 755
>> 756 if ( x1 == x2 ) energyTransfer = x2;
>> 757 else
>> 758 {
>> 759 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
>> 760 else
>> 761 {
>> 762 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
>> 763 }
>> 764 }
>> 765 }
>> 766 return energyTransfer;
>> 767 }
>> 768
>> 769 ///////////////////////////////////////////////////////////////////////
>> 770
>> 771 G4double G4PAIModel::SampleFluctuations( const G4Material* material,
>> 772 const G4DynamicParticle* aParticle,
>> 773 G4double&,
>> 774 G4double& step,
>> 775 G4double&)
>> 776 {
>> 777 size_t jMat = 0;
>> 778 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
>> 779 {
>> 780 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
>> 781 }
>> 782 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
>> 783
>> 784 fPAItransferTable = fPAIxscBank[jMat];
>> 785 fdNdxCutVector = fdNdxCutTable[jMat];
>> 786
>> 787 G4int iTkin, iTransfer, iPlace;
>> 788 G4long numOfCollisions=0;
>> 789
>> 790 //G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
>> 791 //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl;
>> 792 G4double loss = 0.0, charge2 ;
>> 793 G4double stepSum = 0., stepDelta, lambda, omega;
>> 794 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
>> 795 G4bool numb = true;
>> 796 G4double Tkin = aParticle->GetKineticEnergy() ;
>> 797 G4double MassRatio = fMass/aParticle->GetDefinition()->GetPDGMass() ;
>> 798 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
>> 799 charge2 = charge*charge ;
>> 800 G4double TkinScaled = Tkin*MassRatio ;
>> 801
>> 802 for(iTkin=0;iTkin<=fTotBin;iTkin++)
>> 803 {
>> 804 if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 805 }
>> 806 iPlace = iTkin - 1 ;
>> 807 if(iPlace < 0) iPlace = 0;
>> 808 else if(iPlace >= fTotBin) iPlace = fTotBin - 1;
>> 809 //G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
>> 810 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
>> 811 //G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
>> 812
>> 813 if(iTkin == fTotBin) // Fermi plato, try from left
>> 814 {
>> 815 meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2;
>> 816 if(meanNumber < 0.) meanNumber = 0. ;
>> 817 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
>> 818 // numOfCollisions = G4Poisson(meanNumber) ;
>> 819 if( meanNumber > 0.) lambda = step/meanNumber;
>> 820 else lambda = DBL_MAX;
>> 821 while(numb)
>> 822 {
>> 823 stepDelta = CLHEP::RandExponential::shoot(lambda);
>> 824 stepSum += stepDelta;
>> 825 if(stepSum >= step) break;
>> 826 numOfCollisions++;
>> 827 }
>> 828 //G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
>> 829
>> 830 while(numOfCollisions)
>> 831 {
>> 832 position = dNdxCut1+
>> 833 ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
>> 834
>> 835 for( iTransfer = 0;
>> 836 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
>> 837 {
>> 838 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
>> 839 }
>> 840 omega = GetEnergyTransfer(iPlace,position,iTransfer);
>> 841 //G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
>> 842 loss += omega;
>> 843 numOfCollisions-- ;
>> 844 }
>> 845 }
>> 846 else
>> 847 {
>> 848 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
>> 849 //G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl;
>> 850
>> 851 if(iTkin == 0) // Tkin is too small, trying from right only
>> 852 {
>> 853 meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2;
>> 854 if( meanNumber < 0. ) meanNumber = 0. ;
>> 855 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
>> 856 // numOfCollisions = G4Poisson(meanNumber) ;
>> 857 if( meanNumber > 0.) lambda = step/meanNumber;
>> 858 else lambda = DBL_MAX;
>> 859 while(numb)
>> 860 {
>> 861 stepDelta = CLHEP::RandExponential::shoot(lambda);
>> 862 stepSum += stepDelta;
>> 863 if(stepSum >= step) break;
>> 864 numOfCollisions++;
>> 865 }
>> 866
>> 867 //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
>> 868
>> 869 while(numOfCollisions)
>> 870 {
>> 871 position = dNdxCut2+
>> 872 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
>> 873
>> 874 for( iTransfer = 0;
>> 875 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
>> 876 {
>> 877 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
>> 878 }
>> 879 omega = GetEnergyTransfer(iPlace,position,iTransfer);
>> 880 //G4cout<<omega/keV<<"\t";
>> 881 loss += omega;
>> 882 numOfCollisions-- ;
>> 883 }
>> 884 }
>> 885 else // general case: Tkin between two vectors of the material
>> 886 {
>> 887 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
>> 888 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
>> 889 W = 1.0/(E2 - E1) ;
>> 890 W1 = (E2 - TkinScaled)*W ;
>> 891 W2 = (TkinScaled - E1)*W ;
>> 892
>> 893 //G4cout << fPAItransferTable->size() << G4endl;
>> 894 //G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
>> 895 // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
>> 896 //G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
>> 897 // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
>> 898
>> 899 meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 +
>> 900 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2;
>> 901 if(meanNumber<0.0) meanNumber = 0.0;
>> 902 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
>> 903 // numOfCollisions = G4Poisson(meanNumber) ;
>> 904 if( meanNumber > 0.) lambda = step/meanNumber;
>> 905 else lambda = DBL_MAX;
>> 906 while(numb)
>> 907 {
>> 908 stepDelta = CLHEP::RandExponential::shoot(lambda);
>> 909 stepSum += stepDelta;
>> 910 if(stepSum >= step) break;
>> 911 numOfCollisions++;
>> 912 }
>> 913
>> 914 //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
>> 915
>> 916 while(numOfCollisions)
>> 917 {
>> 918 position = dNdxCut1*W1 + dNdxCut2*W2 +
>> 919 ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 +
>> 920 dNdxCut2+
>> 921 ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand();
>> 922
>> 923 // G4cout<<position<<"\t" ;
>> 924
>> 925 for( iTransfer = 0;
>> 926 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
>> 927 {
>> 928 if( position >=
>> 929 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
>> 930 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
>> 931 {
>> 932 break ;
>> 933 }
>> 934 }
>> 935 omega = GetEnergyTransfer(iPlace,position,iTransfer);
>> 936 // G4cout<<omega/keV<<"\t";
>> 937 loss += omega;
>> 938 numOfCollisions-- ;
>> 939 }
>> 940 }
>> 941 }
>> 942 //G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
>> 943 // <<step/mm<<" mm"<<G4endl ;
>> 944 if(loss > Tkin) loss=Tkin;
>> 945 if(loss < 0. ) loss = 0.;
>> 946 return loss ;
331 947
332 } 948 }
333 949
334 ////////////////////////////////////////////// 950 //////////////////////////////////////////////////////////////////////
335 // 951 //
336 // Returns the statistical estimation of the e 952 // Returns the statistical estimation of the energy loss distribution variance
337 // 953 //
338 954
339 955
340 G4double G4PAIModel::Dispersion( const G4Mater 956 G4double G4PAIModel::Dispersion( const G4Material* material,
341 const G4Dynam 957 const G4DynamicParticle* aParticle,
342 const G4double tcut, << 958 G4double& tmax,
343 const G4double tmax, << 959 G4double& step )
344 const G4double step ) << 960 {
345 { << 961 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
346 G4double particleMass = aParticle->GetMass( << 962 for(G4int i = 0 ; i < fMeanNumber; i++)
347 G4double electronDensity = material->GetElec << 963 {
348 G4double kineticEnergy = aParticle->GetKinet << 964 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
349 G4double q = aParticle->GetCharge()/eplus; << 965 sumLoss += loss;
350 G4double etot = kineticEnergy + particleMass << 966 sumLoss2 += loss*loss;
351 G4double beta2 = kineticEnergy*(kineticEnerg << 967 }
352 G4double siga = (tmax/beta2 - 0.5*tcut) * t << 968 meanLoss = sumLoss/fMeanNumber;
353 * electronDensity * q * q; << 969 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
354 << 970 return sigma2;
355 return siga; <<
356 } 971 }
357 972
358 ////////////////////////////////////////////// 973 /////////////////////////////////////////////////////////////////////
359 974
360 G4double G4PAIModel::MaxSecondaryEnergy( const 975 G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
361 G4double kinEnergy) 976 G4double kinEnergy)
362 { 977 {
363 SetParticle(p); <<
364 G4double tmax = kinEnergy; 978 G4double tmax = kinEnergy;
365 if(p == fElectron) { tmax *= 0.5; } << 979 if(p == fElectron) tmax *= 0.5;
366 else if(p != fPositron) { 980 else if(p != fPositron) {
367 G4double ratio= electron_mass_c2/fMass; << 981 G4double mass = p->GetPDGMass();
368 G4double gamma= kinEnergy/fMass + 1.0; << 982 G4double ratio= electron_mass_c2/mass;
>> 983 G4double gamma= kinEnergy/mass + 1.0;
369 tmax = 2.0*electron_mass_c2*(gamma*gamma - 984 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
370 (1. + 2.0*gamma*ratio + rati 985 (1. + 2.0*gamma*ratio + ratio*ratio);
371 } 986 }
372 return tmax; 987 return tmax;
373 } 988 }
374 989
375 ////////////////////////////////////////////// 990 ///////////////////////////////////////////////////////////////
376 991
377 void G4PAIModel::DefineForRegion(const G4Regio 992 void G4PAIModel::DefineForRegion(const G4Region* r)
378 { 993 {
379 fPAIRegionVector.push_back(r); 994 fPAIRegionVector.push_back(r);
380 } 995 }
381 996
382 ////////////////////////////////////////////// << 997 //
>> 998 //
>> 999 /////////////////////////////////////////////////
>> 1000
>> 1001
>> 1002
>> 1003
>> 1004
383 1005
384 1006