Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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: G4PAIPhotonModel.cc 76879 2013-11-18 12:44:54Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 // File name: G4PAIPhotonModel.cc
32 //
33 // Author: Vladimir.Grichine@cern.ch 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 // 12.06.13 V. Grichine Bug fixed in SampleSecondaries for scaled Tkin (fMass -> proton_mass_c2)
43 //
44 //
45 
46 #include "G4PAIPhotonModel.hh"
47 
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 
51 #include "G4Region.hh"
52 #include "G4PhysicsLogVector.hh"
53 #include "G4PhysicsFreeVector.hh"
54 #include "G4PhysicsTable.hh"
55 #include "G4ProductionCutsTable.hh"
56 #include "G4EmCalculator.hh"
57 #include "G4ProductionCuts.hh"
58 #include "G4MaterialCutsCouple.hh"
59 #include "G4MaterialTable.hh"
60 
61 #include "G4SandiaTable.hh"
62 #include "G4PAIxSection.hh"
63 
64 #include "Randomize.hh"
65 #include "G4Electron.hh"
66 #include "G4Positron.hh"
67 #include "G4Gamma.hh"
68 #include "G4Poisson.hh"
69 #include "G4Step.hh"
70 #include "G4Material.hh"
71 #include "G4DynamicParticle.hh"
72 #include "G4ParticleDefinition.hh"
73 #include "G4ParticleChangeForLoss.hh"
74 #include "G4GeometryTolerance.hh"
75 
76 ////////////////////////////////////////////////////////////////////////
77 
78 using namespace std;
79 
80 G4PAIPhotonModel::G4PAIPhotonModel(const G4ParticleDefinition* p, const G4String& nam)
81  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
82  fLowestKineticEnergy(10.0*keV),
83  fHighestKineticEnergy(100.*TeV),
84  fTotBin(200),
85  fMeanNumber(20),
86  fParticle(0),
87  fHighKinEnergy(100.*TeV),
88  fLowKinEnergy(2.0*MeV),
89  fTwoln10(2.0*log(10.0)),
90  fBg2lim(0.0169),
91  fTaulim(8.4146e-3)
92 {
93  fVerbose = 0;
94  fGamma = G4Gamma::Gamma();
95  fElectron = G4Electron::Electron();
96  fPositron = G4Positron::Positron();
97 
98  fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
99  fHighestKineticEnergy,
100  fTotBin);
101  fPAItransferTable = 0;
102  fPAIphotonTable = 0;
103  fPAIplasmonTable = 0;
104 
105  fPAIdEdxTable = 0;
106  fSandiaPhotoAbsCof = 0;
107  fdEdxVector = 0;
108 
109  fLambdaVector = 0;
110  fdNdxCutVector = 0;
111  fdNdxCutPhotonVector = 0;
112  fdNdxCutPlasmonVector = 0;
113 
114  fSandiaIntervalNumber = 0;
115  fMatIndex = 0;
116  fCutCouple = 0;
117  fMaterial = 0;
118 
119  fParticleChange = 0;
120 
121  if(p) { SetParticle(p); }
122  else { SetParticle(fElectron); }
123 
124  isInitialised = false;
125 }
126 
127 ////////////////////////////////////////////////////////////////////////////
128 
129 G4PAIPhotonModel::~G4PAIPhotonModel()
130 {
131  // if(fdEdxVector) delete fdEdxVector;
132  // if ( fLambdaVector) delete fLambdaVector;
133  // if ( fdNdxCutVector) delete fdNdxCutVector;
134 
135  if( fPAItransferTable )
136  {
137  delete fPAItransferTable;
138  }
139  if( fPAIphotonTable )
140  {
141  delete fPAIphotonTable;
142  }
143  if( fPAIplasmonTable )
144  {
145  delete fPAIplasmonTable;
146  }
147 }
148 
149 ///////////////////////////////////////////////////////////////////////////////
150 
151 void G4PAIPhotonModel::SetParticle(const G4ParticleDefinition* p)
152 {
153  fParticle = p;
154  fMass = fParticle->GetPDGMass();
155  fSpin = fParticle->GetPDGSpin();
156  G4double q = fParticle->GetPDGCharge()/eplus;
157  fChargeSquare = q*q;
158  fLowKinEnergy *= fMass/proton_mass_c2;
159  fRatio = electron_mass_c2/fMass;
160  fQc = fMass/fRatio;
161 }
162 
163 ////////////////////////////////////////////////////////////////////////////
164 
165 void G4PAIPhotonModel::Initialise(const G4ParticleDefinition* p,
166  const G4DataVector&)
167 {
168  // G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl;
169  if(isInitialised) { return; }
170  isInitialised = true;
171 
172  if(!fParticle) SetParticle(p);
173 
174  fParticleChange = GetParticleChangeForLoss();
175 
176  //const G4ProductionCutsTable* theCoupleTable =
177  // G4ProductionCutsTable::GetProductionCutsTable();
178  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
179  size_t numOfMat = G4Material::GetNumberOfMaterials();
180  size_t numRegions = fPAIRegionVector.size();
181 
182  for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
183  {
184  const G4Region* curReg = fPAIRegionVector[iReg];
185  G4Region* reg = const_cast<G4Region*>(curReg);
186 
187  for(size_t jMat = 0; jMat < numOfMat; ++jMat) // material loop
188  {
189  G4Material* material = (*theMaterialTable)[jMat];
190  const G4MaterialCutsCouple* couple = reg->FindCouple(material);
191 
192  // theCoupleTable->GetMaterialCutsCouple( material,
193  // curReg->GetProductionCuts() );
194 
195  if(couple) {
196  //G4cout << "Reg <" <<curReg->GetName() << "> mat <"
197  // << material->GetName() << "> fCouple= "
198  // << couple<<" " << p->GetParticleName() <<G4endl;
199 
200  fMaterialCutsCoupleVector.push_back(couple);
201 
202  fMatIndex = jMat;
203  fMaterial = material;
204 
205  // ComputeSandiaPhotoAbsCof();
206  fSandia.Initialize(material);
207  BuildPAIonisationTable();
208  /*
209  if( fSandiaPhotoAbsCof ) // delete SANDIA cofs've been used for pai-xsc
210  {
211  for( G4int i = 0;i < fSandiaIntervalNumber;i++)
212  {
213  delete[] fSandiaPhotoAbsCof[i];
214  }
215  delete[] fSandiaPhotoAbsCof;
216  fSandiaPhotoAbsCof = 0;
217  }
218  */
219  fPAIxscBank.push_back(fPAItransferTable);
220  fPAIphotonBank.push_back(fPAIphotonTable);
221  fPAIplasmonBank.push_back(fPAIplasmonTable);
222  fPAIdEdxBank.push_back(fPAIdEdxTable);
223  fdEdxTable.push_back(fdEdxVector);
224 
225  BuildLambdaVector(couple);
226 
227  fdNdxCutTable.push_back(fdNdxCutVector);
228  fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector);
229  fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector);
230  fLambdaTable.push_back(fLambdaVector);
231  }
232  }
233  }
234 }
235 
236 //////////////////////////////////////////////////////////////////
237 
238 void G4PAIPhotonModel::InitTest(const G4ParticleDefinition* p, G4MaterialCutsCouple* couple,
239  G4double phE, G4double eTkin)
240 {
241  // G4cout<<"G4PAIPhotonModel::InitTest for "<<p->GetParticleName()<<G4endl;
242  if(isInitialised) { return; }
243  isInitialised = true;
244 
245  if( !fParticle ) SetParticle(p);
246 
247  fParticleChange = GetParticleChangeForLoss();
248 
249  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
250 
251  size_t jMat, numOfMat = G4Material::GetNumberOfMaterials();
252 
253 
254  // const G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material,cuts);
255 
256  if( couple )
257  {
258  const G4Material* material = couple->GetMaterial();
259 
260  fMaterialCutsCoupleVector.push_back(couple);
261 
262  for( jMat = 0; jMat < numOfMat; ++jMat ) // material loop
263  {
264  if( material->GetName() == (*theMaterialTable)[jMat]->GetName() ) break;
265  }
266  fMatIndex = jMat;
267  G4Material* mat = (*theMaterialTable)[jMat];
268  fMaterial = mat;
269 
270 
271  // ComputeSandiaPhotoAbsCof();
272  fSandia.Initialize(mat);
273  BuildPAIonisationTable();
274  /*
275  if( fSandiaPhotoAbsCof ) // delete SANDIA cofs've been used for pai-xsc
276  {
277  for( G4int i = 0;i < fSandiaIntervalNumber;i++)
278  {
279  delete[] fSandiaPhotoAbsCof[i];
280  }
281  delete[] fSandiaPhotoAbsCof;
282  fSandiaPhotoAbsCof = 0;
283  }
284  */
285  fPAIxscBank.push_back(fPAItransferTable);
286  fPAIphotonBank.push_back(fPAIphotonTable);
287  fPAIplasmonBank.push_back(fPAIplasmonTable);
288  fPAIdEdxBank.push_back(fPAIdEdxTable);
289  fdEdxTable.push_back(fdEdxVector);
290 
291  BuildLambdaVector(couple,phE,eTkin);
292 
293  fdNdxCutTable.push_back(fdNdxCutVector);
294  fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector);
295  fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector);
296  fLambdaTable.push_back(fLambdaVector);
297  }
298 }
299 
300 //////////////////////////////////////////////////////////////////
301 
302 void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof()
303 {
304  G4int i, j, numberOfElements;
305  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
306 
307  G4SandiaTable thisMaterialSandiaTable(fMatIndex);
308  numberOfElements = (*theMaterialTable)[fMatIndex]->
309  GetNumberOfElements();
310  G4int* thisMaterialZ = new G4int[numberOfElements];
311 
312  for(i=0;i<numberOfElements;i++)
313  {
314  thisMaterialZ[i] =
315  (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ();
316  }
317  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
318  (thisMaterialZ,numberOfElements);
319 
320  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
321  ( thisMaterialZ ,
322  (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
323  numberOfElements,fSandiaIntervalNumber);
324 
325  fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber];
326 
327  for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5];
328 
329  for( i = 0; i < fSandiaIntervalNumber; i++ )
330  {
331  fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0);
332 
333  for( j = 1; j < 5; j++ )
334  {
335  fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
336  GetPhotoAbsorpCof(i+1,j)*
337  (*theMaterialTable)[fMatIndex]->GetDensity();
338  }
339  }
340  delete[] thisMaterialZ;
341 }
342 
343 ////////////////////////////////////////////////////////////////////////////
344 //
345 // Build tables for the ionization energy loss
346 // the tables are built for MATERIALS
347 // *********
348 
349 void
350 G4PAIPhotonModel::BuildPAIonisationTable()
351 {
352  G4double LowEdgeEnergy , ionloss;
353  G4double /*massRatio,*/ tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2;
354  /*
355  if( fPAItransferTable )
356  {
357  fPAItransferTable->clearAndDestroy();
358  delete fPAItransferTable;
359  }
360  */
361  fPAItransferTable = new G4PhysicsTable(fTotBin);
362  /*
363  if( fPAIratioTable )
364  {
365  fPAIratioTable->clearAndDestroy();
366  delete fPAIratioTable;
367  }
368  */
369  fPAIphotonTable = new G4PhysicsTable(fTotBin);
370  fPAIplasmonTable = new G4PhysicsTable(fTotBin);
371  /*
372  if( fPAIdEdxTable )
373  {
374  fPAIdEdxTable->clearAndDestroy();
375  delete fPAIdEdxTable;
376  }
377  */
378  fPAIdEdxTable = new G4PhysicsTable(fTotBin);
379 
380  // if(fdEdxVector) delete fdEdxVector;
381  fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
382  fHighestKineticEnergy,
383  fTotBin );
384  // Tmin = fSandiaPhotoAbsCof[0][0]; // low energy Sandia interval
385  Tmin = fSandia.GetSandiaMatTablePAI(0,0); // low energy Sandia interval
386  deltaLow = 100.*eV; // 0.5*eV;
387 
388  for (G4int i = 0; i <= fTotBin; i++) //The loop for the kinetic energy
389  {
390  LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i);
391  tau = LowEdgeEnergy/proton_mass_c2;
392  // if(tau < 0.01) tau = 0.01;
393  //gamma = tau +1.;
394  // G4cout<<"gamma = "<<gamma<<endl;
395  bg2 = tau*(tau + 2. );
396  //massRatio = electron_mass_c2/proton_mass_c2;
397  Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
398  // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
399  // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
400  // Tkin = DeltaCutInKineticEnergyNow;
401 
402  // if ( DeltaCutInKineticEnergyNow > Tmax) // was <
403  Tkin = Tmax;
404  if ( Tkin < Tmin + deltaLow ) // low energy safety
405  {
406  Tkin = Tmin + deltaLow;
407  }
408  fPAIxSection.Initialize(fMaterial, Tkin, bg2,
409  &fSandia);
410  /*
411  G4PAIxSection protonPAI( fMatIndex,
412  Tkin,
413  bg2,
414  fSandiaPhotoAbsCof,
415  fSandiaIntervalNumber );
416 
417  */
418  // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl;
419  // G4cout<<"n1 = "<<fPAIxSection.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl;
420  // G4cout<<"fPAIxSection.GetSplineSize() = "<<
421  // fPAIxSection.GetSplineSize()<<G4endl<<G4endl;
422 
423  G4PhysicsFreeVector* transferVector = new
424  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
425  G4PhysicsFreeVector* photonVector = new
426  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
427  G4PhysicsFreeVector* plasmonVector = new
428  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
429  G4PhysicsFreeVector* dEdxVector = new
430  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
431 
432  for( G4int k = 0; k < fPAIxSection.GetSplineSize(); k++ )
433  {
434  transferVector->PutValue( k ,
435  fPAIxSection.GetSplineEnergy(k+1),
436  fPAIxSection.GetIntegralPAIxSection(k+1) );
437  photonVector->PutValue( k ,
438  fPAIxSection.GetSplineEnergy(k+1),
439  fPAIxSection.GetIntegralCerenkov(k+1) );
440  plasmonVector->PutValue( k ,
441  fPAIxSection.GetSplineEnergy(k+1),
442  fPAIxSection.GetIntegralPlasmon(k+1) );
443  dEdxVector->PutValue( k ,
444  fPAIxSection.GetSplineEnergy(k+1),
445  fPAIxSection.GetIntegralPAIdEdx(k+1) );
446  }
447  ionloss = fPAIxSection.GetMeanEnergyLoss(); // total <dE/dx>
448  if ( ionloss <= 0.) ionloss = DBL_MIN;
449  fdEdxVector->PutValue(i,ionloss);
450 
451  fPAItransferTable->insertAt(i,transferVector);
452  fPAIphotonTable->insertAt(i,photonVector);
453  fPAIplasmonTable->insertAt(i,plasmonVector);
454  fPAIdEdxTable->insertAt(i,dEdxVector);
455 
456  } // end of Tkin loop
457  // theLossTable->insert(fdEdxVector);
458  // end of material loop
459  // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl;
460  // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl;
461 }
462 
463 ///////////////////////////////////////////////////////////////////////
464 //
465 // Build mean free path tables for the delta ray production process
466 // tables are built for MATERIALS
467 //
468 
469 void
470 G4PAIPhotonModel::BuildLambdaVector(const G4MaterialCutsCouple* matCutsCouple)
471 {
472  G4int i;
473  G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda, deltaCutInKineticEnergyNow, photonCutInKineticEnergyNow;
474  G4double kCarTolerance = G4GeometryTolerance::GetInstance()
475  ->GetSurfaceTolerance();
476 
477  G4ProductionCutsTable* theCoupleTable=
478  G4ProductionCutsTable::GetProductionCutsTable();
479 
480  // G4EmCalculator converter;
481 
482  // const G4Material* material = matCutsCouple->GetMaterial();
483 
484  // G4double rangeGamma = matCutsCouple->GetProductionCuts()->GetProductionCut(0);
485  // G4double rangeElectron = matCutsCouple->GetProductionCuts()->GetProductionCut(1);
486 
487  size_t numOfCouples = theCoupleTable->GetTableSize();
488  size_t jMatCC;
489 
490  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
491  {
492  if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
493  }
494  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
495 
496  const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
497  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
498 
499  if (fLambdaVector) delete fLambdaVector;
500  if (fdNdxCutVector) delete fdNdxCutVector;
501  if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector;
502  if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector;
503 
504  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
505  fHighestKineticEnergy,
506  fTotBin );
507  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
508  fHighestKineticEnergy,
509  fTotBin );
510  fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
511  fHighestKineticEnergy,
512  fTotBin );
513  fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
514  fHighestKineticEnergy,
515  fTotBin );
516 
517  deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
518  photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
519 
520  // deltaCutInKineticEnergyNow = theCoupleTable->ConvertRangeToEnergy(fElectron,material,rangeElectron);
521  // photonCutInKineticEnergyNow = theCoupleTable->ConvertRangeToEnergy(fGamma,material,rangeGamma);
522 
523  // photonCutInKineticEnergyNow = converter.GetKinEnergy(rangeGamma, fGamma,material);
524  // deltaCutInKineticEnergyNow = converter.GetKinEnergy(rangeElectron, fElectron,material);
525 
526  if(fVerbose > 0)
527  {
528  G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = "
529  <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
530  G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = "
531  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
532  }
533  for ( i = 0; i <= fTotBin; i++ )
534  {
535  dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
536  dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow);
537 
538  dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
539  lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut;
540 
541  if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ???
542 
543  fLambdaVector->PutValue(i, lambda);
544 
545  fdNdxCutVector->PutValue(i, dNdxCut);
546  fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut);
547  fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut);
548  }
549 }
550 
551 ///////////////////////////////////////////////////////////////////////
552 //
553 // Build mean free path tables for the delta ray production process
554 // tables are built for MATERIALS
555 //
556 
557 void
558 G4PAIPhotonModel::BuildLambdaVector(const G4MaterialCutsCouple* matCutsCouple, G4double photEnergy, G4double eTkin)
559 {
560  G4int i;
561  G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda, deltaCutInKineticEnergyNow, photonCutInKineticEnergyNow;
562  G4double kCarTolerance = G4GeometryTolerance::GetInstance()
563  ->GetSurfaceTolerance();
564 
565  G4ProductionCutsTable* theCoupleTable=
566  G4ProductionCutsTable::GetProductionCutsTable();
567 
568 
569  // const G4Material* material = matCutsCouple->GetMaterial();
570 
571 
572  size_t numOfCouples = theCoupleTable->GetTableSize();
573  size_t jMatCC;
574 
575  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
576  {
577  if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
578  }
579  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
580 
581  if (fLambdaVector) delete fLambdaVector;
582  if (fdNdxCutVector) delete fdNdxCutVector;
583  if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector;
584  if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector;
585 
586  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
587  fHighestKineticEnergy,
588  fTotBin );
589  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
590  fHighestKineticEnergy,
591  fTotBin );
592  fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
593  fHighestKineticEnergy,
594  fTotBin );
595  fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
596  fHighestKineticEnergy,
597  fTotBin );
598 
599 
600  photonCutInKineticEnergyNow = photEnergy;
601  deltaCutInKineticEnergyNow = eTkin;
602 
603  if(fVerbose > 0)
604  {
605  G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = "
606  <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
607  G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = "
608  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
609  }
610  for ( i = 0; i <= fTotBin; i++ )
611  {
612  dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
613  dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow);
614 
615  dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
616  lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut;
617 
618  if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ???
619 
620  fLambdaVector->PutValue(i, lambda);
621 
622  fdNdxCutVector->PutValue(i, dNdxCut);
623  fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut);
624  fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut);
625  }
626 }
627 
628 ///////////////////////////////////////////////////////////////////////
629 //
630 // Returns integral PAI cross section for energy transfers >= transferCut
631 
632 G4double
633 G4PAIPhotonModel::GetdNdxCut( G4int iPlace, G4double transferCut)
634 {
635  G4int iTransfer;
636  G4double x1, x2, y1, y2, dNdxCut;
637  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
638  // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
639  // <<G4endl;
640  for( iTransfer = 0;
641  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
642  iTransfer++)
643  {
644  if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
645  {
646  break;
647  }
648  }
649  if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
650  {
651  iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1;
652  }
653  if (iTransfer == 0) return (*(*fPAItransferTable)(iPlace))(iTransfer);
654  y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
655  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
656  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
657  x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
658  x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
659  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
660 
661  if ( y1 == y2 ) dNdxCut = y2;
662  else
663  {
664  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
665  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
666  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5;
667  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
668  }
669  // G4cout<<""<<dNdxCut<<G4endl;
670  return dNdxCut;
671 }
672 
673 ///////////////////////////////////////////////////////////////////////
674 //
675 // Returns integral PAI cherenkovcross section for energy transfers >= transferCut
676 
677 G4double
678 G4PAIPhotonModel::GetdNdxPhotonCut( G4int iPlace, G4double transferCut)
679 {
680  G4int iTransfer;
681  G4double x1, x2, y1, y2, dNdxCut;
682  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
683  // G4cout<<"size = "<<G4int((*fPAIphotonTable)(iPlace)->GetVectorLength())
684  // <<G4endl;
685  for( iTransfer = 0;
686  iTransfer < G4int((*fPAIphotonTable)(iPlace)->GetVectorLength());
687  iTransfer++)
688  {
689  if(transferCut <= (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
690  {
691  break;
692  }
693  }
694  if ( iTransfer >= G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) )
695  {
696  iTransfer = (*fPAIphotonTable)(iPlace)->GetVectorLength() - 1;
697  }
698  if (iTransfer == 0) return (*(*fPAIphotonTable)(iPlace))(iTransfer);
699  y1 = (*(*fPAIphotonTable)(iPlace))(iTransfer-1);
700  y2 = (*(*fPAIphotonTable)(iPlace))(iTransfer);
701  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
702  x1 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
703  x2 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
704  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
705 
706  if ( y1 == y2 ) dNdxCut = y2;
707  else
708  {
709  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
710  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
711  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5;
712  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
713  }
714  // G4cout<<""<<dNdxPhotonCut<<G4endl;
715  return dNdxCut;
716 }
717 
718 ///////////////////////////////////////////////////////////////////////
719 //
720 // Returns integral PAI cross section for energy transfers >= transferCut
721 
722 G4double
723 G4PAIPhotonModel::GetdNdxPlasmonCut( G4int iPlace, G4double transferCut)
724 {
725  G4int iTransfer;
726  G4double x1, x2, y1, y2, dNdxCut;
727 
728  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
729  // G4cout<<"size = "<<G4int((*fPAIPlasmonTable)(iPlace)->GetVectorLength())
730  // <<G4endl;
731  for( iTransfer = 0;
732  iTransfer < G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength());
733  iTransfer++)
734  {
735  if(transferCut <= (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
736  {
737  break;
738  }
739  }
740  if ( iTransfer >= G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) )
741  {
742  iTransfer = (*fPAIplasmonTable)(iPlace)->GetVectorLength() - 1;
743  }
744  if (iTransfer == 0) return (*(*fPAIplasmonTable)(iPlace))(iTransfer);
745  y1 = (*(*fPAIplasmonTable)(iPlace))(iTransfer-1);
746  y2 = (*(*fPAIplasmonTable)(iPlace))(iTransfer);
747  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
748  x1 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
749  x2 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
750  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
751 
752  if ( y1 == y2 ) dNdxCut = y2;
753  else
754  {
755  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
756  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
757  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5;
758  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
759  }
760  // G4cout<<""<<dNdxPlasmonCut<<G4endl;
761  return dNdxCut;
762 }
763 
764 ///////////////////////////////////////////////////////////////////////
765 //
766 // Returns integral dEdx for energy transfers >= transferCut
767 
768 G4double
769 G4PAIPhotonModel::GetdEdxCut( G4int iPlace, G4double transferCut)
770 {
771  G4int iTransfer;
772  G4double x1, x2, y1, y2, dEdxCut;
773  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
774  // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
775  // <<G4endl;
776  for( iTransfer = 0;
777  iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength());
778  iTransfer++)
779  {
780  if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
781  {
782  break;
783  }
784  }
785  if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
786  {
787  iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1;
788  }
789  if (iTransfer == 0) return (*(*fPAIdEdxTable)(iPlace))(iTransfer);
790  y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1);
791  y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer);
792  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
793  x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
794  x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
795  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
796 
797  if ( y1 == y2 ) dEdxCut = y2;
798  else
799  {
800  // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand();
801  // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand();
802  if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5;
803  else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
804  }
805  // G4cout<<""<<dEdxCut<<G4endl;
806  return dEdxCut;
807 }
808 
809 //////////////////////////////////////////////////////////////////////////////
810 
811 G4double G4PAIPhotonModel::ComputeDEDXPerVolume(const G4Material*,
812  const G4ParticleDefinition* p,
813  G4double kineticEnergy,
814  G4double cutEnergy)
815 {
816  G4int iTkin,iPlace;
817  size_t jMat;
818 
819  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
820  G4double cut = cutEnergy;
821 
822  G4double particleMass = p->GetPDGMass();
823  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
824  G4double charge = p->GetPDGCharge()/eplus;
825  G4double charge2 = charge*charge;
826  G4double dEdx = 0.;
827  const G4MaterialCutsCouple* matCC = CurrentCouple();
828 
829  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
830  {
831  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
832  }
833  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
834  /*
835  G4cout << "G4PAIPhotonModel::ComputeDEDXPerVolume: jMat= " << jMat
836  << " jMax= " << fMaterialCutsCoupleVector.size()
837  << " matCC: " << matCC;
838  if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
839  G4cout << G4endl;
840  G4cout << fPAIdEdxTable << " " << fdEdxVector << " "
841  << fProtonEnergyVector << G4endl;
842  */
843  fPAIdEdxTable = fPAIdEdxBank[jMat];
844  fdEdxVector = fdEdxTable[jMat];
845  for(iTkin = 0; iTkin <= fTotBin; iTkin++)
846  {
847  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
848  }
849  iPlace = iTkin - 1;
850  if(iPlace < 0) iPlace = 0;
851  else if(iPlace > fTotBin) iPlace = fTotBin;
852  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
853 
854  if( dEdx < 0.) dEdx = 0.;
855  return dEdx;
856 }
857 
858 //////////////////////////////////////////////////////////////////////
859 //
860 // Return xsc for mean interaction length
861 
862 G4double G4PAIPhotonModel::CrossSectionPerVolume( const G4Material*,
863  const G4ParticleDefinition* p,
864  G4double kineticEnergy,
865  G4double cutEnergy,
866  G4double maxEnergy )
867 {
868  G4int iTkin,iPlace;
869  size_t jMat, jMatCC;
870  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
871  if(cutEnergy >= tmax) return 0.0;
872  G4double particleMass = p->GetPDGMass();
873  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
874  G4double charge = p->GetPDGCharge();
875  G4double charge2 = charge*charge, cross, cross1, cross2;
876  G4double photon1, photon2, plasmon1, plasmon2;
877 
878  const G4MaterialCutsCouple* matCC = CurrentCouple();
879 
880  const G4ProductionCutsTable* theCoupleTable=
881  G4ProductionCutsTable::GetProductionCutsTable();
882 
883  size_t numOfCouples = theCoupleTable->GetTableSize();
884 
885  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
886  {
887  if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
888  }
889  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
890 
891  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
892  GetEnergyCutsVector(idxG4GammaCut);
893 
894  G4double photonCut = (*photonCutInKineticEnergy)[jMatCC];
895 
896  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
897  {
898  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
899  }
900  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
901 
902  fPAItransferTable = fPAIxscBank[jMat];
903  fPAIphotonTable = fPAIphotonBank[jMat];
904  fPAIplasmonTable = fPAIplasmonBank[jMat];
905 
906  for(iTkin = 0; iTkin <= fTotBin; iTkin++)
907  {
908  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
909  }
910  iPlace = iTkin - 1;
911  if(iPlace < 0) iPlace = 0;
912 
913  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
914  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
915  photon1 = GetdNdxPhotonCut(iPlace,tmax);
916  photon2 = GetdNdxPhotonCut(iPlace,photonCut);
917 
918  plasmon1 = GetdNdxPlasmonCut(iPlace,tmax);
919  plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy);
920 
921  cross1 = photon1 + plasmon1;
922  // G4cout<<"cross1 = "<<cross1<<G4endl;
923  cross2 = photon2 + plasmon2;
924  // G4cout<<"cross2 = "<<cross2<<G4endl;
925  cross = (cross2 - cross1)*charge2;
926  // G4cout<<"cross = "<<cross<<G4endl;
927 
928  if( cross < 0. ) cross = 0.;
929  return cross;
930 }
931 //////////////////////////////////////////////////////////////////////
932 //
933 // Return xsc for mean interaction length in test
934 
935 G4double G4PAIPhotonModel::GetXscPerVolume( const G4Material*,
936  const G4ParticleDefinition* p,
937  G4double kineticEnergy,
938  G4double photonCut,
939  G4double cutEnergy,
940  G4double maxEnergy )
941 {
942  G4int iTkin,iPlace;
943  size_t jMat, jMatCC;
944  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
945  if(cutEnergy >= tmax) return 0.0;
946  G4double particleMass = p->GetPDGMass();
947  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
948  G4double charge = p->GetPDGCharge();
949  G4double charge2 = charge*charge, cross, cross1, cross2;
950  G4double photon1, photon2, plasmon1, plasmon2;
951 
952  const G4MaterialCutsCouple* matCC = CurrentCouple();
953 
954  const G4ProductionCutsTable* theCoupleTable=
955  G4ProductionCutsTable::GetProductionCutsTable();
956 
957  size_t numOfCouples = theCoupleTable->GetTableSize();
958 
959  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
960  {
961  if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
962  }
963  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
964 
965  // const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
966  // GetEnergyCutsVector(idxG4GammaCut);
967  // G4double photonCut = (*photonCutInKineticEnergy)[jMatCC];
968 
969  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
970  {
971  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
972  }
973  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
974 
975  fPAItransferTable = fPAIxscBank[jMat];
976  fPAIphotonTable = fPAIphotonBank[jMat];
977  fPAIplasmonTable = fPAIplasmonBank[jMat];
978 
979  for(iTkin = 0; iTkin <= fTotBin; iTkin++)
980  {
981  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
982  }
983  iPlace = iTkin - 1;
984  if(iPlace < 0) iPlace = 0;
985 
986  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
987  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
988  photon1 = GetdNdxPhotonCut(iPlace,tmax);
989  photon2 = GetdNdxPhotonCut(iPlace,photonCut);
990 
991  plasmon1 = GetdNdxPlasmonCut(iPlace,tmax);
992  plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy);
993 
994  cross1 = photon1 + plasmon1;
995  // G4cout<<"cross1 = "<<cross1<<G4endl;
996  cross2 = photon2 + plasmon2;
997  // G4cout<<"cross2 = "<<cross2<<G4endl;
998  cross = (cross2 - cross1)*charge2;
999  // G4cout<<"cross = "<<cross<<G4endl;
1000 
1001  if( cross < 0. ) cross = 0.;
1002  return cross;
1003 }
1004 
1005 ///////////////////////////////////////////////////////////////////////////
1006 //
1007 // It is analog of PostStepDoIt in terms of secondary electron or photon to
1008 // be returned as G4Dynamicparticle*.
1009 //
1010 
1011 void G4PAIPhotonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
1012  const G4MaterialCutsCouple* matCC,
1013  const G4DynamicParticle* dp,
1014  G4double tmin,
1015  G4double maxEnergy)
1016 {
1017  size_t jMat;
1018  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
1019  {
1020  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
1021  }
1022  if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--;
1023 
1024  fPAItransferTable = fPAIxscBank[jMat];
1025  fPAIphotonTable = fPAIphotonBank[jMat];
1026  fPAIplasmonTable = fPAIplasmonBank[jMat];
1027 
1028  fdNdxCutVector = fdNdxCutTable[jMat];
1029  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1030  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1031 
1032  G4double tmax = std::min(MaxSecondaryEnergy(dp->GetDefinition(),dp->GetKineticEnergy()), maxEnergy);
1033  if( tmin >= tmax && fVerbose > 0)
1034  {
1035  G4cout<<"G4PAIPhotonModel::SampleSecondary: tmin >= tmax "<<G4endl;
1036  }
1037 
1038  G4ThreeVector direction = dp->GetMomentumDirection();
1039  G4double particleMass = dp->GetMass();
1040  G4double kineticEnergy = dp->GetKineticEnergy();
1041  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; // fMass
1042  G4double totalEnergy = kineticEnergy + particleMass;
1043  G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
1044 
1045  G4int iTkin;
1046  for(iTkin=0;iTkin<=fTotBin;iTkin++)
1047  {
1048  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
1049  }
1050  G4int iPlace = iTkin - 1;
1051  if(iPlace < 0) iPlace = 0;
1052 
1053  G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace);
1054  G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace);
1055  G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
1056 
1057  G4double ratio;
1058  if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut;
1059  else return; // ratio = 0.;
1060 
1061  if(ratio < G4UniformRand() ) // secondary e-
1062  {
1063  G4double deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector,
1064  iPlace, scaledTkin);
1065 
1066 // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1067 
1068  if( deltaTkin <= 0. )
1069  {
1070  G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
1071  }
1072  if( deltaTkin <= 0.) return;
1073 
1074  if( deltaTkin >= kineticEnergy ) // stop primary
1075  {
1076  deltaTkin = kineticEnergy;
1077  kineticEnergy = 0.0;
1078  }
1079  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
1080  G4double totalMomentum = sqrt(pSquare);
1081  G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
1082  /(deltaTotalMomentum * totalMomentum);
1083 
1084  if( costheta > 0.99999 ) costheta = 0.99999;
1085  G4double sintheta = 0.0;
1086  G4double sin2 = 1. - costheta*costheta;
1087  if( sin2 > 0.) sintheta = sqrt(sin2);
1088 
1089  // direction of the delta electron
1090 
1091  G4double phi = twopi*G4UniformRand();
1092  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1093 
1094  G4ThreeVector deltaDirection(dirx,diry,dirz);
1095  deltaDirection.rotateUz(direction);
1096 
1097  if( kineticEnergy > 0.) // primary change
1098  {
1099  kineticEnergy -= deltaTkin;
1100  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
1101  direction = dir.unit();
1102  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
1103  fParticleChange->SetProposedMomentumDirection(direction);
1104  }
1105  else // stop primary
1106  {
1107  fParticleChange->ProposeTrackStatus(fStopAndKill);
1108  fParticleChange->SetProposedKineticEnergy(0.0);
1109  }
1110 
1111  // create G4DynamicParticle object for e- delta ray
1112 
1113  G4DynamicParticle* deltaRay = new G4DynamicParticle;
1114  deltaRay->SetDefinition(G4Electron::Electron());
1115  deltaRay->SetKineticEnergy( deltaTkin );
1116  deltaRay->SetMomentumDirection(deltaDirection);
1117  vdp->push_back(deltaRay);
1118 
1119  }
1120  else // secondary 'Cherenkov' photon
1121  {
1122  G4double deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector,
1123  iPlace,scaledTkin);
1124 
1125  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1126 
1127  if( deltaTkin <= 0. )
1128  {
1129  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
1130  }
1131  if( deltaTkin <= 0.) return;
1132 
1133  if( deltaTkin >= kineticEnergy ) // stop primary
1134  {
1135  deltaTkin = kineticEnergy;
1136  kineticEnergy = 0.0;
1137  }
1138  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
1139  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
1140 
1141  // direction of the 'Cherenkov' photon
1142  G4double phi = twopi*G4UniformRand();
1143  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1144 
1145  G4ThreeVector deltaDirection(dirx,diry,dirz);
1146  deltaDirection.rotateUz(direction);
1147 
1148  if( kineticEnergy > 0.) // primary change
1149  {
1150  kineticEnergy -= deltaTkin;
1151  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
1152  }
1153  else // stop primary
1154  {
1155  fParticleChange->ProposeTrackStatus(fStopAndKill);
1156  fParticleChange->SetProposedKineticEnergy(0.0);
1157  }
1158  // create G4DynamicParticle object for photon ray
1159 
1160  G4DynamicParticle* photonRay = new G4DynamicParticle;
1161  photonRay->SetDefinition( G4Gamma::Gamma() );
1162  photonRay->SetKineticEnergy( deltaTkin );
1163  photonRay->SetMomentumDirection(deltaDirection);
1164 
1165  vdp->push_back(photonRay);
1166  }
1167 }
1168 
1169 ///////////////////////////////////////////////////////////////////////////////
1170 //
1171 // test function for losses more than cut
1172 
1173 G4double G4PAIPhotonModel::TestSecondaries( G4MaterialCutsCouple* matCC, G4DynamicParticle* dp,
1174  G4double tmin,
1175  G4double maxEnergy)
1176 {
1177  size_t jMat;
1178  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
1179  {
1180  if( matCC->GetMaterial()->GetName() == fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName() ) break;
1181  }
1182  if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--;
1183 
1184  fPAItransferTable = fPAIxscBank[jMat];
1185  fPAIphotonTable = fPAIphotonBank[jMat];
1186  fPAIplasmonTable = fPAIplasmonBank[jMat];
1187 
1188  fdNdxCutVector = fdNdxCutTable[jMat];
1189  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1190  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1191 
1192  G4double tmax = std::min(MaxSecondaryEnergy(dp->GetDefinition(),dp->GetKineticEnergy()), maxEnergy);
1193  if( tmin >= tmax && fVerbose > 0)
1194  {
1195  G4cout<<"G4PAIPhotonModel::TestSecondaries: tmin >= tmax "<<G4endl;
1196  }
1197 
1198  G4ThreeVector direction = dp->GetMomentumDirection();
1199  G4double particleMass = dp->GetMass();
1200  G4double kineticEnergy = dp->GetKineticEnergy();
1201  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; // fMass
1202  G4double totalEnergy = kineticEnergy + particleMass;
1203  G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
1204 
1205  G4int iTkin;
1206  for(iTkin=0;iTkin<=fTotBin;iTkin++)
1207  {
1208  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
1209  }
1210  G4int iPlace = iTkin - 1;
1211  if(iPlace < 0) iPlace = 0;
1212 
1213  G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace);
1214  G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace);
1215  G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
1216 
1217  G4double ratio, deltaTkin;
1218  if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut;
1219  else ratio = 0.;
1220 
1221  if(ratio < G4UniformRand() ) // secondary e-
1222  {
1223  deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector,
1224  iPlace, scaledTkin);
1225 
1226 // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1227 
1228  if( deltaTkin <= 0. )
1229  {
1230  G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
1231  }
1232  if( deltaTkin <= 0.) return 0.;
1233 
1234  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
1235  G4double totalMomentum = sqrt(pSquare);
1236  G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
1237  /(deltaTotalMomentum * totalMomentum);
1238 
1239  if( costheta > 0.99999 ) costheta = 0.99999;
1240  G4double sintheta = 0.0;
1241  G4double sin2 = 1. - costheta*costheta;
1242  if( sin2 > 0.) sintheta = sqrt(sin2);
1243 
1244  // direction of the delta electron
1245 
1246  G4double phi = twopi*G4UniformRand();
1247  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1248 
1249  G4ThreeVector deltaDirection(dirx,diry,dirz);
1250  deltaDirection.rotateUz(direction);
1251 
1252  // primary change
1253 
1254  kineticEnergy -= deltaTkin;
1255  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
1256  direction = dir.unit();
1257  fParticleChange->SetProposedMomentumDirection(direction);
1258 
1259  // create G4DynamicParticle object for e- delta ray
1260 
1261  G4DynamicParticle* deltaRay = new G4DynamicParticle;
1262  deltaRay->SetDefinition(G4Electron::Electron());
1263  deltaRay->SetKineticEnergy( deltaTkin );
1264  deltaRay->SetMomentumDirection(deltaDirection);
1265 
1266  }
1267  else // secondary 'Cherenkov' photon
1268  {
1269  deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector,
1270  iPlace,scaledTkin);
1271 
1272  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1273 
1274  if( deltaTkin <= 0. )
1275  {
1276  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
1277  }
1278  if( deltaTkin <= 0.) return 0.;
1279 
1280  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
1281  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
1282 
1283  // direction of the 'Cherenkov' photon
1284  G4double phi = twopi*G4UniformRand();
1285  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1286 
1287  G4ThreeVector deltaDirection(dirx,diry,dirz);
1288  deltaDirection.rotateUz(direction);
1289 
1290  // primary change
1291  kineticEnergy -= deltaTkin;
1292 
1293  // create G4DynamicParticle object for photon ray
1294 
1295  G4DynamicParticle* photonRay = new G4DynamicParticle;
1296  photonRay->SetDefinition( G4Gamma::Gamma() );
1297  photonRay->SetKineticEnergy( deltaTkin );
1298  photonRay->SetMomentumDirection(deltaDirection);
1299 
1300  }
1301  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
1302 
1303  return deltaTkin;
1304 }
1305 
1306 
1307 ///////////////////////////////////////////////////////////////////////
1308 //
1309 // Returns post step PAI energy transfer > cut electron/photon energy according to passed
1310 // scaled kinetic energy of particle
1311 
1312 G4double
1313 G4PAIPhotonModel::GetPostStepTransfer( G4PhysicsTable* pTable,
1314  G4PhysicsLogVector* pVector,
1315  G4int iPlace, G4double scaledTkin )
1316 {
1317  // G4cout<<"G4PAIPhotonModel::GetPostStepTransfer"<<G4endl;
1318 
1319  G4int iTkin = iPlace+1, iTransfer;
1320  G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
1321 
1322  dNdxCut1 = (*pVector)(iPlace);
1323 
1324  // G4cout<<"iPlace = "<<iPlace<<endl;
1325 
1326  if(iTkin == fTotBin) // Fermi plato, try from left
1327  {
1328  position = dNdxCut1*G4UniformRand();
1329 
1330  for( iTransfer = 0;
1331  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1332  {
1333  if(position >= (*(*pTable)(iPlace))(iTransfer)) break;
1334  }
1335  transfer = GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1336  }
1337  else
1338  {
1339  dNdxCut2 = (*pVector)(iPlace+1);
1340  if(iTkin == 0) // Tkin is too small, trying from right only
1341  {
1342  position = dNdxCut2*G4UniformRand();
1343 
1344  for( iTransfer = 0;
1345  iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
1346  {
1347  if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break;
1348  }
1349  transfer = GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
1350  }
1351  else // general case: Tkin between two vectors of the material
1352  {
1353  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
1354  E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1355  W = 1.0/(E2 - E1);
1356  W1 = (E2 - scaledTkin)*W;
1357  W2 = (scaledTkin - E1)*W;
1358 
1359  position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand();
1360 
1361  // G4cout<<position<<"\t";
1362 
1363  G4int iTrMax1, iTrMax2, iTrMax;
1364 
1365  iTrMax1 = G4int((*pTable)(iPlace)->GetVectorLength());
1366  iTrMax2 = G4int((*pTable)(iPlace+1)->GetVectorLength());
1367 
1368  if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
1369  else iTrMax = iTrMax1;
1370 
1371  for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
1372  {
1373  if( position >=
1374  ( (*(*pTable)(iPlace))(iTransfer)*W1 +
1375  (*(*pTable)(iPlace+1))(iTransfer)*W2) ) break;
1376  }
1377  transfer = GetEnergyTransfer(pTable, iPlace, position, iTransfer);
1378  }
1379  }
1380  // G4cout<<"PAIPhotonModel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
1381  if( transfer < 0.0 ) transfer = 0.0;
1382  return transfer;
1383 }
1384 
1385 ///////////////////////////////////////////////////////////////////////
1386 //
1387 // Returns random PAI energy transfer according to passed
1388 // indexes of particle
1389 
1390 G4double
1391 G4PAIPhotonModel::GetEnergyTransfer( G4PhysicsTable* pTable, G4int iPlace,
1392  G4double position, G4int iTransfer )
1393 {
1394  G4int iTransferMax;
1395  G4double x1, x2, y1, y2, energyTransfer;
1396 
1397  if(iTransfer == 0)
1398  {
1399  energyTransfer = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1400  }
1401  else
1402  {
1403  iTransferMax = G4int((*pTable)(iPlace)->GetVectorLength());
1404 
1405  if ( iTransfer >= iTransferMax) iTransfer = iTransferMax - 1;
1406 
1407  y1 = (*(*pTable)(iPlace))(iTransfer-1);
1408  y2 = (*(*pTable)(iPlace))(iTransfer);
1409 
1410  x1 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1411  x2 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1412 
1413  if ( x1 == x2 ) energyTransfer = x2;
1414  else
1415  {
1416  if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
1417  else
1418  {
1419  energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1420  }
1421  }
1422  }
1423  return energyTransfer;
1424 }
1425 
1426 ///////////////////////////////////////////////////////////////////////
1427 //
1428 // Works like AlongStepDoIt method of process family
1429 
1430 G4double
1431 G4PAIPhotonModel::SampleFluctuations( const G4MaterialCutsCouple* matCC,
1432  const G4DynamicParticle* aParticle,
1433  G4double,
1434  G4double step,
1435  G4double eloss)
1436 {
1437  size_t jMat = 0;
1438  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
1439  {
1440  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
1441  }
1442  if(jMat == fMaterialCutsCoupleVector.size()) { return eloss; }
1443 
1444  fPAItransferTable = fPAIxscBank[jMat];
1445  fPAIphotonTable = fPAIphotonBank[jMat];
1446  fPAIplasmonTable = fPAIplasmonBank[jMat];
1447 
1448  fdNdxCutVector = fdNdxCutTable[jMat];
1449  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1450  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1451 
1452  G4int iTkin, iPlace ;
1453 
1454  // G4cout<<"G4PAIPhotonModel::SampleFluctuations"<<G4endl;
1455 
1456  G4double loss, photonLoss, plasmonLoss, charge2;
1457 
1458 
1459  G4double Tkin = aParticle->GetKineticEnergy();
1460  G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass();
1461  G4double charge = aParticle->GetDefinition()->GetPDGCharge();
1462  charge2 = charge*charge;
1463  G4double scaledTkin = Tkin*MassRatio;
1464  G4double cof = step*charge2;
1465 
1466  for( iTkin = 0; iTkin <= fTotBin; iTkin++)
1467  {
1468  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
1469  }
1470  iPlace = iTkin - 1;
1471  if( iPlace < 0 ) iPlace = 0;
1472 
1473  photonLoss = GetAlongStepTransfer(fPAIphotonTable,fdNdxCutPhotonVector,
1474 iPlace,scaledTkin,step,cof);
1475 
1476  // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl;
1477 
1478  plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector,
1479 iPlace,scaledTkin,step,cof);
1480 
1481  // G4cout<<"PAIPhotonModel AlongStepPlasmonLoss = "<<plasmonLoss/keV<<" keV"<<G4endl;
1482 
1483  loss = photonLoss + plasmonLoss;
1484 
1485  // G4cout<<"PAIPhotonModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl;
1486 
1487  return loss;
1488 }
1489 
1490 ///////////////////////////////////////////////////////////////////////
1491 //
1492 // Returns along step PAI energy transfer < cut electron/photon energy according to passed
1493 // scaled kinetic energy of particle and cof = step*charge*charge
1494 
1495 G4double
1496 G4PAIPhotonModel::GetAlongStepTransfer( G4PhysicsTable* pTable,
1497  G4PhysicsLogVector* pVector,
1498  G4int iPlace, G4double scaledTkin,G4double step,
1499  G4double cof )
1500 {
1501  G4int iTkin = iPlace + 1, iTransfer;
1502  G4double loss = 0., position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
1503  G4double lambda, stepDelta, stepSum=0.;
1504  G4long numOfCollisions=0;
1505  G4bool numb = true;
1506 
1507  dNdxCut1 = (*pVector)(iPlace);
1508 
1509  // G4cout<<"iPlace = "<<iPlace<<endl;
1510 
1511  if(iTkin == fTotBin) // Fermi plato, try from left
1512  {
1513  meanNumber = ((*(*pTable)(iPlace))(0) - dNdxCut1)*cof;
1514  if(meanNumber < 0.) meanNumber = 0.;
1515  // numOfCollisions = RandPoisson::shoot(meanNumber);
1516  if( meanNumber > 0.) lambda = step/meanNumber;
1517  else lambda = DBL_MAX;
1518  while(numb)
1519  {
1520  stepDelta = G4RandExponential::shoot(lambda);
1521  stepSum += stepDelta;
1522  if(stepSum >= step) break;
1523  numOfCollisions++;
1524  }
1525 
1526  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1527 
1528  while(numOfCollisions)
1529  {
1530  position = dNdxCut1+
1531  ((*(*pTable)(iPlace))(0) - dNdxCut1)*G4UniformRand();
1532 
1533  for( iTransfer = 0;
1534  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1535  {
1536  if(position >= (*(*pTable)(iPlace))(iTransfer)) break;
1537  }
1538  loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1539  numOfCollisions--;
1540  }
1541  }
1542  else
1543  {
1544  dNdxCut2 = (*pVector)(iPlace+1);
1545 
1546  if(iTkin == 0) // Tkin is too small, trying from right only
1547  {
1548  meanNumber = ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*cof;
1549  if( meanNumber < 0. ) meanNumber = 0.;
1550  // numOfCollisions = G4RandPoisson::shoot(meanNumber);
1551  if( meanNumber > 0.) lambda = step/meanNumber;
1552  else lambda = DBL_MAX;
1553  while(numb)
1554  {
1555  stepDelta = G4RandExponential::shoot(lambda);
1556  stepSum += stepDelta;
1557  if(stepSum >= step) break;
1558  numOfCollisions++;
1559  }
1560 
1561  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1562 
1563  while(numOfCollisions)
1564  {
1565  position = dNdxCut2+
1566  ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*G4UniformRand();
1567 
1568  for( iTransfer = 0;
1569  iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
1570  {
1571  if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break;
1572  }
1573  loss += GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
1574  numOfCollisions--;
1575  }
1576  }
1577  else // general case: Tkin between two vectors of the material
1578  {
1579  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
1580  E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1581  W = 1.0/(E2 - E1);
1582  W1 = (E2 - scaledTkin)*W;
1583  W2 = (scaledTkin - E1)*W;
1584 
1585  // G4cout<<"(*(*pTable)(iPlace))(0) = "<<
1586  // (*(*pTable)(iPlace))(0)<<G4endl;
1587  // G4cout<<"(*(*pTable)(iPlace+1))(0) = "<<
1588  // (*(*pTable)(iPlace+1))(0)<<G4endl;
1589 
1590  meanNumber=( ((*(*pTable)(iPlace))(0)-dNdxCut1)*W1 +
1591  ((*(*pTable)(iPlace+1))(0)-dNdxCut2)*W2 )*cof;
1592  if(meanNumber<0.0) meanNumber = 0.0;
1593  // numOfCollisions = G4RandPoisson::shoot(meanNumber);
1594  if( meanNumber > 0.) lambda = step/meanNumber;
1595  else lambda = DBL_MAX;
1596  while(numb)
1597  {
1598  stepDelta = G4RandExponential::shoot(lambda);
1599  stepSum += stepDelta;
1600  if(stepSum >= step) break;
1601  numOfCollisions++;
1602  }
1603 
1604  // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl;
1605 
1606  while(numOfCollisions)
1607  {
1608  position = dNdxCut1*W1 + dNdxCut2*W2 +
1609  ( ( (*(*pTable)(iPlace ))(0) - dNdxCut1)*W1 +
1610 
1611  ( (*(*pTable)(iPlace+1))(0) - dNdxCut2)*W2 )*G4UniformRand();
1612 
1613  // G4cout<<position<<"\t";
1614 
1615  for( iTransfer = 0;
1616  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1617  {
1618  if( position >=
1619  ( (*(*pTable)(iPlace))(iTransfer)*W1 +
1620  (*(*pTable)(iPlace+1))(iTransfer)*W2) )
1621  {
1622  break;
1623  }
1624  }
1625  // loss += (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1626  loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1627  numOfCollisions--;
1628  }
1629  }
1630  }
1631 
1632  return loss;
1633 
1634 }
1635 
1636 //////////////////////////////////////////////////////////////////////
1637 //
1638 // Returns the statistical estimation of the energy loss distribution variance
1639 //
1640 
1641 
1642 G4double G4PAIPhotonModel::Dispersion( const G4Material* material,
1643  const G4DynamicParticle* aParticle,
1644  G4double tmax,
1645  G4double step )
1646 {
1647  G4double particleMass = aParticle->GetMass();
1648  G4double electronDensity = material->GetElectronDensity();
1649  G4double kineticEnergy = aParticle->GetKineticEnergy();
1650  G4double q = aParticle->GetCharge()/eplus;
1651  G4double etot = kineticEnergy + particleMass;
1652  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
1653  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
1654  * electronDensity * q * q;
1655 
1656  return siga;
1657 
1658  /*
1659  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
1660  for(G4int i = 0; i < fMeanNumber; i++)
1661  {
1662  loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
1663  sumLoss += loss;
1664  sumLoss2 += loss*loss;
1665  }
1666  meanLoss = sumLoss/fMeanNumber;
1667  sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
1668  return sigma2;
1669  */
1670 }
1671 
1672 /////////////////////////////////////////////////////////////////////
1673 
1674 G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
1675  G4double kinEnergy)
1676 {
1677  G4double tmax = kinEnergy;
1678  if(p == fElectron) tmax *= 0.5;
1679  else if(p != fPositron)
1680  {
1681  G4double mass = p->GetPDGMass();
1682  G4double ratio= electron_mass_c2/mass;
1683  G4double gamma= kinEnergy/mass + 1.0;
1684  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
1685  (1. + 2.0*gamma*ratio + ratio*ratio);
1686  }
1687  return tmax;
1688 }
1689 
1690 ///////////////////////////////////////////////////////////////
1691 
1692 void G4PAIPhotonModel::DefineForRegion(const G4Region* r)
1693 {
1694  fPAIRegionVector.push_back(r);
1695 }
1696 
1697 
1698 //
1699 //
1700 /////////////////////////////////////////////////
1701 
1702 
1703 
1704 
1705 
1706 
G4ParticleChangeForLoss.hh
G4PAIPhotonModel::TestSecondaries
G4double TestSecondaries(G4MaterialCutsCouple *, G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4PAIPhotonModel.cc:1173
G4PAIPhotonModel::GetdEdxCut
G4double GetdEdxCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:769
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250
G4VEmModel::GetParticleChangeForLoss
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:107
G4PhysicsFreeVector::PutValue
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Definition: G4PhysicsFreeVector.cc:94
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ProductionCuts.hh
G4PAIPhotonModel::GetXscPerVolume
G4double GetXscPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double photonCut, G4double cutEnergy, G4double maxEnergy)
Definition: G4PAIPhotonModel.cc:935
G4PAIPhotonModel::DefineForRegion
void DefineForRegion(const G4Region *r)
Definition: G4PAIPhotonModel.cc:1692
G4MaterialCutsCouple.hh
p
const char * p
Definition: xmltok.h:285
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:176
G4PAIPhotonModel::GetEnergyTransfer
G4double GetEnergyTransfer(G4PhysicsTable *, G4int iPlace, G4double position, G4int iTransfer)
Definition: G4PAIPhotonModel.cc:1391
G4PAIxSection::GetIntegralPlasmon
G4double GetIntegralPlasmon(G4int i) const
Definition: G4PAIxSection.hh:336
G4PAIPhotonModel::ComputeDEDXPerVolume
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
Definition: G4PAIPhotonModel.cc:811
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:87
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4long
long G4long
Definition: G4Types.hh:80
G4DynamicParticle.hh
G4Step.hh
G4PAIxSection::GetIntegralPAIxSection
G4double GetIntegralPAIxSection(G4int i) const
Definition: G4PAIxSection.hh:312
idxG4GammaCut
Definition: G4ProductionCuts.hh:53
G4PAIPhotonModel::GetdNdxPlasmonCut
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:723
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int)
Definition: G4SandiaTable.cc:986
G4PhysicsLogVector.hh
G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4PAIPhotonModel.hh
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4PAIxSection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition: G4PAIxSection.hh:318
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:157
G4int
int G4int
Definition: G4Types.hh:78
G4PAIPhotonModel::SampleSecondaries
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Definition: G4PAIPhotonModel.cc:1011
idxG4ElectronCut
Definition: G4ProductionCuts.hh:54
G4SandiaTable::SandiaMixing
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
Definition: G4SandiaTable.cc:810
G4Region.hh
G4PAIPhotonModel::BuildLambdaVector
void BuildLambdaVector(const G4MaterialCutsCouple *matCutsCouple)
Definition: G4PAIPhotonModel.cc:470
eplot.material
string material
Definition: eplot.py:19
G4PAIxSection.hh
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:426
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4PAIPhotonModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Definition: G4PAIPhotonModel.cc:862
G4PhysicsTable.hh
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4PAIPhotonModel::InitTest
virtual void InitTest(const G4ParticleDefinition *, G4MaterialCutsCouple *, G4double, G4double)
Definition: G4PAIPhotonModel.cc:238
G4DynamicParticle::GetMass
G4double GetMass() const
G4PAIxSection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIxSection.hh:187
G4bool
bool G4bool
Definition: G4Types.hh:79
G4PAIPhotonModel::GetdNdxPhotonCut
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:678
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4Poisson.hh
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1674
G4Material.hh
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
G4DynamicParticle::GetCharge
G4double GetCharge() const
G4PAIxSection::GetIntegralCerenkov
G4double GetIntegralCerenkov(G4int i) const
Definition: G4PAIxSection.hh:324
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4ParticleChangeForLoss::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:150
G4Positron.hh
G4EmCalculator.hh
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4ParticleDefinition.hh
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
Randomize.hh
G4Region::FindCouple
G4MaterialCutsCouple * FindCouple(G4Material *mat)
G4PAIPhotonModel::SampleFluctuations
virtual G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double, G4double)
Definition: G4PAIPhotonModel.cc:1431
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:571
G4ParticleChangeForLoss::SetProposedMomentumDirection
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:194
G4PAIPhotonModel::GetAlongStepTransfer
G4double GetAlongStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin, G4double step, G4double cof)
Definition: G4PAIPhotonModel.cc:1496
G4DynamicParticle::SetKineticEnergy
void SetKineticEnergy(G4double aEnergy)
G4SandiaTable::SandiaIntervals
G4int SandiaIntervals(G4int Z[], G4int el)
Definition: G4SandiaTable.cc:714
G4PhysicalConstants.hh
G4PAIPhotonModel::Dispersion
virtual G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double, G4double)
Definition: G4PAIPhotonModel.cc:1642
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron.hh
G4ProductionCutsTable.hh
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4PAIPhotonModel::~G4PAIPhotonModel
virtual ~G4PAIPhotonModel()
Definition: G4PAIPhotonModel.cc:129
G4PAIPhotonModel::Initialise
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4PAIPhotonModel.cc:165
G4PAIxSection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIxSection.cc:606
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
position
int position
Definition: filter.cc:7
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
fStopAndKill
Definition: G4TrackStatus.hh:56
G4SandiaTable::GetPhotoAbsorpCof
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
Definition: G4SandiaTable.cc:687
G4PAIPhotonModel::GetdNdxCut
G4double GetdNdxCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:633
G4PhysicsTable::insertAt
void insertAt(size_t, G4PhysicsVector *)
G4ParticleDefinition::GetPDGSpin
G4double GetPDGSpin() const
Definition: G4ParticleDefinition.hh:165
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4PhysicsFreeVector.hh
G4MaterialTable.hh
G4SandiaTable.hh
G4PAIPhotonModel::G4PAIPhotonModel
G4PAIPhotonModel(const G4ParticleDefinition *p=0, const G4String &nam="PAIPhoton")
Definition: G4PAIPhotonModel.cc:80
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:273
G4SystemOfUnits.hh
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4GeometryTolerance.hh
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4PAIPhotonModel::GetPostStepTransfer
G4double GetPostStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin)
Definition: G4PAIPhotonModel.cc:1313
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:74
G4PAIxSection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition: G4PAIxSection.hh:199
G4SandiaTable::Initialize
void Initialize(G4Material *)
Definition: G4SandiaTable.cc:625
G4PAIxSection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition: G4PAIxSection.hh:306
G4Gamma.hh
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
Generated on Wed Apr 30 2014 15:54:58 for Geant4.10 by   doxygen 1.8.7