Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PAIPhotonModel.cc
Go to the documentation of this file.
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$
27//
28// -------------------------------------------------------------------
29//
30// GEANT4 Class
31// File name: G4PAIPhotonModel.cc
32//
33// Author: [email protected] based on G4PAIModel class
34//
35// Creation date: 20.05.2004
36//
37// Modifications:
38//
39// 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary
40// 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary
41// 11.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
42//
43
44#include "G4PAIPhotonModel.hh"
45
46#include "G4PhysicalConstants.hh"
47#include "G4SystemOfUnits.hh"
48
49#include "G4Region.hh"
50#include "G4PhysicsLogVector.hh"
51#include "G4PhysicsFreeVector.hh"
52#include "G4PhysicsTable.hh"
53#include "G4ProductionCutsTable.hh"
54#include "G4MaterialCutsCouple.hh"
55#include "G4MaterialTable.hh"
56#include "G4SandiaTable.hh"
57#include "G4PAIxSection.hh"
58
59#include "Randomize.hh"
60#include "G4Electron.hh"
61#include "G4Positron.hh"
62#include "G4Gamma.hh"
63#include "G4Poisson.hh"
64#include "G4Step.hh"
65#include "G4Material.hh"
66#include "G4DynamicParticle.hh"
67#include "G4ParticleDefinition.hh"
68#include "G4ParticleChangeForLoss.hh"
69#include "G4GeometryTolerance.hh"
70
71////////////////////////////////////////////////////////////////////////
72
73using namespace std;
74
75G4PAIPhotonModel::G4PAIPhotonModel(const G4ParticleDefinition* p, const G4String& nam)
76 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
77 fLowestKineticEnergy(10.0*keV),
78 fHighestKineticEnergy(100.*TeV),
79 fTotBin(200),
80 fMeanNumber(20),
81 fParticle(0),
82 fHighKinEnergy(100.*TeV),
83 fLowKinEnergy(2.0*MeV),
84 fTwoln10(2.0*log(10.0)),
85 fBg2lim(0.0169),
86 fTaulim(8.4146e-3)
87{
88 fVerbose = 0;
89 fElectron = G4Electron::Electron();
90 fPositron = G4Positron::Positron();
91
92 fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
93 fHighestKineticEnergy,
94 fTotBin);
95 fPAItransferTable = 0;
96 fPAIphotonTable = 0;
97 fPAIplasmonTable = 0;
98
99 fPAIdEdxTable = 0;
100 fSandiaPhotoAbsCof = 0;
101 fdEdxVector = 0;
102
103 fLambdaVector = 0;
104 fdNdxCutVector = 0;
105 fdNdxCutPhotonVector = 0;
106 fdNdxCutPlasmonVector = 0;
107
108 fSandiaIntervalNumber = 0;
109 fMatIndex = 0;
110
111 fParticleChange = 0;
112
113 if(p) { SetParticle(p); }
114 else { SetParticle(fElectron); }
115
116 isInitialised = false;
117}
118
119////////////////////////////////////////////////////////////////////////////
120
121G4PAIPhotonModel::~G4PAIPhotonModel()
122{
123 if(fProtonEnergyVector) delete fProtonEnergyVector;
124 if(fdEdxVector) delete fdEdxVector ;
125 if ( fLambdaVector) delete fLambdaVector;
126 if ( fdNdxCutVector) delete fdNdxCutVector;
127
128 if( fPAItransferTable )
129 {
130 fPAItransferTable->clearAndDestroy();
131 delete fPAItransferTable ;
132 }
133 if( fPAIphotonTable )
134 {
135 fPAIphotonTable->clearAndDestroy();
136 delete fPAIphotonTable ;
137 }
138 if( fPAIplasmonTable )
139 {
140 fPAIplasmonTable->clearAndDestroy();
141 delete fPAIplasmonTable ;
142 }
143 if(fSandiaPhotoAbsCof)
144 {
145 for(G4int i=0;i<fSandiaIntervalNumber;i++)
146 {
147 delete[] fSandiaPhotoAbsCof[i];
148 }
149 delete[] fSandiaPhotoAbsCof;
150 }
151}
152
153///////////////////////////////////////////////////////////////////////////////
154
155void G4PAIPhotonModel::SetParticle(const G4ParticleDefinition* p)
156{
157 fParticle = p;
158 fMass = fParticle->GetPDGMass();
159 fSpin = fParticle->GetPDGSpin();
160 G4double q = fParticle->GetPDGCharge()/eplus;
161 fChargeSquare = q*q;
162 fLowKinEnergy *= fMass/proton_mass_c2;
163 fRatio = electron_mass_c2/fMass;
164 fQc = fMass/fRatio;
165}
166
167////////////////////////////////////////////////////////////////////////////
168
169void G4PAIPhotonModel::Initialise(const G4ParticleDefinition* p,
170 const G4DataVector&)
171{
172 // G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl;
173 if(isInitialised) { return; }
174 isInitialised = true;
175
176 if(!fParticle) SetParticle(p);
177
178 fParticleChange = GetParticleChangeForLoss();
179
180 const G4ProductionCutsTable* theCoupleTable =
181 G4ProductionCutsTable::GetProductionCutsTable();
182
183 for(size_t iReg = 0; iReg < fPAIRegionVector.size();++iReg) // region loop
184 {
185 const G4Region* curReg = fPAIRegionVector[iReg];
186
187 vector<G4Material*>::const_iterator matIter = curReg->GetMaterialIterator();
188 size_t jMat;
189 size_t numOfMat = curReg->GetNumberOfMaterials();
190
191 // for(size_t jMat = 0; jMat < curReg->GetNumberOfMaterials();++jMat){}
192 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
193 size_t numberOfMat = G4Material::GetNumberOfMaterials();
194
195 for(jMat = 0 ; jMat < numOfMat; ++jMat) // region material loop
196 {
197 const G4MaterialCutsCouple* matCouple = theCoupleTable->
198 GetMaterialCutsCouple( *matIter, curReg->GetProductionCuts() );
199 fMaterialCutsCoupleVector.push_back(matCouple);
200
201 size_t iMatGlob;
202 for(iMatGlob = 0 ; iMatGlob < numberOfMat ; iMatGlob++ )
203 {
204 if( *matIter == (*theMaterialTable)[iMatGlob]) break ;
205 }
206 fMatIndex = iMatGlob;
207
208 ComputeSandiaPhotoAbsCof();
209 BuildPAIonisationTable();
210
211 fPAIxscBank.push_back(fPAItransferTable);
212 fPAIphotonBank.push_back(fPAIphotonTable);
213 fPAIplasmonBank.push_back(fPAIplasmonTable);
214 fPAIdEdxBank.push_back(fPAIdEdxTable);
215 fdEdxTable.push_back(fdEdxVector);
216
217 BuildLambdaVector(matCouple);
218
219 fdNdxCutTable.push_back(fdNdxCutVector);
220 fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector);
221 fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector);
222 fLambdaTable.push_back(fLambdaVector);
223
224
225 matIter++;
226 }
227 }
228}
229
230//////////////////////////////////////////////////////////////////
231
232void G4PAIPhotonModel::InitialiseMe(const G4ParticleDefinition*)
233{}
234
235//////////////////////////////////////////////////////////////////
236
237void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof()
238{
239 G4int i, j, numberOfElements ;
240 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
241
242 G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
243 numberOfElements = (*theMaterialTable)[fMatIndex]->
244 GetNumberOfElements();
245 G4int* thisMaterialZ = new G4int[numberOfElements] ;
246
247 for(i=0;i<numberOfElements;i++)
248 {
249 thisMaterialZ[i] =
250 (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
251 }
252 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
253 (thisMaterialZ,numberOfElements) ;
254
255 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
256 ( thisMaterialZ ,
257 (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
258 numberOfElements,fSandiaIntervalNumber) ;
259
260 fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
261
262 for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ;
263
264 for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
265 {
266 fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ;
267
268 for( j = 1; j < 5 ; j++ )
269 {
270 fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
271 GetPhotoAbsorpCof(i+1,j)*
272 (*theMaterialTable)[fMatIndex]->GetDensity() ;
273 }
274 }
275 delete[] thisMaterialZ ;
276}
277
278////////////////////////////////////////////////////////////////////////////
279//
280// Build tables for the ionization energy loss
281// the tables are built for MATERIALS
282// *********
283
284void
285G4PAIPhotonModel::BuildPAIonisationTable()
286{
287 G4double LowEdgeEnergy , ionloss ;
288 G4double /*massRatio,*/ tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2 ;
289 /*
290 if( fPAItransferTable )
291 {
292 fPAItransferTable->clearAndDestroy() ;
293 delete fPAItransferTable ;
294 }
295 */
296 fPAItransferTable = new G4PhysicsTable(fTotBin);
297 /*
298 if( fPAIratioTable )
299 {
300 fPAIratioTable->clearAndDestroy() ;
301 delete fPAIratioTable ;
302 }
303 */
304 fPAIphotonTable = new G4PhysicsTable(fTotBin);
305 fPAIplasmonTable = new G4PhysicsTable(fTotBin);
306 /*
307 if( fPAIdEdxTable )
308 {
309 fPAIdEdxTable->clearAndDestroy() ;
310 delete fPAIdEdxTable ;
311 }
312 */
313 fPAIdEdxTable = new G4PhysicsTable(fTotBin);
314
315 // if(fdEdxVector) delete fdEdxVector ;
316 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
317 fHighestKineticEnergy,
318 fTotBin ) ;
319 Tmin = fSandiaPhotoAbsCof[0][0] ; // low energy Sandia interval
320 deltaLow = 100.*eV; // 0.5*eV ;
321
322 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy
323 {
324 LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ;
325 tau = LowEdgeEnergy/proton_mass_c2 ;
326 // if(tau < 0.01) tau = 0.01 ;
327 //gamma = tau +1. ;
328 // G4cout<<"gamma = "<<gamma<<endl ;
329 bg2 = tau*(tau + 2. ) ;
330 //massRatio = electron_mass_c2/proton_mass_c2 ;
331 Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
332 // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
333 // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
334 // Tkin = DeltaCutInKineticEnergyNow ;
335
336 // if ( DeltaCutInKineticEnergyNow > Tmax) // was <
337 Tkin = Tmax ;
338 if ( Tkin < Tmin + deltaLow ) // low energy safety
339 {
340 Tkin = Tmin + deltaLow ;
341 }
342 G4PAIxSection protonPAI( fMatIndex,
343 Tkin,
344 bg2,
345 fSandiaPhotoAbsCof,
346 fSandiaIntervalNumber ) ;
347
348
349 // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
350 // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
351 // G4cout<<"protonPAI.GetSplineSize() = "<<
352 // protonPAI.GetSplineSize()<<G4endl<<G4endl ;
353
354 G4PhysicsFreeVector* transferVector = new
355 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
356 G4PhysicsFreeVector* photonVector = new
357 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
358 G4PhysicsFreeVector* plasmonVector = new
359 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
360 G4PhysicsFreeVector* dEdxVector = new
361 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
362
363 for( G4int k = 0 ; k < protonPAI.GetSplineSize() ; k++ )
364 {
365 transferVector->PutValue( k ,
366 protonPAI.GetSplineEnergy(k+1),
367 protonPAI.GetIntegralPAIxSection(k+1) ) ;
368 photonVector->PutValue( k ,
369 protonPAI.GetSplineEnergy(k+1),
370 protonPAI.GetIntegralCerenkov(k+1) ) ;
371 plasmonVector->PutValue( k ,
372 protonPAI.GetSplineEnergy(k+1),
373 protonPAI.GetIntegralPlasmon(k+1) ) ;
374 dEdxVector->PutValue( k ,
375 protonPAI.GetSplineEnergy(k+1),
376 protonPAI.GetIntegralPAIdEdx(k+1) ) ;
377 }
378 ionloss = protonPAI.GetMeanEnergyLoss() ; // total <dE/dx>
379 if ( ionloss <= 0.) ionloss = DBL_MIN ;
380 fdEdxVector->PutValue(i,ionloss) ;
381
382 fPAItransferTable->insertAt(i,transferVector) ;
383 fPAIphotonTable->insertAt(i,photonVector) ;
384 fPAIplasmonTable->insertAt(i,plasmonVector) ;
385 fPAIdEdxTable->insertAt(i,dEdxVector) ;
386
387 } // end of Tkin loop
388 // theLossTable->insert(fdEdxVector);
389 // end of material loop
390 // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;
391 // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;
392}
393
394///////////////////////////////////////////////////////////////////////
395//
396// Build mean free path tables for the delta ray production process
397// tables are built for MATERIALS
398//
399
400void
401G4PAIPhotonModel::BuildLambdaVector(const G4MaterialCutsCouple* matCutsCouple)
402{
403 G4int i ;
404 G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda;
405 G4double kCarTolerance = G4GeometryTolerance::GetInstance()
406 ->GetSurfaceTolerance();
407
408 const G4ProductionCutsTable* theCoupleTable=
409 G4ProductionCutsTable::GetProductionCutsTable();
410
411 size_t numOfCouples = theCoupleTable->GetTableSize();
412 size_t jMatCC;
413
414 for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ )
415 {
416 if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
417 }
418 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
419
420 const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->
421 GetEnergyCutsVector(idxG4ElectronCut);
422 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
423 GetEnergyCutsVector(idxG4GammaCut);
424
425 if (fLambdaVector) delete fLambdaVector;
426 if (fdNdxCutVector) delete fdNdxCutVector;
427 if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector;
428 if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector;
429
430 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
431 fHighestKineticEnergy,
432 fTotBin );
433 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
434 fHighestKineticEnergy,
435 fTotBin );
436 fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
437 fHighestKineticEnergy,
438 fTotBin );
439 fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
440 fHighestKineticEnergy,
441 fTotBin );
442
443 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
444 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
445
446 if(fVerbose > 0)
447 {
448 G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = "
449 <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
450 G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = "
451 <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
452 }
453 for ( i = 0 ; i <= fTotBin ; i++ )
454 {
455 dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
456 dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow);
457
458 dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
459 lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut;
460
461 if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ???
462
463 fLambdaVector->PutValue(i, lambda);
464
465 fdNdxCutVector->PutValue(i, dNdxCut);
466 fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut);
467 fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut);
468 }
469}
470
471///////////////////////////////////////////////////////////////////////
472//
473// Returns integral PAI cross section for energy transfers >= transferCut
474
475G4double
476G4PAIPhotonModel::GetdNdxCut( G4int iPlace, G4double transferCut)
477{
478 G4int iTransfer;
479 G4double x1, x2, y1, y2, dNdxCut;
480 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
481 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
482 // <<G4endl;
483 for( iTransfer = 0 ;
484 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
485 iTransfer++)
486 {
487 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
488 {
489 break ;
490 }
491 }
492 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
493 {
494 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
495 }
496 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
497 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
498 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
499 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
500 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
501 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
502
503 if ( y1 == y2 ) dNdxCut = y2 ;
504 else
505 {
506 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
507 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
508 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
509 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
510 }
511 // G4cout<<""<<dNdxCut<<G4endl;
512 return dNdxCut ;
513}
514
515///////////////////////////////////////////////////////////////////////
516//
517// Returns integral PAI cherenkovcross section for energy transfers >= transferCut
518
519G4double
520G4PAIPhotonModel::GetdNdxPhotonCut( G4int iPlace, G4double transferCut)
521{
522 G4int iTransfer;
523 G4double x1, x2, y1, y2, dNdxCut;
524 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
525 // G4cout<<"size = "<<G4int((*fPAIphotonTable)(iPlace)->GetVectorLength())
526 // <<G4endl;
527 for( iTransfer = 0 ;
528 iTransfer < G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) ;
529 iTransfer++)
530 {
531 if(transferCut <= (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
532 {
533 break ;
534 }
535 }
536 if ( iTransfer >= G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) )
537 {
538 iTransfer = (*fPAIphotonTable)(iPlace)->GetVectorLength() - 1 ;
539 }
540 y1 = (*(*fPAIphotonTable)(iPlace))(iTransfer-1) ;
541 y2 = (*(*fPAIphotonTable)(iPlace))(iTransfer) ;
542 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
543 x1 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
544 x2 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
545 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
546
547 if ( y1 == y2 ) dNdxCut = y2 ;
548 else
549 {
550 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
551 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
552 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
553 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
554 }
555 // G4cout<<""<<dNdxPhotonCut<<G4endl;
556 return dNdxCut ;
557}
558
559///////////////////////////////////////////////////////////////////////
560//
561// Returns integral PAI cross section for energy transfers >= transferCut
562
563G4double
564G4PAIPhotonModel::GetdNdxPlasmonCut( G4int iPlace, G4double transferCut)
565{
566 G4int iTransfer;
567 G4double x1, x2, y1, y2, dNdxCut;
568
569 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
570 // G4cout<<"size = "<<G4int((*fPAIPlasmonTable)(iPlace)->GetVectorLength())
571 // <<G4endl;
572 for( iTransfer = 0 ;
573 iTransfer < G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) ;
574 iTransfer++)
575 {
576 if(transferCut <= (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
577 {
578 break ;
579 }
580 }
581 if ( iTransfer >= G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) )
582 {
583 iTransfer = (*fPAIplasmonTable)(iPlace)->GetVectorLength() - 1 ;
584 }
585 y1 = (*(*fPAIplasmonTable)(iPlace))(iTransfer-1) ;
586 y2 = (*(*fPAIplasmonTable)(iPlace))(iTransfer) ;
587 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
588 x1 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
589 x2 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
590 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
591
592 if ( y1 == y2 ) dNdxCut = y2 ;
593 else
594 {
595 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
596 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
597 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
598 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
599 }
600 // G4cout<<""<<dNdxPlasmonCut<<G4endl;
601 return dNdxCut ;
602}
603
604///////////////////////////////////////////////////////////////////////
605//
606// Returns integral dEdx for energy transfers >= transferCut
607
608G4double
609G4PAIPhotonModel::GetdEdxCut( G4int iPlace, G4double transferCut)
610{
611 G4int iTransfer;
612 G4double x1, x2, y1, y2, dEdxCut;
613 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
614 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
615 // <<G4endl;
616 for( iTransfer = 0 ;
617 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
618 iTransfer++)
619 {
620 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
621 {
622 break ;
623 }
624 }
625 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
626 {
627 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
628 }
629 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
630 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
631 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
632 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
633 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
634 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
635
636 if ( y1 == y2 ) dEdxCut = y2 ;
637 else
638 {
639 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
640 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
641 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
642 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
643 }
644 // G4cout<<""<<dEdxCut<<G4endl;
645 return dEdxCut ;
646}
647
648//////////////////////////////////////////////////////////////////////////////
649
650G4double G4PAIPhotonModel::ComputeDEDXPerVolume(const G4Material*,
651 const G4ParticleDefinition* p,
652 G4double kineticEnergy,
653 G4double cutEnergy)
654{
655 G4int iTkin,iPlace;
656 size_t jMat;
657
658 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
659 G4double cut = cutEnergy;
660
661 G4double particleMass = p->GetPDGMass();
662 G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
663 G4double charge = p->GetPDGCharge()/eplus;
664 G4double charge2 = charge*charge;
665 G4double dEdx = 0.;
666 const G4MaterialCutsCouple* matCC = CurrentCouple();
667
668 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
669 {
670 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
671 }
672 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
673
674 fPAIdEdxTable = fPAIdEdxBank[jMat];
675 fdEdxVector = fdEdxTable[jMat];
676 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
677 {
678 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
679 }
680 iPlace = iTkin - 1;
681 if(iPlace < 0) iPlace = 0;
682 dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ) ;
683
684 if( dEdx < 0.) dEdx = 0.;
685 return dEdx;
686}
687
688/////////////////////////////////////////////////////////////////////////
689
690G4double G4PAIPhotonModel::CrossSectionPerVolume( const G4Material*,
691 const G4ParticleDefinition* p,
692 G4double kineticEnergy,
693 G4double cutEnergy,
694 G4double maxEnergy )
695{
696 G4int iTkin,iPlace;
697 size_t jMat, jMatCC;
698 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
699 if(cutEnergy >= tmax) return 0.0;
700 G4double particleMass = p->GetPDGMass();
701 G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
702 G4double charge = p->GetPDGCharge();
703 G4double charge2 = charge*charge, cross, cross1, cross2;
704 G4double photon1, photon2, plasmon1, plasmon2;
705
706 const G4MaterialCutsCouple* matCC = CurrentCouple();
707
708 const G4ProductionCutsTable* theCoupleTable=
709 G4ProductionCutsTable::GetProductionCutsTable();
710
711 size_t numOfCouples = theCoupleTable->GetTableSize();
712
713 for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ )
714 {
715 if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
716 }
717 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
718
719 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
720 GetEnergyCutsVector(idxG4GammaCut);
721
722 G4double photonCut = (*photonCutInKineticEnergy)[jMatCC] ;
723
724 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
725 {
726 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
727 }
728 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
729
730 fPAItransferTable = fPAIxscBank[jMat];
731 fPAIphotonTable = fPAIphotonBank[jMat];
732 fPAIplasmonTable = fPAIplasmonBank[jMat];
733
734 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
735 {
736 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
737 }
738 iPlace = iTkin - 1;
739 if(iPlace < 0) iPlace = 0;
740
741 // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
742 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
743 photon1 = GetdNdxPhotonCut(iPlace,tmax);
744 photon2 = GetdNdxPhotonCut(iPlace,photonCut);
745
746 plasmon1 = GetdNdxPlasmonCut(iPlace,tmax);
747 plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy);
748
749 cross1 = photon1 + plasmon1;
750 // G4cout<<"cross1 = "<<cross1<<G4endl;
751 cross2 = photon2 + plasmon2;
752 // G4cout<<"cross2 = "<<cross2<<G4endl;
753 cross = (cross2 - cross1)*charge2;
754 // G4cout<<"cross = "<<cross<<G4endl;
755
756 if( cross < 0. ) cross = 0.;
757 return cross;
758}
759
760///////////////////////////////////////////////////////////////////////////
761//
762// It is analog of PostStepDoIt in terms of secondary electron or photon to
763// be returned as G4Dynamicparticle*.
764//
765
766void G4PAIPhotonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
767 const G4MaterialCutsCouple* matCC,
768 const G4DynamicParticle* dp,
769 G4double tmin,
770 G4double maxEnergy)
771{
772 size_t jMat;
773 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
774 {
775 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
776 }
777 if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--;
778
779 fPAItransferTable = fPAIxscBank[jMat];
780 fPAIphotonTable = fPAIphotonBank[jMat];
781 fPAIplasmonTable = fPAIplasmonBank[jMat];
782
783 fdNdxCutVector = fdNdxCutTable[jMat];
784 fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
785 fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
786
787 G4double tmax = min(MaxSecondaryKinEnergy(dp), maxEnergy);
788 if( tmin >= tmax && fVerbose > 0)
789 {
790 G4cout<<"G4PAIPhotonModel::SampleSecondary: tmin >= tmax "<<G4endl;
791 }
792
793 G4ThreeVector direction = dp->GetMomentumDirection();
794 G4double particleMass = dp->GetMass();
795 G4double kineticEnergy = dp->GetKineticEnergy();
796 G4double scaledTkin = kineticEnergy*fMass/particleMass;
797 G4double totalEnergy = kineticEnergy + particleMass;
798 G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
799
800 G4int iTkin;
801 for(iTkin=0;iTkin<=fTotBin;iTkin++)
802 {
803 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
804 }
805 G4int iPlace = iTkin - 1 ;
806 if(iPlace < 0) iPlace = 0;
807
808 G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace) ;
809 G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace) ;
810 G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
811
812 G4double ratio;
813 if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut;
814 else ratio = 0.;
815
816 if(ratio < G4UniformRand() ) // secondary e-
817 {
818 G4double deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector,
819 iPlace, scaledTkin);
820
821// G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ;
822
823 if( deltaTkin <= 0. )
824 {
825 G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
826 }
827 if( deltaTkin <= 0.) return;
828
829 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
830 G4double totalMomentum = sqrt(pSquare);
831 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
832 /(deltaTotalMomentum * totalMomentum);
833
834 if( costheta > 0.99999 ) costheta = 0.99999;
835 G4double sintheta = 0.0;
836 G4double sin2 = 1. - costheta*costheta;
837 if( sin2 > 0.) sintheta = sqrt(sin2);
838
839 // direction of the delta electron
840
841 G4double phi = twopi*G4UniformRand();
842 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
843
844 G4ThreeVector deltaDirection(dirx,diry,dirz);
845 deltaDirection.rotateUz(direction);
846
847 // primary change
848
849 kineticEnergy -= deltaTkin;
850 G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
851 direction = dir.unit();
852 fParticleChange->SetProposedMomentumDirection(direction);
853
854 // create G4DynamicParticle object for e- delta ray
855
856 G4DynamicParticle* deltaRay = new G4DynamicParticle;
857 deltaRay->SetDefinition(G4Electron::Electron());
858 deltaRay->SetKineticEnergy( deltaTkin );
859 deltaRay->SetMomentumDirection(deltaDirection);
860 vdp->push_back(deltaRay);
861
862 }
863 else // secondary 'Cherenkov' photon
864 {
865 G4double deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector,
866 iPlace,scaledTkin);
867
868 // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ;
869
870 if( deltaTkin <= 0. )
871 {
872 G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
873 }
874 if( deltaTkin <= 0.) return;
875
876 G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
877 G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
878
879 // direction of the 'Cherenkov' photon
880 G4double phi = twopi*G4UniformRand();
881 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
882
883 G4ThreeVector deltaDirection(dirx,diry,dirz);
884 deltaDirection.rotateUz(direction);
885
886 // primary change
887 kineticEnergy -= deltaTkin;
888
889 // create G4DynamicParticle object for photon ray
890
891 G4DynamicParticle* photonRay = new G4DynamicParticle;
892 photonRay->SetDefinition( G4Gamma::Gamma() );
893 photonRay->SetKineticEnergy( deltaTkin );
894 photonRay->SetMomentumDirection(deltaDirection);
895
896 vdp->push_back(photonRay);
897 }
898
899 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
900}
901
902
903///////////////////////////////////////////////////////////////////////
904//
905// Returns post step PAI energy transfer > cut electron/photon energy according to passed
906// scaled kinetic energy of particle
907
908G4double
909G4PAIPhotonModel::GetPostStepTransfer( G4PhysicsTable* pTable,
910 G4PhysicsLogVector* pVector,
911 G4int iPlace, G4double scaledTkin )
912{
913 // G4cout<<"G4PAIPhotonModel::GetPostStepTransfer"<<G4endl ;
914
915 G4int iTkin = iPlace+1, iTransfer;
916 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
917
918 dNdxCut1 = (*pVector)(iPlace) ;
919
920 // G4cout<<"iPlace = "<<iPlace<<endl ;
921
922 if(iTkin == fTotBin) // Fermi plato, try from left
923 {
924 position = dNdxCut1*G4UniformRand() ;
925
926 for( iTransfer = 0;
927 iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
928 {
929 if(position >= (*(*pTable)(iPlace))(iTransfer)) break ;
930 }
931 transfer = GetEnergyTransfer(pTable,iPlace,position,iTransfer);
932 }
933 else
934 {
935 dNdxCut2 = (*pVector)(iPlace+1) ;
936 if(iTkin == 0) // Tkin is too small, trying from right only
937 {
938 position = dNdxCut2*G4UniformRand() ;
939
940 for( iTransfer = 0;
941 iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
942 {
943 if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ;
944 }
945 transfer = GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
946 }
947 else // general case: Tkin between two vectors of the material
948 {
949 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
950 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
951 W = 1.0/(E2 - E1) ;
952 W1 = (E2 - scaledTkin)*W ;
953 W2 = (scaledTkin - E1)*W ;
954
955 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
956
957 // G4cout<<position<<"\t" ;
958
959 G4int iTrMax1, iTrMax2, iTrMax;
960
961 iTrMax1 = G4int((*pTable)(iPlace)->GetVectorLength());
962 iTrMax2 = G4int((*pTable)(iPlace+1)->GetVectorLength());
963
964 if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
965 else iTrMax = iTrMax1;
966
967 for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
968 {
969 if( position >=
970 ( (*(*pTable)(iPlace))(iTransfer)*W1 +
971 (*(*pTable)(iPlace+1))(iTransfer)*W2) ) break ;
972 }
973 transfer = GetEnergyTransfer(pTable, iPlace, position, iTransfer);
974 }
975 }
976 // G4cout<<"PAIPhotonModel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
977 if(transfer < 0.0 ) transfer = 0.0 ;
978 return transfer ;
979}
980
981///////////////////////////////////////////////////////////////////////
982//
983// Returns random PAI energy transfer according to passed
984// indexes of particle
985
986G4double
987G4PAIPhotonModel::GetEnergyTransfer( G4PhysicsTable* pTable, G4int iPlace,
988 G4double position, G4int iTransfer )
989{
990 G4int iTransferMax;
991 G4double x1, x2, y1, y2, energyTransfer;
992
993 if(iTransfer == 0)
994 {
995 energyTransfer = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
996 }
997 else
998 {
999 iTransferMax = G4int((*pTable)(iPlace)->GetVectorLength());
1000
1001 if ( iTransfer >= iTransferMax) iTransfer = iTransferMax - 1;
1002
1003 y1 = (*(*pTable)(iPlace))(iTransfer-1);
1004 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
1005
1006 x1 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1007 x2 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1008
1009 if ( x1 == x2 ) energyTransfer = x2;
1010 else
1011 {
1012 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
1013 else
1014 {
1015 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1016 }
1017 }
1018 }
1019 return energyTransfer;
1020}
1021
1022///////////////////////////////////////////////////////////////////////
1023//
1024// Works like AlongStepDoIt method of process family
1025
1026
1027
1028
1029G4double G4PAIPhotonModel::SampleFluctuations( const G4Material* material,
1030 const G4DynamicParticle* aParticle,
1031 G4double&,
1032 G4double& step,
1033 G4double&)
1034{
1035 size_t jMat;
1036 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
1037 {
1038 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
1039 }
1040 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
1041
1042 fPAItransferTable = fPAIxscBank[jMat];
1043 fPAIphotonTable = fPAIphotonBank[jMat];
1044 fPAIplasmonTable = fPAIplasmonBank[jMat];
1045
1046 fdNdxCutVector = fdNdxCutTable[jMat];
1047 fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1048 fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1049
1050 G4int iTkin, iPlace ;
1051
1052 // G4cout<<"G4PAIPhotonModel::SampleFluctuations"<<G4endl ;
1053
1054 G4double loss, photonLoss, plasmonLoss, charge2 ;
1055
1056
1057 G4double Tkin = aParticle->GetKineticEnergy() ;
1058 G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ;
1059 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
1060 charge2 = charge*charge ;
1061 G4double scaledTkin = Tkin*MassRatio ;
1062 G4double cof = step*charge2;
1063
1064 for( iTkin = 0; iTkin <= fTotBin; iTkin++)
1065 {
1066 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
1067 }
1068 iPlace = iTkin - 1 ;
1069 if( iPlace < 0 ) iPlace = 0;
1070
1071 photonLoss = GetAlongStepTransfer(fPAIphotonTable,fdNdxCutPhotonVector,
1072iPlace,scaledTkin,step,cof);
1073
1074 // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl ;
1075
1076 plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector,
1077iPlace,scaledTkin,step,cof);
1078
1079 // G4cout<<"PAIPhotonModel AlongStepPlasmonLoss = "<<plasmonLoss/keV<<" keV"<<G4endl ;
1080
1081 loss = photonLoss + plasmonLoss;
1082
1083 // G4cout<<"PAIPhotonModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl ;
1084
1085 return loss;
1086}
1087
1088///////////////////////////////////////////////////////////////////////
1089//
1090// Returns along step PAI energy transfer < cut electron/photon energy according to passed
1091// scaled kinetic energy of particle and cof = step*charge*charge
1092
1093G4double
1094G4PAIPhotonModel::GetAlongStepTransfer( G4PhysicsTable* pTable,
1095 G4PhysicsLogVector* pVector,
1096 G4int iPlace, G4double scaledTkin,G4double step,
1097 G4double cof )
1098{
1099 G4int iTkin = iPlace + 1, iTransfer;
1100 G4double loss = 0., position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
1101 G4double lambda, stepDelta, stepSum=0. ;
1102 G4long numOfCollisions=0;
1103 G4bool numb = true;
1104
1105 dNdxCut1 = (*pVector)(iPlace) ;
1106
1107 // G4cout<<"iPlace = "<<iPlace<<endl ;
1108
1109 if(iTkin == fTotBin) // Fermi plato, try from left
1110 {
1111 meanNumber = ((*(*pTable)(iPlace))(0) - dNdxCut1)*cof;
1112 if(meanNumber < 0.) meanNumber = 0. ;
1113 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
1114 if( meanNumber > 0.) lambda = step/meanNumber;
1115 else lambda = DBL_MAX;
1116 while(numb)
1117 {
1118 stepDelta = CLHEP::RandExponential::shoot(lambda);
1119 stepSum += stepDelta;
1120 if(stepSum >= step) break;
1121 numOfCollisions++;
1122 }
1123
1124 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1125
1126 while(numOfCollisions)
1127 {
1128 position = dNdxCut1+
1129 ((*(*pTable)(iPlace))(0) - dNdxCut1)*G4UniformRand() ;
1130
1131 for( iTransfer = 0;
1132 iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1133 {
1134 if(position >= (*(*pTable)(iPlace))(iTransfer)) break ;
1135 }
1136 loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1137 numOfCollisions-- ;
1138 }
1139 }
1140 else
1141 {
1142 dNdxCut2 = (*pVector)(iPlace+1) ;
1143
1144 if(iTkin == 0) // Tkin is too small, trying from right only
1145 {
1146 meanNumber = ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*cof;
1147 if( meanNumber < 0. ) meanNumber = 0. ;
1148 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
1149 if( meanNumber > 0.) lambda = step/meanNumber;
1150 else lambda = DBL_MAX;
1151 while(numb)
1152 {
1153 stepDelta = CLHEP::RandExponential::shoot(lambda);
1154 stepSum += stepDelta;
1155 if(stepSum >= step) break;
1156 numOfCollisions++;
1157 }
1158
1159 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1160
1161 while(numOfCollisions)
1162 {
1163 position = dNdxCut2+
1164 ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*G4UniformRand();
1165
1166 for( iTransfer = 0;
1167 iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
1168 {
1169 if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ;
1170 }
1171 loss += GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
1172 numOfCollisions-- ;
1173 }
1174 }
1175 else // general case: Tkin between two vectors of the material
1176 {
1177 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
1178 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1179 W = 1.0/(E2 - E1) ;
1180 W1 = (E2 - scaledTkin)*W ;
1181 W2 = (scaledTkin - E1)*W ;
1182
1183 // G4cout<<"(*(*pTable)(iPlace))(0) = "<<
1184 // (*(*pTable)(iPlace))(0)<<G4endl ;
1185 // G4cout<<"(*(*pTable)(iPlace+1))(0) = "<<
1186 // (*(*pTable)(iPlace+1))(0)<<G4endl ;
1187
1188 meanNumber=( ((*(*pTable)(iPlace))(0)-dNdxCut1)*W1 +
1189 ((*(*pTable)(iPlace+1))(0)-dNdxCut2)*W2 )*cof;
1190 if(meanNumber<0.0) meanNumber = 0.0;
1191 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
1192 if( meanNumber > 0.) lambda = step/meanNumber;
1193 else lambda = DBL_MAX;
1194 while(numb)
1195 {
1196 stepDelta = CLHEP::RandExponential::shoot(lambda);
1197 stepSum += stepDelta;
1198 if(stepSum >= step) break;
1199 numOfCollisions++;
1200 }
1201
1202 // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl ;
1203
1204 while(numOfCollisions)
1205 {
1206 position = dNdxCut1*W1 + dNdxCut2*W2 +
1207 ( ( (*(*pTable)(iPlace ))(0) - dNdxCut1)*W1 +
1208
1209 ( (*(*pTable)(iPlace+1))(0) - dNdxCut2)*W2 )*G4UniformRand();
1210
1211 // G4cout<<position<<"\t" ;
1212
1213 for( iTransfer = 0;
1214 iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1215 {
1216 if( position >=
1217 ( (*(*pTable)(iPlace))(iTransfer)*W1 +
1218 (*(*pTable)(iPlace+1))(iTransfer)*W2) )
1219 {
1220 break ;
1221 }
1222 }
1223 // loss += (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1224 loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1225 numOfCollisions-- ;
1226 }
1227 }
1228 }
1229
1230 return loss ;
1231
1232}
1233
1234//////////////////////////////////////////////////////////////////////
1235//
1236// Returns the statistical estimation of the energy loss distribution variance
1237//
1238
1239
1240G4double G4PAIPhotonModel::Dispersion( const G4Material* material,
1241 const G4DynamicParticle* aParticle,
1242 G4double& tmax,
1243 G4double& step )
1244{
1245 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
1246 for(G4int i = 0 ; i < fMeanNumber; i++)
1247 {
1248 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
1249 sumLoss += loss;
1250 sumLoss2 += loss*loss;
1251 }
1252 meanLoss = sumLoss/fMeanNumber;
1253 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
1254 return sigma2;
1255}
1256
1257/////////////////////////////////////////////////////////////////////
1258
1259G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
1260 G4double kinEnergy)
1261{
1262 G4double tmax = kinEnergy;
1263 if(p == fElectron) tmax *= 0.5;
1264 else if(p != fPositron) {
1265 G4double mass = p->GetPDGMass();
1266 G4double ratio= electron_mass_c2/mass;
1267 G4double gamma= kinEnergy/mass + 1.0;
1268 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
1269 (1. + 2.0*gamma*ratio + ratio*ratio);
1270 }
1271 return tmax;
1272}
1273
1274///////////////////////////////////////////////////////////////
1275
1276void G4PAIPhotonModel::DefineForRegion(const G4Region* r)
1277{
1278 fPAIRegionVector.push_back(r);
1279}
1280
1281
1282//
1283//
1284/////////////////////////////////////////////////
1285
1286
1287
1288
1289
1290
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:43
idxG4ElectronCut
@ idxG4ElectronCut
Definition: G4ProductionCuts.hh:54
idxG4GammaCut
@ idxG4GammaCut
Definition: G4ProductionCuts.hh:53
G4double
double G4double
Definition: G4Types.hh:64
G4long
long G4long
Definition: G4Types.hh:68
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4DynamicParticle::GetMass
G4double GetMass() const
G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:271
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::SetKineticEnergy
void SetKineticEnergy(G4double aEnergy)
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:87
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:73
G4Material::GetMaterialTable
static const G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:562
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:569
G4PAIPhotonModel::BuildLambdaVector
void BuildLambdaVector(const G4MaterialCutsCouple *matCutsCouple)
Definition: G4PAIPhotonModel.cc:401
G4PAIPhotonModel::GetdNdxCut
G4double GetdNdxCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:476
G4PAIPhotonModel::DefineForRegion
void DefineForRegion(const G4Region *r)
Definition: G4PAIPhotonModel.cc:1276
G4PAIPhotonModel::G4PAIPhotonModel
G4PAIPhotonModel(const G4ParticleDefinition *p=0, const G4String &nam="PAIPhoton")
Definition: G4PAIPhotonModel.cc:75
G4PAIPhotonModel::GetdNdxPlasmonCut
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:564
G4PAIPhotonModel::SampleFluctuations
virtual G4double SampleFluctuations(const G4Material *, const G4DynamicParticle *, G4double &, G4double &, G4double &)
Definition: G4PAIPhotonModel.cc:1029
G4PAIPhotonModel::GetdNdxPhotonCut
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:520
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1259
G4PAIPhotonModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4PAIPhotonModel.cc:169
G4PAIPhotonModel::InitialiseMe
virtual void InitialiseMe(const G4ParticleDefinition *)
Definition: G4PAIPhotonModel.cc:232
G4PAIPhotonModel::~G4PAIPhotonModel
virtual ~G4PAIPhotonModel()
Definition: G4PAIPhotonModel.cc:121
G4PAIPhotonModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Definition: G4PAIPhotonModel.cc:690
G4PAIPhotonModel::GetPostStepTransfer
G4double GetPostStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin)
Definition: G4PAIPhotonModel.cc:909
G4PAIPhotonModel::GetAlongStepTransfer
G4double GetAlongStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin, G4double step, G4double cof)
Definition: G4PAIPhotonModel.cc:1094
G4PAIPhotonModel::GetEnergyTransfer
G4double GetEnergyTransfer(G4PhysicsTable *, G4int iPlace, G4double position, G4int iTransfer)
Definition: G4PAIPhotonModel.cc:987
G4PAIPhotonModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4PAIPhotonModel.cc:766
G4PAIPhotonModel::Dispersion
virtual G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double &, G4double &)
Definition: G4PAIPhotonModel.cc:1240
G4PAIPhotonModel::ComputeDEDXPerVolume
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
Definition: G4PAIPhotonModel.cc:650
G4PAIPhotonModel::GetdEdxCut
G4double GetdEdxCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:609
G4PAIxSection::GetIntegralCerenkov
G4double GetIntegralCerenkov(G4int i) const
Definition: G4PAIxSection.hh:316
G4PAIxSection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIxSection.hh:179
G4PAIxSection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition: G4PAIxSection.hh:310
G4PAIxSection::GetIntegralPlasmon
G4double GetIntegralPlasmon(G4int i) const
Definition: G4PAIxSection.hh:328
G4PAIxSection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition: G4PAIxSection.hh:298
G4PAIxSection::GetIntegralPAIxSection
G4double GetIntegralPAIxSection(G4int i) const
Definition: G4PAIxSection.hh:304
G4PAIxSection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition: G4PAIxSection.hh:191
G4ParticleChangeForLoss::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:150
G4ParticleChangeForLoss::SetProposedMomentumDirection
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:194
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119
G4PhysicsFreeVector::PutValue
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Definition: G4PhysicsFreeVector.cc:94
G4PhysicsTable::insertAt
void insertAt(size_t, G4PhysicsVector *)
G4PhysicsTable::clearAndDestroy
void clearAndDestroy()
G4PhysicsVector::GetLowEdgeEnergy
virtual G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:158
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4Region::GetProductionCuts
G4ProductionCuts * GetProductionCuts() const
G4Region::GetNumberOfMaterials
size_t GetNumberOfMaterials() const
G4Region::GetMaterialIterator
std::vector< G4Material * >::const_iterator GetMaterialIterator() const
G4SandiaTable::SandiaMixing
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
Definition: G4SandiaTable.cc:666
G4SandiaTable::SandiaIntervals
G4int SandiaIntervals(G4int Z[], G4int el)
Definition: G4SandiaTable.cc:577
G4SandiaTable::GetPhotoAbsorpCof
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
Definition: G4SandiaTable.hh:211
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:377
G4VEmModel::MaxSecondaryKinEnergy
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:399
G4VEmModel::GetParticleChangeForLoss
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:95
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83