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