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

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 32.2 KB
Line 
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// $Id: G4PAIModel.cc,v 1.52 2010/06/03 07:28:39 grichine Exp $
27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class
32// File name: G4PAIModel.cc
33//
34// Author: Vladimir.Grichine@cern.ch on base of Vladimir Ivanchenko model interface
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// 26.07.09 Fixed logic to work with several materials (V.Ivantchenko)
44//
45
46#include "G4Region.hh"
47#include "G4PhysicsLogVector.hh"
48#include "G4PhysicsFreeVector.hh"
49#include "G4PhysicsTable.hh"
50#include "G4ProductionCutsTable.hh"
51#include "G4MaterialCutsCouple.hh"
52#include "G4MaterialTable.hh"
53#include "G4SandiaTable.hh"
54#include "G4OrderedTable.hh"
55
56#include "G4PAIModel.hh"
57#include "Randomize.hh"
58#include "G4Electron.hh"
59#include "G4Positron.hh"
60#include "G4Poisson.hh"
61#include "G4Step.hh"
62#include "G4Material.hh"
63#include "G4DynamicParticle.hh"
64#include "G4ParticleDefinition.hh"
65#include "G4ParticleChangeForLoss.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 G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
151 if(isInitialised) return;
152 isInitialised = true;
153
154 SetParticle(p);
155 fLowestKineticEnergy = fMass*(fLowestGamma - 1.0);
156 fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
157
158 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
159 fHighestKineticEnergy,
160 fTotBin);
161
162 fParticleChange = GetParticleChangeForLoss();
163
164 // Prepare initialization
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 //G4cout << "Reg <" <<curReg->GetName() << "> mat <"
180 // << fMaterial->GetName() << "> fCouple= "
181 // << fCutCouple<<" " << p->GetParticleName() <<G4endl;
182 if( fCutCouple ) {
183 fMaterialCutsCoupleVector.push_back(fCutCouple);
184
185 fPAItransferTable = new G4PhysicsTable(fTotBin+1);
186 fPAIdEdxTable = new G4PhysicsTable(fTotBin+1);
187
188 fDeltaCutInKinEnergy =
189 (*theCoupleTable->GetEnergyCutsVector(1))[fCutCouple->GetIndex()];
190
191 //ComputeSandiaPhotoAbsCof();
192 BuildPAIonisationTable();
193
194 fPAIxscBank.push_back(fPAItransferTable);
195 fPAIdEdxBank.push_back(fPAIdEdxTable);
196 fdEdxTable.push_back(fdEdxVector);
197
198 BuildLambdaVector();
199 fdNdxCutTable.push_back(fdNdxCutVector);
200 fLambdaTable.push_back(fLambdaVector);
201 }
202 }
203 }
204}
205
206//////////////////////////////////////////////////////////////////
207
208void G4PAIModel::InitialiseMe(const G4ParticleDefinition*)
209{}
210
211//////////////////////////////////////////////////////////////////
212
213void G4PAIModel::ComputeSandiaPhotoAbsCof()
214{
215 G4int i, j, numberOfElements ;
216 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
217
218 G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
219 numberOfElements = (*theMaterialTable)[fMatIndex]->
220 GetNumberOfElements();
221 G4int* thisMaterialZ = new G4int[numberOfElements] ;
222
223 for(i=0;i<numberOfElements;i++)
224 {
225 thisMaterialZ[i] =
226 (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
227 }
228 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
229 (thisMaterialZ,numberOfElements) ;
230
231 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
232 ( thisMaterialZ ,
233 (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
234 numberOfElements,fSandiaIntervalNumber) ;
235
236 fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
237
238 for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ;
239
240 for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
241 {
242 fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ;
243
244 for( j = 1; j < 5 ; j++ )
245 {
246 fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
247 GetPhotoAbsorpCof(i+1,j)*
248 (*theMaterialTable)[fMatIndex]->GetDensity() ;
249 }
250 }
251 // delete[] thisMaterialZ ;
252}
253
254////////////////////////////////////////////////////////////////////////////
255//
256// Build tables for the ionization energy loss
257// the tables are built for MATERIALS
258// *********
259
260void G4PAIModel::BuildPAIonisationTable()
261{
262 G4double LowEdgeEnergy , ionloss ;
263 G4double tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ;
264
265 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
266 fHighestKineticEnergy,
267 fTotBin);
268 G4SandiaTable* sandia = fMaterial->GetSandiaTable();
269
270 Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV;
271 deltaLow = 100.*eV; // 0.5*eV ;
272
273 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy
274 {
275 LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
276 tau = LowEdgeEnergy/fMass ;
277 gamma = tau +1. ;
278 // G4cout<<"gamma = "<<gamma<<endl ;
279 bg2 = tau*( tau + 2. );
280 Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
281 // Tmax = std::min(fDeltaCutInKinEnergy, Tmax);
282 Tkin = Tmax ;
283
284 if ( Tmax < Tmin + deltaLow ) // low energy safety
285 Tkin = Tmin + deltaLow ;
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}
317
318///////////////////////////////////////////////////////////////////////
319//
320// Build mean free path tables for the delta ray production process
321// tables are built for MATERIALS
322//
323
324void G4PAIModel::BuildLambdaVector()
325{
326 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
327 fHighestKineticEnergy,
328 fTotBin ) ;
329 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
330 fHighestKineticEnergy,
331 fTotBin ) ;
332 if(fVerbose > 1)
333 {
334 G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
335 <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
336 }
337 for (G4int i = 0 ; i <= fTotBin ; i++ )
338 {
339 G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
340 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
341 if(dNdxCut < 0.0) dNdxCut = 0.0;
342 // G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
343 G4double lambda = DBL_MAX;
344 if(dNdxCut > 0.0) lambda = 1.0/dNdxCut;
345
346 fLambdaVector->PutValue(i, lambda) ;
347 fdNdxCutVector->PutValue(i, dNdxCut) ;
348 }
349}
350
351///////////////////////////////////////////////////////////////////////
352//
353// Returns integral PAI cross section for energy transfers >= transferCut
354
355G4double
356G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
357{
358 G4int iTransfer;
359 G4double x1, x2, y1, y2, dNdxCut;
360 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
361 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
362 // <<G4endl;
363 for( iTransfer = 0 ;
364 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
365 iTransfer++)
366 {
367 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
368 {
369 break ;
370 }
371 }
372 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
373 {
374 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
375 }
376 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
377 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
378 if(y1 < 0.0) y1 = 0.0;
379 if(y2 < 0.0) y2 = 0.0;
380 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
381 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
382 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
383 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
384
385 if ( y1 == y2 ) dNdxCut = y2 ;
386 else
387 {
388 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
389 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
390 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
391 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
392 }
393 // G4cout<<""<<dNdxCut<<G4endl;
394 if(dNdxCut < 0.0) dNdxCut = 0.0;
395 return dNdxCut ;
396}
397///////////////////////////////////////////////////////////////////////
398//
399// Returns integral dEdx for energy transfers >= transferCut
400
401G4double
402G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
403{
404 G4int iTransfer;
405 G4double x1, x2, y1, y2, dEdxCut;
406 //G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
407 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
408 // <<G4endl;
409 for( iTransfer = 0 ;
410 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
411 iTransfer++)
412 {
413 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
414 {
415 break ;
416 }
417 }
418 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
419 {
420 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
421 }
422 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
423 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
424 if(y1 < 0.0) y1 = 0.0;
425 if(y2 < 0.0) y2 = 0.0;
426 //G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
427 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
428 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
429 //G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
430
431 if ( y1 == y2 ) dEdxCut = y2 ;
432 else
433 {
434 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
435 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
436 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
437 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
438 }
439 //G4cout<<""<<dEdxCut<<G4endl;
440 if(dEdxCut < 0.0) dEdxCut = 0.0;
441 return dEdxCut ;
442}
443
444//////////////////////////////////////////////////////////////////////////////
445
446G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
447 const G4ParticleDefinition* p,
448 G4double kineticEnergy,
449 G4double cutEnergy)
450{
451 G4int iTkin,iPlace;
452
453 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
454 G4double cut = cutEnergy;
455
456 G4double massRatio = fMass/p->GetPDGMass();
457 G4double scaledTkin = kineticEnergy*massRatio;
458 G4double charge = p->GetPDGCharge();
459 G4double charge2 = charge*charge;
460 const G4MaterialCutsCouple* matCC = CurrentCouple();
461
462 size_t jMat = 0;
463 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
464 {
465 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
466 }
467 //G4cout << "jMat= " << jMat
468 // << " jMax= " << fMaterialCutsCoupleVector.size()
469 // << " matCC: " << matCC;
470 //if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
471 // G4cout << G4endl;
472 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
473
474 fPAIdEdxTable = fPAIdEdxBank[jMat];
475 fdEdxVector = fdEdxTable[jMat];
476 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
477 {
478 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
479 }
480 //G4cout << "E= " << scaledTkin << " iTkin= " << iTkin
481 // << " " << fdEdxVector->GetVectorLength()<<G4endl;
482 iPlace = iTkin - 1;
483 if(iPlace < 0) iPlace = 0;
484 else if(iPlace >= fTotBin) iPlace = fTotBin-1;
485 G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
486 if( dEdx < 0.) dEdx = 0.;
487 return dEdx;
488}
489
490/////////////////////////////////////////////////////////////////////////
491
492G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
493 const G4ParticleDefinition* p,
494 G4double kineticEnergy,
495 G4double cutEnergy,
496 G4double maxEnergy )
497{
498 G4int iTkin,iPlace;
499 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
500 if(tmax <= cutEnergy) return 0.0;
501 G4double massRatio = fMass/p->GetPDGMass();
502 G4double scaledTkin = kineticEnergy*massRatio;
503 G4double charge = p->GetPDGCharge();
504 G4double charge2 = charge*charge, cross, cross1, cross2;
505 const G4MaterialCutsCouple* matCC = CurrentCouple();
506
507 size_t jMat = 0;
508 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
509 {
510 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
511 }
512 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
513
514 fPAItransferTable = fPAIxscBank[jMat];
515
516 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
517 {
518 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
519 }
520 iPlace = iTkin - 1;
521 if(iPlace < 0) iPlace = 0;
522 else if(iPlace >= fTotBin) iPlace = fTotBin-1;
523
524 //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
525 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
526 cross1 = GetdNdxCut(iPlace,tmax) ;
527 //G4cout<<"cross1 = "<<cross1<<G4endl;
528 cross2 = GetdNdxCut(iPlace,cutEnergy) ;
529 //G4cout<<"cross2 = "<<cross2<<G4endl;
530 cross = (cross2-cross1)*charge2;
531 //G4cout<<"cross = "<<cross<<G4endl;
532 if( cross < DBL_MIN) cross = 0.0;
533 // if( cross2 < DBL_MIN) cross2 = DBL_MIN;
534
535 // return cross2;
536 return cross;
537}
538
539///////////////////////////////////////////////////////////////////////////
540//
541// It is analog of PostStepDoIt in terms of secondary electron.
542//
543
544void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
545 const G4MaterialCutsCouple* matCC,
546 const G4DynamicParticle* dp,
547 G4double tmin,
548 G4double maxEnergy)
549{
550 size_t jMat = 0;
551 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
552 {
553 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
554 }
555 if(jMat == fMaterialCutsCoupleVector.size()) return;
556
557 fPAItransferTable = fPAIxscBank[jMat];
558 fdNdxCutVector = fdNdxCutTable[jMat];
559
560 G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy);
561 if( tmin >= tmax && fVerbose > 0)
562 {
563 G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
564 }
565 G4ThreeVector direction= dp->GetMomentumDirection();
566 G4double particleMass = dp->GetMass();
567 G4double kineticEnergy = dp->GetKineticEnergy();
568
569 G4double massRatio = fMass/particleMass;
570 G4double scaledTkin = kineticEnergy*massRatio;
571 G4double totalEnergy = kineticEnergy + particleMass;
572 G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
573
574 G4double deltaTkin = GetPostStepTransfer(scaledTkin);
575
576 // G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl;
577
578 if( deltaTkin <= 0. && fVerbose > 0)
579 {
580 G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
581 }
582 if( deltaTkin <= 0.) return;
583
584 if( deltaTkin > tmax) deltaTkin = tmax;
585
586 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
587 G4double totalMomentum = sqrt(pSquare);
588 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
589 /(deltaTotalMomentum * totalMomentum);
590
591 if( costheta > 0.99999 ) costheta = 0.99999;
592 G4double sintheta = 0.0;
593 G4double sin2 = 1. - costheta*costheta;
594 if( sin2 > 0.) sintheta = sqrt(sin2);
595
596 // direction of the delta electron
597 G4double phi = twopi*G4UniformRand();
598 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
599
600 G4ThreeVector deltaDirection(dirx,diry,dirz);
601 deltaDirection.rotateUz(direction);
602 deltaDirection.unit();
603
604 // primary change
605 kineticEnergy -= deltaTkin;
606 G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
607 direction = dir.unit();
608 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
609 fParticleChange->SetProposedMomentumDirection(direction);
610
611 // create G4DynamicParticle object for e- delta ray
612 G4DynamicParticle* deltaRay = new G4DynamicParticle;
613 deltaRay->SetDefinition(G4Electron::Electron());
614 deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ???
615 deltaRay->SetMomentumDirection(deltaDirection);
616
617 vdp->push_back(deltaRay);
618}
619
620
621///////////////////////////////////////////////////////////////////////
622//
623// Returns post step PAI energy transfer > cut electron energy according to passed
624// scaled kinetic energy of particle
625
626G4double
627G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
628{
629 // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
630
631 G4int iTkin, iTransfer, iPlace ;
632 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
633
634 for(iTkin=0;iTkin<=fTotBin;iTkin++)
635 {
636 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
637 }
638 iPlace = iTkin - 1 ;
639 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
640 if(iPlace < 0) iPlace = 0;
641 else if(iPlace >= fTotBin) iPlace = fTotBin-1;
642 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
643 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
644
645
646 if(iTkin == fTotBin) // Fermi plato, try from left
647 {
648 position = dNdxCut1*G4UniformRand() ;
649
650 for( iTransfer = 0;
651 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
652 {
653 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
654 }
655 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
656 }
657 else
658 {
659 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
660 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
661 if(iTkin == 0) // Tkin is too small, trying from right only
662 {
663 position = dNdxCut2*G4UniformRand() ;
664
665 for( iTransfer = 0;
666 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
667 {
668 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
669 }
670 transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
671 }
672 else // general case: Tkin between two vectors of the material
673 {
674 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
675 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
676 W = 1.0/(E2 - E1) ;
677 W1 = (E2 - scaledTkin)*W ;
678 W2 = (scaledTkin - E1)*W ;
679
680 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
681
682 // G4cout<<position<<"\t" ;
683
684 G4int iTrMax1, iTrMax2, iTrMax;
685
686 iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
687 iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength());
688
689 if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
690 else iTrMax = iTrMax1;
691
692
693 for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
694 {
695 if( position >=
696 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
697 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
698 }
699 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
700 }
701 }
702 // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
703 if(transfer < 0.0 ) transfer = 0.0 ;
704 // if(transfer < DBL_MIN ) transfer = DBL_MIN;
705
706 return transfer ;
707}
708
709///////////////////////////////////////////////////////////////////////
710//
711// Returns random PAI energy transfer according to passed
712// indexes of particle kinetic
713
714G4double
715G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
716{
717 G4int iTransferMax;
718 G4double x1, x2, y1, y2, energyTransfer;
719
720 if(iTransfer == 0)
721 {
722 energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
723 }
724 else
725 {
726 iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
727
728 if ( iTransfer >= iTransferMax ) iTransfer = iTransferMax - 1;
729
730 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
731 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
732
733 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
734 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
735
736 if ( x1 == x2 ) energyTransfer = x2;
737 else
738 {
739 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
740 else
741 {
742 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
743 }
744 }
745 }
746 return energyTransfer;
747}
748
749///////////////////////////////////////////////////////////////////////
750
751G4double G4PAIModel::SampleFluctuations( const G4Material* material,
752 const G4DynamicParticle* aParticle,
753 G4double&,
754 G4double& step,
755 G4double&)
756{
757 size_t jMat = 0;
758 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
759 {
760 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
761 }
762 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
763
764 fPAItransferTable = fPAIxscBank[jMat];
765 fdNdxCutVector = fdNdxCutTable[jMat];
766
767 G4int iTkin, iTransfer, iPlace;
768 G4long numOfCollisions=0;
769
770 //G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
771 //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl;
772 G4double loss = 0.0, charge2 ;
773 G4double stepSum = 0., stepDelta, lambda, omega;
774 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
775 G4bool numb = true;
776 G4double Tkin = aParticle->GetKineticEnergy() ;
777 G4double MassRatio = fMass/aParticle->GetDefinition()->GetPDGMass() ;
778 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
779 charge2 = charge*charge ;
780 G4double TkinScaled = Tkin*MassRatio ;
781
782 for(iTkin=0;iTkin<=fTotBin;iTkin++)
783 {
784 if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
785 }
786 iPlace = iTkin - 1 ;
787 if(iPlace < 0) iPlace = 0;
788 else if(iPlace >= fTotBin) iPlace = fTotBin - 1;
789 //G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
790 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
791 //G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
792
793 if(iTkin == fTotBin) // Fermi plato, try from left
794 {
795 meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2;
796 if(meanNumber < 0.) meanNumber = 0. ;
797 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
798 // numOfCollisions = G4Poisson(meanNumber) ;
799 if( meanNumber > 0.) lambda = step/meanNumber;
800 else lambda = DBL_MAX;
801 while(numb)
802 {
803 stepDelta = CLHEP::RandExponential::shoot(lambda);
804 stepSum += stepDelta;
805 if(stepSum >= step) break;
806 numOfCollisions++;
807 }
808 //G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
809
810 while(numOfCollisions)
811 {
812 position = dNdxCut1+
813 ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
814
815 for( iTransfer = 0;
816 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
817 {
818 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
819 }
820 omega = GetEnergyTransfer(iPlace,position,iTransfer);
821 //G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
822 loss += omega;
823 numOfCollisions-- ;
824 }
825 }
826 else
827 {
828 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
829 //G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl;
830
831 if(iTkin == 0) // Tkin is too small, trying from right only
832 {
833 meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2;
834 if( meanNumber < 0. ) meanNumber = 0. ;
835 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
836 // numOfCollisions = G4Poisson(meanNumber) ;
837 if( meanNumber > 0.) lambda = step/meanNumber;
838 else lambda = DBL_MAX;
839 while(numb)
840 {
841 stepDelta = CLHEP::RandExponential::shoot(lambda);
842 stepSum += stepDelta;
843 if(stepSum >= step) break;
844 numOfCollisions++;
845 }
846
847 //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
848
849 while(numOfCollisions)
850 {
851 position = dNdxCut2+
852 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
853
854 for( iTransfer = 0;
855 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
856 {
857 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
858 }
859 omega = GetEnergyTransfer(iPlace,position,iTransfer);
860 //G4cout<<omega/keV<<"\t";
861 loss += omega;
862 numOfCollisions-- ;
863 }
864 }
865 else // general case: Tkin between two vectors of the material
866 {
867 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
868 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
869 W = 1.0/(E2 - E1) ;
870 W1 = (E2 - TkinScaled)*W ;
871 W2 = (TkinScaled - E1)*W ;
872
873 //G4cout << fPAItransferTable->size() << G4endl;
874 //G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
875 // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
876 //G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
877 // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
878
879 meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 +
880 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2;
881 if(meanNumber<0.0) meanNumber = 0.0;
882 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
883 // numOfCollisions = G4Poisson(meanNumber) ;
884 if( meanNumber > 0.) lambda = step/meanNumber;
885 else lambda = DBL_MAX;
886 while(numb)
887 {
888 stepDelta = CLHEP::RandExponential::shoot(lambda);
889 stepSum += stepDelta;
890 if(stepSum >= step) break;
891 numOfCollisions++;
892 }
893
894 //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
895
896 while(numOfCollisions)
897 {
898 position = dNdxCut1*W1 + dNdxCut2*W2 +
899 ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 +
900 dNdxCut2+
901 ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand();
902
903 // G4cout<<position<<"\t" ;
904
905 for( iTransfer = 0;
906 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
907 {
908 if( position >=
909 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
910 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
911 {
912 break ;
913 }
914 }
915 omega = GetEnergyTransfer(iPlace,position,iTransfer);
916 // G4cout<<omega/keV<<"\t";
917 loss += omega;
918 numOfCollisions-- ;
919 }
920 }
921 }
922 //G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
923 // <<step/mm<<" mm"<<G4endl ;
924 if(loss > Tkin) loss=Tkin;
925 if(loss < 0. ) loss = 0.;
926 return loss ;
927
928}
929
930//////////////////////////////////////////////////////////////////////
931//
932// Returns the statistical estimation of the energy loss distribution variance
933//
934
935
936G4double G4PAIModel::Dispersion( const G4Material* material,
937 const G4DynamicParticle* aParticle,
938 G4double& tmax,
939 G4double& step )
940{
941 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
942 for(G4int i = 0 ; i < fMeanNumber; i++)
943 {
944 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
945 sumLoss += loss;
946 sumLoss2 += loss*loss;
947 }
948 meanLoss = sumLoss/fMeanNumber;
949 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
950 return sigma2;
951}
952
953/////////////////////////////////////////////////////////////////////
954
955G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
956 G4double kinEnergy)
957{
958 G4double tmax = kinEnergy;
959 if(p == fElectron) tmax *= 0.5;
960 else if(p != fPositron) {
961 G4double mass = p->GetPDGMass();
962 G4double ratio= electron_mass_c2/mass;
963 G4double gamma= kinEnergy/mass + 1.0;
964 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
965 (1. + 2.0*gamma*ratio + ratio*ratio);
966 }
967 return tmax;
968}
969
970///////////////////////////////////////////////////////////////
971
972void G4PAIModel::DefineForRegion(const G4Region* r)
973{
974 fPAIRegionVector.push_back(r);
975}
976
977//
978//
979/////////////////////////////////////////////////
980
981
982
983
984
985
Note: See TracBrowser for help on using the repository browser.