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

Last change on this file since 1123 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 32.1 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.47 2009/04/09 18:41:18 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-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 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 fParticleChange = GetParticleChangeForLoss();
162
163 // Prepare initialization
164
165 fPAItransferTable = new G4PhysicsTable(fTotBin);
166 fPAIdEdxTable = new G4PhysicsTable(fTotBin);
167
168 const G4ProductionCutsTable* theCoupleTable =
169 G4ProductionCutsTable::GetProductionCutsTable();
170 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
171 size_t numOfMat = G4Material::GetNumberOfMaterials();
172 size_t numRegions = fPAIRegionVector.size();
173
174 for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
175 {
176 const G4Region* curReg = fPAIRegionVector[iReg];
177 for(size_t jMat = 0; jMat < numOfMat; ++jMat) // region material loop
178 {
179 fMaterial = (*theMaterialTable)[jMat];
180 fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial,
181 curReg->GetProductionCuts() );
182 //G4cout << "Reg <" <<curReg->GetName() << "> mat <"
183 // << fMaterial->GetName() << "> fCouple= "
184 // << fCutCouple<<G4endl;
185 if( fCutCouple ) {
186 fMaterialCutsCoupleVector.push_back(fCutCouple);
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 if(fdEdxVector) delete fdEdxVector;
266 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
267 fHighestKineticEnergy,
268 fTotBin);
269 G4SandiaTable* sandia = fMaterial->GetSandiaTable();
270
271 Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV;
272 deltaLow = 100.*eV; // 0.5*eV ;
273
274 for (G4int i = 0 ; i < fTotBin ; i++) //The loop for the kinetic energy
275 {
276 LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
277 tau = LowEdgeEnergy/fMass ;
278 gamma = tau +1. ;
279 // G4cout<<"gamma = "<<gamma<<endl ;
280 bg2 = tau*( tau + 2. );
281 Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
282 // Tmax = std::min(fDeltaCutInKinEnergy, Tmax);
283 Tkin = Tmax ;
284
285 // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
286 // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
287
288 if ( Tmax < Tmin + deltaLow ) // low energy safety
289 Tkin = Tmin + deltaLow ;
290
291 /*
292 G4PAIxSection protonPAI( fMatIndex,
293 Tkin,
294 bg2,
295 fSandiaPhotoAbsCof,
296 fSandiaIntervalNumber ) ;
297 */
298 fPAIySection.Initialize(fMaterial, Tkin, bg2);
299
300 // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
301 // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
302 // G4cout<<"protonPAI.GetSplineSize() = "<<
303 // protonPAI.GetSplineSize()<<G4endl<<G4endl ;
304
305 G4int n = fPAIySection.GetSplineSize();
306 G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n) ;
307 G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
308
309 for( G4int k = 0 ; k < n; k++ )
310 {
311 transferVector->PutValue( k ,
312 fPAIySection.GetSplineEnergy(k+1),
313 fPAIySection.GetIntegralPAIySection(k+1) ) ;
314 dEdxVector->PutValue( k ,
315 fPAIySection.GetSplineEnergy(k+1),
316 fPAIySection.GetIntegralPAIdEdx(k+1) ) ;
317 }
318 ionloss = fPAIySection.GetMeanEnergyLoss() ; // total <dE/dx>
319
320 if ( ionloss < DBL_MIN) ionloss = DBL_MIN;
321 fdEdxVector->PutValue(i,ionloss) ;
322
323 fPAItransferTable->insertAt(i,transferVector) ;
324 fPAIdEdxTable->insertAt(i,dEdxVector) ;
325
326 } // end of Tkin loop
327 // theLossTable->insert(fdEdxVector);
328 // end of material loop
329 // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;
330 // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;
331}
332
333///////////////////////////////////////////////////////////////////////
334//
335// Build mean free path tables for the delta ray production process
336// tables are built for MATERIALS
337//
338
339void G4PAIModel::BuildLambdaVector()
340{
341 //G4double kCarTolerance = G4GeometryTolerance::GetInstance()
342 // ->GetSurfaceTolerance();
343
344 if (fLambdaVector) delete fLambdaVector;
345 if (fdNdxCutVector) delete fdNdxCutVector;
346
347 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
348 fHighestKineticEnergy,
349 fTotBin ) ;
350 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
351 fHighestKineticEnergy,
352 fTotBin ) ;
353 if(fVerbose > 1)
354 {
355 G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
356 <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
357 }
358 for (G4int i = 0 ; i < fTotBin ; i++ )
359 {
360 G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
361 G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
362
363 // if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance ; // Mmm ???
364
365 fLambdaVector->PutValue(i, lambda) ;
366 fdNdxCutVector->PutValue(i, dNdxCut) ;
367 }
368}
369
370///////////////////////////////////////////////////////////////////////
371//
372// Returns integral PAI cross section for energy transfers >= transferCut
373
374G4double
375G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
376{
377 G4int iTransfer;
378 G4double x1, x2, y1, y2, dNdxCut;
379 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
380 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
381 // <<G4endl;
382 for( iTransfer = 0 ;
383 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
384 iTransfer++)
385 {
386 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
387 {
388 break ;
389 }
390 }
391 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
392 {
393 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
394 }
395 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
396 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
397 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
398 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
399 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
400 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
401
402 if ( y1 == y2 ) dNdxCut = y2 ;
403 else
404 {
405 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
406 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
407 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
408 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
409 }
410 // G4cout<<""<<dNdxCut<<G4endl;
411 return dNdxCut ;
412}
413///////////////////////////////////////////////////////////////////////
414//
415// Returns integral dEdx for energy transfers >= transferCut
416
417G4double
418G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
419{
420 G4int iTransfer;
421 G4double x1, x2, y1, y2, dEdxCut;
422 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
423 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
424 // <<G4endl;
425 for( iTransfer = 0 ;
426 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
427 iTransfer++)
428 {
429 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
430 {
431 break ;
432 }
433 }
434 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
435 {
436 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
437 }
438 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
439 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
440 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
441 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
442 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
443 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
444
445 if ( y1 == y2 ) dEdxCut = y2 ;
446 else
447 {
448 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
449 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
450 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
451 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
452 }
453 // G4cout<<""<<dEdxCut<<G4endl;
454 return dEdxCut ;
455}
456
457//////////////////////////////////////////////////////////////////////////////
458
459G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
460 const G4ParticleDefinition* p,
461 G4double kineticEnergy,
462 G4double cutEnergy)
463{
464 G4int iTkin,iPlace;
465 size_t jMat;
466
467 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
468 G4double cut = cutEnergy;
469
470 G4double massRatio = fMass/p->GetPDGMass();
471 G4double scaledTkin = kineticEnergy*massRatio;
472 G4double charge = p->GetPDGCharge();
473 G4double charge2 = charge*charge;
474 const G4MaterialCutsCouple* matCC = CurrentCouple();
475
476 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
477 {
478 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
479 }
480 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
481
482 fPAIdEdxTable = fPAIdEdxBank[jMat];
483 fdEdxVector = fdEdxTable[jMat];
484 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
485 {
486 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
487 }
488 iPlace = iTkin - 1;
489 if(iPlace < 0) iPlace = 0;
490 G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
491 if( dEdx < 0.) dEdx = 0.;
492 return dEdx;
493}
494
495/////////////////////////////////////////////////////////////////////////
496
497G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
498 const G4ParticleDefinition* p,
499 G4double kineticEnergy,
500 G4double cutEnergy,
501 G4double maxEnergy )
502{
503 G4int iTkin,iPlace;
504 size_t jMat;
505 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
506 if(tmax <= cutEnergy) return 0.0;
507 G4double massRatio = fMass/p->GetPDGMass();
508 G4double scaledTkin = kineticEnergy*massRatio;
509 G4double charge = p->GetPDGCharge();
510 G4double charge2 = charge*charge, cross, cross1, cross2;
511 const G4MaterialCutsCouple* matCC = CurrentCouple();
512
513 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
514 {
515 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
516 }
517 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
518
519 fPAItransferTable = fPAIxscBank[jMat];
520
521 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
522 {
523 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
524 }
525 iPlace = iTkin - 1;
526 if(iPlace < 0) iPlace = 0;
527
528 // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
529 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
530 cross1 = GetdNdxCut(iPlace,tmax) ;
531 // G4cout<<"cross1 = "<<cross1<<G4endl;
532 cross2 = GetdNdxCut(iPlace,cutEnergy) ;
533 // G4cout<<"cross2 = "<<cross2<<G4endl;
534 cross = (cross2-cross1)*charge2;
535 // G4cout<<"cross = "<<cross<<G4endl;
536 if( cross < DBL_MIN) cross = DBL_MIN;
537 // if( cross2 < DBL_MIN) cross2 = DBL_MIN;
538
539 // return cross2;
540 return cross;
541}
542
543///////////////////////////////////////////////////////////////////////////
544//
545// It is analog of PostStepDoIt in terms of secondary electron.
546//
547
548void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
549 const G4MaterialCutsCouple* matCC,
550 const G4DynamicParticle* dp,
551 G4double tmin,
552 G4double maxEnergy)
553{
554 size_t jMat;
555 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
556 {
557 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
558 }
559 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
560
561 fPAItransferTable = fPAIxscBank[jMat];
562 fdNdxCutVector = fdNdxCutTable[jMat];
563
564 G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy);
565 if( tmin >= tmax && fVerbose > 0)
566 {
567 G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
568 }
569 G4ThreeVector direction= dp->GetMomentumDirection();
570 G4double particleMass = dp->GetMass();
571 G4double kineticEnergy = dp->GetKineticEnergy();
572
573 G4double massRatio = fMass/particleMass;
574 G4double scaledTkin = kineticEnergy*massRatio;
575 G4double totalEnergy = kineticEnergy + particleMass;
576 G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
577
578 G4double deltaTkin = GetPostStepTransfer(scaledTkin);
579
580 // G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl;
581
582 if( deltaTkin <= 0. && fVerbose > 0)
583 {
584 G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
585 }
586 if( deltaTkin <= 0.) return;
587
588 if( deltaTkin > tmax) deltaTkin = tmax;
589
590 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
591 G4double totalMomentum = sqrt(pSquare);
592 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
593 /(deltaTotalMomentum * totalMomentum);
594
595 if( costheta > 0.99999 ) costheta = 0.99999;
596 G4double sintheta = 0.0;
597 G4double sin2 = 1. - costheta*costheta;
598 if( sin2 > 0.) sintheta = sqrt(sin2);
599
600 // direction of the delta electron
601 G4double phi = twopi*G4UniformRand();
602 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
603
604 G4ThreeVector deltaDirection(dirx,diry,dirz);
605 deltaDirection.rotateUz(direction);
606 deltaDirection.unit();
607
608 // primary change
609 kineticEnergy -= deltaTkin;
610 G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
611 direction = dir.unit();
612 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
613 fParticleChange->SetProposedMomentumDirection(direction);
614
615 // create G4DynamicParticle object for e- delta ray
616 G4DynamicParticle* deltaRay = new G4DynamicParticle;
617 deltaRay->SetDefinition(G4Electron::Electron());
618 deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ???
619 deltaRay->SetMomentumDirection(deltaDirection);
620
621 vdp->push_back(deltaRay);
622}
623
624
625///////////////////////////////////////////////////////////////////////
626//
627// Returns post step PAI energy transfer > cut electron energy according to passed
628// scaled kinetic energy of particle
629
630G4double
631G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
632{
633 // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
634
635 G4int iTkin, iTransfer, iPlace ;
636 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
637
638 for(iTkin=0;iTkin<fTotBin;iTkin++)
639 {
640 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
641 }
642 iPlace = iTkin - 1 ;
643 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
644 if(iPlace < 0) iPlace = 0;
645 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
646 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
647
648
649 if(iTkin == fTotBin) // Fermi plato, try from left
650 {
651 position = dNdxCut1*G4UniformRand() ;
652
653 for( iTransfer = 0;
654 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
655 {
656 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
657 }
658 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
659 }
660 else
661 {
662 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
663 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
664 if(iTkin == 0) // Tkin is too small, trying from right only
665 {
666 position = dNdxCut2*G4UniformRand() ;
667
668 for( iTransfer = 0;
669 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
670 {
671 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
672 }
673 transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
674 }
675 else // general case: Tkin between two vectors of the material
676 {
677 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
678 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
679 W = 1.0/(E2 - E1) ;
680 W1 = (E2 - scaledTkin)*W ;
681 W2 = (scaledTkin - E1)*W ;
682
683 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
684
685 // G4cout<<position<<"\t" ;
686
687 G4int iTrMax1, iTrMax2, iTrMax;
688
689 iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
690 iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength());
691
692 if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
693 else iTrMax = iTrMax1;
694
695
696 for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
697 {
698 if( position >=
699 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
700 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
701 }
702 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
703 }
704 }
705 // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
706 if(transfer < 0.0 ) transfer = 0.0 ;
707 // if(transfer < DBL_MIN ) transfer = DBL_MIN;
708
709 return transfer ;
710}
711
712///////////////////////////////////////////////////////////////////////
713//
714// Returns random PAI energy transfer according to passed
715// indexes of particle kinetic
716
717G4double
718G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
719{
720 G4int iTransferMax;
721 G4double x1, x2, y1, y2, energyTransfer;
722
723 if(iTransfer == 0)
724 {
725 energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
726 }
727 else
728 {
729 iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
730
731 if ( iTransfer >= iTransferMax ) iTransfer = iTransferMax - 1;
732
733 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
734 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
735
736 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
737 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
738
739 if ( x1 == x2 ) energyTransfer = x2;
740 else
741 {
742 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
743 else
744 {
745 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
746 }
747 }
748 }
749 return energyTransfer;
750}
751
752///////////////////////////////////////////////////////////////////////
753
754G4double G4PAIModel::SampleFluctuations( const G4Material* material,
755 const G4DynamicParticle* aParticle,
756 G4double&,
757 G4double& step,
758 G4double&)
759{
760 size_t jMat;
761 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
762 {
763 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
764 }
765 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
766
767 fPAItransferTable = fPAIxscBank[jMat];
768 fdNdxCutVector = fdNdxCutTable[jMat];
769
770 G4int iTkin, iTransfer, iPlace ;
771 G4long numOfCollisions=0;
772
773 // G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
774 //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl ;
775
776 G4double loss = 0.0, charge2 ;
777 G4double stepSum = 0., stepDelta, lambda, omega;
778 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
779 G4bool numb = true;
780 G4double Tkin = aParticle->GetKineticEnergy() ;
781 G4double MassRatio = fMass/aParticle->GetDefinition()->GetPDGMass() ;
782 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
783 charge2 = charge*charge ;
784 G4double TkinScaled = Tkin*MassRatio ;
785
786 for(iTkin=0;iTkin<fTotBin;iTkin++)
787 {
788 if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
789 }
790 iPlace = iTkin - 1 ;
791 if(iPlace < 0) iPlace = 0;
792 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
793 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
794 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
795
796
797 if(iTkin == fTotBin) // Fermi plato, try from left
798 {
799 meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2;
800 if(meanNumber < 0.) meanNumber = 0. ;
801 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
802 // numOfCollisions = G4Poisson(meanNumber) ;
803 if( meanNumber > 0.) lambda = step/meanNumber;
804 else lambda = DBL_MAX;
805 while(numb)
806 {
807 stepDelta = CLHEP::RandExponential::shoot(lambda);
808 stepSum += stepDelta;
809 if(stepSum >= step) break;
810 numOfCollisions++;
811 }
812 // G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
813
814 while(numOfCollisions)
815 {
816 position = dNdxCut1+
817 ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
818
819 for( iTransfer = 0;
820 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
821 {
822 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
823 }
824 omega = GetEnergyTransfer(iPlace,position,iTransfer);
825 // G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
826 loss += omega;
827 numOfCollisions-- ;
828 }
829 }
830 else
831 {
832 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
833 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
834
835 if(iTkin == 0) // Tkin is too small, trying from right only
836 {
837 meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2;
838 if( meanNumber < 0. ) meanNumber = 0. ;
839 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
840 // numOfCollisions = G4Poisson(meanNumber) ;
841 if( meanNumber > 0.) lambda = step/meanNumber;
842 else lambda = DBL_MAX;
843 while(numb)
844 {
845 stepDelta = CLHEP::RandExponential::shoot(lambda);
846 stepSum += stepDelta;
847 if(stepSum >= step) break;
848 numOfCollisions++;
849 }
850
851 //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
852
853 while(numOfCollisions)
854 {
855 position = dNdxCut2+
856 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
857
858 for( iTransfer = 0;
859 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
860 {
861 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
862 }
863 omega = GetEnergyTransfer(iPlace,position,iTransfer);
864 // G4cout<<omega/keV<<"\t";
865 loss += omega;
866 numOfCollisions-- ;
867 }
868 }
869 else // general case: Tkin between two vectors of the material
870 {
871 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
872 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
873 W = 1.0/(E2 - E1) ;
874 W1 = (E2 - TkinScaled)*W ;
875 W2 = (TkinScaled - E1)*W ;
876
877 // G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
878 // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
879 // G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
880 // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
881
882 meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 +
883 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2;
884 if(meanNumber<0.0) meanNumber = 0.0;
885 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
886 // numOfCollisions = G4Poisson(meanNumber) ;
887 if( meanNumber > 0.) lambda = step/meanNumber;
888 else lambda = DBL_MAX;
889 while(numb)
890 {
891 stepDelta = CLHEP::RandExponential::shoot(lambda);
892 stepSum += stepDelta;
893 if(stepSum >= step) break;
894 numOfCollisions++;
895 }
896
897 //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
898
899 while(numOfCollisions)
900 {
901 position = dNdxCut1*W1 + dNdxCut2*W2 +
902 ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 +
903 dNdxCut2+
904 ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand();
905
906 // G4cout<<position<<"\t" ;
907
908 for( iTransfer = 0;
909 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
910 {
911 if( position >=
912 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
913 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
914 {
915 break ;
916 }
917 }
918 omega = GetEnergyTransfer(iPlace,position,iTransfer);
919 // G4cout<<omega/keV<<"\t";
920 loss += omega;
921 numOfCollisions-- ;
922 }
923 }
924 }
925 // G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
926 // <<step/mm<<" mm"<<G4endl ;
927 if(loss > Tkin) loss=Tkin;
928 if(loss < 0. ) loss = 0.;
929 return loss ;
930
931}
932
933//////////////////////////////////////////////////////////////////////
934//
935// Returns the statistical estimation of the energy loss distribution variance
936//
937
938
939G4double G4PAIModel::Dispersion( const G4Material* material,
940 const G4DynamicParticle* aParticle,
941 G4double& tmax,
942 G4double& step )
943{
944 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
945 for(G4int i = 0 ; i < fMeanNumber; i++)
946 {
947 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
948 sumLoss += loss;
949 sumLoss2 += loss*loss;
950 }
951 meanLoss = sumLoss/fMeanNumber;
952 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
953 return sigma2;
954}
955
956/////////////////////////////////////////////////////////////////////
957
958G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
959 G4double kinEnergy)
960{
961 G4double tmax = kinEnergy;
962 if(p == fElectron) tmax *= 0.5;
963 else if(p != fPositron) {
964 G4double mass = p->GetPDGMass();
965 G4double ratio= electron_mass_c2/mass;
966 G4double gamma= kinEnergy/mass + 1.0;
967 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
968 (1. + 2.0*gamma*ratio + ratio*ratio);
969 }
970 return tmax;
971}
972
973///////////////////////////////////////////////////////////////
974
975void G4PAIModel::DefineForRegion(const G4Region* r)
976{
977 fPAIRegionVector.push_back(r);
978}
979
980//
981//
982/////////////////////////////////////////////////
983
984
985
986
987
988
Note: See TracBrowser for help on using the repository browser.