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