source: trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc@ 1344

Last change on this file since 1344 was 1340, checked in by garnier, 15 years ago

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