source: trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

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