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

Last change on this file since 980 was 961, checked in by garnier, 17 years ago

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