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