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