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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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