source: trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc @ 924

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

import all except CVS
File size: 30.9 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:     G4PAIModel.cc
27//
28// Author: Vladimir.Grichine@cern.ch on base of Vladimir Ivanchenko code
29//
30// Creation date: 05.10.2003
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// 08.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 "G4OrderedTable.hh"
48
49#include "G4PAIModel.hh"
50#include "Randomize.hh"
51#include "G4Electron.hh"
52#include "G4Positron.hh"
53#include "G4Poisson.hh"
54#include "G4Step.hh"
55#include "G4Material.hh"
56#include "G4DynamicParticle.hh"
57#include "G4ParticleDefinition.hh"
58#include "G4ParticleChangeForLoss.hh"
59#include "G4GeometryTolerance.hh"
60
61////////////////////////////////////////////////////////////////////////
62
63using namespace std;
64
65G4PAIModel::G4PAIModel(const G4ParticleDefinition* p, const G4String& nam)
66  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
67  fVerbose(0),
68  fLowestGamma(1.005),
69  fHighestGamma(10000.),
70  fTotBin(200),
71  fMeanNumber(20),
72  fParticle(0),
73  fHighKinEnergy(100.*TeV),
74  fTwoln10(2.0*log(10.0)),
75  fBg2lim(0.0169),
76  fTaulim(8.4146e-3)
77{
78  if(p) SetParticle(p);
79 
80  fElectron = G4Electron::Electron();
81  fPositron = G4Positron::Positron();
82
83  fPAItransferTable  = 0;
84  fPAIdEdxTable      = 0;
85  fSandiaPhotoAbsCof = 0;
86  fdEdxVector        = 0;
87  fLambdaVector      = 0;
88  fdNdxCutVector     = 0;
89
90  isInitialised      = false;
91}
92
93////////////////////////////////////////////////////////////////////////////
94
95G4PAIModel::~G4PAIModel()
96{
97  //  G4cout << "PAI: start destruction" << G4endl;
98  if(fParticleEnergyVector) delete fParticleEnergyVector;
99  if(fdEdxVector)           delete fdEdxVector ;
100  if(fLambdaVector)         delete fLambdaVector;
101  if(fdNdxCutVector)        delete fdNdxCutVector;
102
103  if( fPAItransferTable )
104    {
105      fPAItransferTable->clearAndDestroy();
106      delete fPAItransferTable ;
107    }
108  if( fPAIdEdxTable )
109    {
110      fPAIdEdxTable->clearAndDestroy();
111      delete fPAIdEdxTable ;
112    }
113  if(fSandiaPhotoAbsCof)
114    {
115      for(G4int i=0;i<fSandiaIntervalNumber;i++)
116        {
117          delete[] fSandiaPhotoAbsCof[i];
118        }
119      delete[] fSandiaPhotoAbsCof;
120    }
121  //G4cout << "PAI: end destruction" << G4endl;
122}
123
124///////////////////////////////////////////////////////////////////////////////
125
126void G4PAIModel::SetParticle(const G4ParticleDefinition* p)
127{
128  if(fParticle == p) return;
129  fParticle = p;
130  fMass = fParticle->GetPDGMass();
131  fSpin = fParticle->GetPDGSpin();
132  G4double q = fParticle->GetPDGCharge()/eplus;
133  fChargeSquare = q*q;
134  fLowKinEnergy = 0.2*MeV*fMass/proton_mass_c2;
135  fRatio = electron_mass_c2/fMass;
136  fQc = fMass/fRatio;
137}
138
139////////////////////////////////////////////////////////////////////////////
140
141void G4PAIModel::Initialise(const G4ParticleDefinition* p,
142                            const G4DataVector&)
143{
144  if(isInitialised) return;
145  isInitialised = true;
146
147  SetParticle(p);
148  fLowestKineticEnergy  = fMass*(fLowestGamma  - 1.0);
149  fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
150
151  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
152                                                 fHighestKineticEnergy,
153                                                 fTotBin);
154
155  if(pParticleChange)
156    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
157  else
158    fParticleChange = new G4ParticleChangeForLoss();
159
160  // Prepare initialization
161
162  fPAItransferTable = new G4PhysicsTable(fTotBin);
163  fPAIdEdxTable = new G4PhysicsTable(fTotBin);
164
165  const G4ProductionCutsTable* theCoupleTable =
166        G4ProductionCutsTable::GetProductionCutsTable();
167  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
168  size_t numOfMat   = G4Material::GetNumberOfMaterials();
169  size_t numRegions = fPAIRegionVector.size();
170
171  for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
172  {
173    const G4Region* curReg = fPAIRegionVector[iReg];
174    for(size_t jMat = 0; jMat < numOfMat; ++jMat) // region material loop
175    {
176      fMaterial  = (*theMaterialTable)[jMat];
177      fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial, 
178                                          curReg->GetProductionCuts() );
179      if( fCutCouple ) {
180        fMaterialCutsCoupleVector.push_back(fCutCouple);
181
182        fDeltaCutInKinEnergy = 
183          (*theCoupleTable->GetEnergyCutsVector(1))[fCutCouple->GetIndex()];
184     
185        //ComputeSandiaPhotoAbsCof();
186        BuildPAIonisationTable();
187
188        fPAIxscBank.push_back(fPAItransferTable);
189        fPAIdEdxBank.push_back(fPAIdEdxTable);
190        fdEdxTable.push_back(fdEdxVector);
191
192        BuildLambdaVector();
193        fdNdxCutTable.push_back(fdNdxCutVector);
194        fLambdaTable.push_back(fLambdaVector);
195      }
196    }
197  }
198}
199
200//////////////////////////////////////////////////////////////////
201
202void G4PAIModel::ComputeSandiaPhotoAbsCof()
203{ 
204  G4int i, j, numberOfElements ;
205  static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
206
207  G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
208  numberOfElements = (*theMaterialTable)[fMatIndex]->
209                                              GetNumberOfElements();
210  G4int* thisMaterialZ = new G4int[numberOfElements] ;
211
212  for(i=0;i<numberOfElements;i++) 
213  {
214    thisMaterialZ[i] = 
215    (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
216  } 
217  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
218                           (thisMaterialZ,numberOfElements) ;
219
220  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
221                           ( thisMaterialZ ,
222                             (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
223                             numberOfElements,fSandiaIntervalNumber) ;
224   
225  fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
226
227  for(i=0;i<fSandiaIntervalNumber;i++)  fSandiaPhotoAbsCof[i] = new G4double[5] ;
228   
229  for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
230  {
231    fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ; 
232
233    for( j = 1; j < 5 ; j++ )
234    {
235      fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
236                                      GetPhotoAbsorpCof(i+1,j)*
237                 (*theMaterialTable)[fMatIndex]->GetDensity() ;
238    }
239  }
240  // delete[] thisMaterialZ ;
241}
242
243////////////////////////////////////////////////////////////////////////////
244//
245// Build tables for the ionization energy loss
246//  the tables are built for MATERIALS
247//                           *********
248
249void G4PAIModel::BuildPAIonisationTable()
250{
251  G4double LowEdgeEnergy , ionloss ;
252  G4double tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ;
253
254  if(fdEdxVector) delete fdEdxVector;
255  fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
256                                        fHighestKineticEnergy,
257                                        fTotBin);
258  G4SandiaTable* sandia = fMaterial->GetSandiaTable();
259
260  Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV;
261  deltaLow = 100.*eV; // 0.5*eV ;
262
263  for (G4int i = 0 ; i < fTotBin ; i++)  //The loop for the kinetic energy
264  {
265    LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
266    tau = LowEdgeEnergy/fMass ;
267    gamma = tau +1. ;
268    // G4cout<<"gamma = "<<gamma<<endl ;
269    bg2 = tau*( tau + 2. );
270    Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy); 
271    //    Tmax = std::min(fDeltaCutInKinEnergy, Tmax);
272    Tkin = Tmax ;
273
274    // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
275    // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
276 
277    if ( Tmax < Tmin + deltaLow )  // low energy safety
278      Tkin = Tmin + deltaLow ;
279
280    /*
281    G4PAIxSection protonPAI( fMatIndex,
282                             Tkin,
283                             bg2,
284                             fSandiaPhotoAbsCof,
285                             fSandiaIntervalNumber  ) ;
286    */
287    fPAIySection.Initialize(fMaterial, Tkin, bg2);
288
289    // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
290    // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
291    // G4cout<<"protonPAI.GetSplineSize() = "<<
292    //    protonPAI.GetSplineSize()<<G4endl<<G4endl ;
293
294    G4int n = fPAIySection.GetSplineSize();
295    G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n) ;
296    G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
297
298    for( G4int k = 0 ; k < n; k++ )
299    {
300      transferVector->PutValue( k ,
301                                fPAIySection.GetSplineEnergy(k+1),
302                                fPAIySection.GetIntegralPAIySection(k+1) ) ;
303      dEdxVector->PutValue( k ,
304                                fPAIySection.GetSplineEnergy(k+1),
305                                fPAIySection.GetIntegralPAIdEdx(k+1) ) ;
306    }
307    ionloss = fPAIySection.GetMeanEnergyLoss() ;   //  total <dE/dx>
308
309    if ( ionloss < DBL_MIN)  ionloss = DBL_MIN;
310    fdEdxVector->PutValue(i,ionloss) ;
311
312    fPAItransferTable->insertAt(i,transferVector) ;
313    fPAIdEdxTable->insertAt(i,dEdxVector) ;
314
315  }                                        // end of Tkin loop
316  //  theLossTable->insert(fdEdxVector);
317  // end of material loop
318  // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;
319  // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;
320}
321
322///////////////////////////////////////////////////////////////////////
323//
324// Build mean free path tables for the delta ray production process
325//     tables are built for MATERIALS
326//
327
328void G4PAIModel::BuildLambdaVector()
329{
330  //G4double kCarTolerance = G4GeometryTolerance::GetInstance()
331  //                         ->GetSurfaceTolerance();
332
333  if (fLambdaVector)   delete fLambdaVector;
334  if (fdNdxCutVector)  delete fdNdxCutVector;
335
336  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
337                                          fHighestKineticEnergy,
338                                          fTotBin                ) ;
339  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
340                                          fHighestKineticEnergy,
341                                          fTotBin                ) ;
342  if(fVerbose > 1)
343  {
344    G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
345          <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
346  }
347  for (G4int i = 0 ; i < fTotBin ; i++ )
348  {
349    G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
350    G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
351
352    //    if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance ; // Mmm ???
353
354    fLambdaVector->PutValue(i, lambda) ;
355    fdNdxCutVector->PutValue(i, dNdxCut) ;
356  }
357}
358
359///////////////////////////////////////////////////////////////////////
360//
361// Returns integral PAI cross section for energy transfers >= transferCut
362
363G4double 
364G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
365{ 
366  G4int iTransfer;
367  G4double x1, x2, y1, y2, dNdxCut;
368  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
369  // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
370  //           <<G4endl; 
371  for( iTransfer = 0 ; 
372       iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ; 
373       iTransfer++)
374  {
375    if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
376    {
377      break ;
378    }
379  } 
380  if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
381  {
382      iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
383  }
384  y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
385  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
386  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
387  x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
388  x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
389  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
390
391  if ( y1 == y2 )    dNdxCut = y2 ;
392  else
393  {
394    //  if ( x1 == x2  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
395    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
396    else             dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;     
397  }
398  //  G4cout<<""<<dNdxCut<<G4endl;
399  return dNdxCut ;
400}
401///////////////////////////////////////////////////////////////////////
402//
403// Returns integral dEdx for energy transfers >= transferCut
404
405G4double 
406G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
407{ 
408  G4int iTransfer;
409  G4double x1, x2, y1, y2, dEdxCut;
410  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
411  // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
412  //           <<G4endl; 
413  for( iTransfer = 0 ; 
414       iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ; 
415       iTransfer++)
416  {
417    if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
418    {
419      break ;
420    }
421  } 
422  if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
423  {
424      iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
425  }
426  y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
427  y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
428  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
429  x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
430  x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
431  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
432
433  if ( y1 == y2 )    dEdxCut = y2 ;
434  else
435  {
436    //  if ( x1 == x2  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
437    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
438    else             dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;     
439  }
440  //  G4cout<<""<<dEdxCut<<G4endl;
441  return dEdxCut ;
442}
443
444//////////////////////////////////////////////////////////////////////////////
445
446G4double G4PAIModel::ComputeDEDX(const G4MaterialCutsCouple* matCC,
447                                 const G4ParticleDefinition* p,
448                                       G4double kineticEnergy,
449                                       G4double cutEnergy)
450{
451  G4int iTkin,iPlace;
452  size_t jMat;
453  G4double massRatio  = fMass/p->GetPDGMass();
454  G4double scaledTkin = kineticEnergy*massRatio;
455  G4double charge     = p->GetPDGCharge();
456  G4double charge2    = charge*charge, dEdx;
457
458  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
459  {
460    if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
461  }
462  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
463
464  fPAIdEdxTable = fPAIdEdxBank[jMat];
465  fdEdxVector = fdEdxTable[jMat];
466  for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
467  {
468    if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;   
469  }
470  iPlace = iTkin - 1;
471  if(iPlace < 0) iPlace = 0;
472  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cutEnergy) ) ; 
473
474  if( dEdx < 0.) dEdx = 0.;
475  return dEdx;
476}
477
478/////////////////////////////////////////////////////////////////////////
479
480G4double G4PAIModel::CrossSection( const G4MaterialCutsCouple* matCC,
481                                   const G4ParticleDefinition* p,
482                                         G4double kineticEnergy,
483                                         G4double cutEnergy,
484                                         G4double maxEnergy  ) 
485{
486  G4int iTkin,iPlace;
487  size_t jMat;
488  G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
489  G4double massRatio  = fMass/p->GetPDGMass();
490  G4double scaledTkin = kineticEnergy*massRatio;
491  G4double charge     = p->GetPDGCharge();
492  G4double charge2    = charge*charge, cross, cross1, cross2;
493
494  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
495  {
496    if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
497  }
498  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
499
500  fPAItransferTable = fPAIxscBank[jMat];
501
502  for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
503  {
504    if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;   
505  }
506  iPlace = iTkin - 1;
507  if(iPlace < 0) iPlace = 0;
508
509  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
510  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl; 
511  cross1 = GetdNdxCut(iPlace,tmax) ;
512  // G4cout<<"cross1 = "<<cross1<<G4endl; 
513  cross2 = GetdNdxCut(iPlace,cutEnergy) ; 
514  // G4cout<<"cross2 = "<<cross2<<G4endl; 
515  cross  = (cross2-cross1)*charge2;
516  // G4cout<<"cross = "<<cross<<G4endl; 
517  if( cross < DBL_MIN) cross = DBL_MIN;
518  //  if( cross2 < DBL_MIN) cross2 = DBL_MIN;
519
520  // return cross2;
521  return cross;
522}
523
524///////////////////////////////////////////////////////////////////////////
525//
526// It is analog of PostStepDoIt in terms of secondary electron.
527//
528 
529void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
530                                   const G4MaterialCutsCouple* matCC,
531                                   const G4DynamicParticle* dp,
532                                   G4double tmin,
533                                   G4double maxEnergy)
534{
535  size_t jMat;
536  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
537  {
538    if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
539  }
540  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
541
542  fPAItransferTable = fPAIxscBank[jMat];
543  fdNdxCutVector    = fdNdxCutTable[jMat];
544
545  G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy);
546  if( tmin >= tmax && fVerbose > 0)
547  {
548    G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
549  }
550  G4ThreeVector direction= dp->GetMomentumDirection();
551  G4double particleMass  = dp->GetMass();
552  G4double kineticEnergy = dp->GetKineticEnergy();
553
554  G4double massRatio     = fMass/particleMass;
555  G4double scaledTkin    = kineticEnergy*massRatio;
556  G4double totalEnergy   = kineticEnergy + particleMass;
557  G4double pSquare       = kineticEnergy*(totalEnergy+particleMass);
558 
559  G4double deltaTkin     = GetPostStepTransfer(scaledTkin);
560
561  // G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl;
562
563  if( deltaTkin <= 0. && fVerbose > 0) 
564  {
565    G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
566  }
567  if( deltaTkin <= 0.) return;
568
569  if( deltaTkin > tmax) deltaTkin = tmax;
570
571  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
572  G4double totalMomentum      = sqrt(pSquare);
573  G4double costheta           = deltaTkin*(totalEnergy + electron_mass_c2)
574                                /(deltaTotalMomentum * totalMomentum);
575
576  if( costheta > 0.99999 ) costheta = 0.99999;
577  G4double sintheta = 0.0;
578  G4double sin2 = 1. - costheta*costheta;
579  if( sin2 > 0.) sintheta = sqrt(sin2);
580
581  //  direction of the delta electron
582  G4double phi  = twopi*G4UniformRand(); 
583  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
584
585  G4ThreeVector deltaDirection(dirx,diry,dirz);
586  deltaDirection.rotateUz(direction);
587  deltaDirection.unit();
588
589  // primary change
590  kineticEnergy -= deltaTkin;
591  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
592  direction = dir.unit();
593  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
594  fParticleChange->SetProposedMomentumDirection(direction);
595
596  // create G4DynamicParticle object for e- delta ray
597  G4DynamicParticle* deltaRay = new G4DynamicParticle;
598  deltaRay->SetDefinition(G4Electron::Electron());
599  deltaRay->SetKineticEnergy( deltaTkin );  //  !!! trick for last steps /2.0 ???
600  deltaRay->SetMomentumDirection(deltaDirection); 
601
602  vdp->push_back(deltaRay);
603}
604
605
606///////////////////////////////////////////////////////////////////////
607//
608// Returns post step PAI energy transfer > cut electron energy according to passed
609// scaled kinetic energy of particle
610
611G4double 
612G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
613{ 
614  // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
615
616  G4int iTkin, iTransfer, iPlace  ;
617  G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
618
619  for(iTkin=0;iTkin<fTotBin;iTkin++)
620  {
621    if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin))  break ;
622  }
623  iPlace = iTkin - 1 ;
624  // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
625  if(iPlace < 0) iPlace = 0;
626  dNdxCut1 = (*fdNdxCutVector)(iPlace) ; 
627  // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
628
629
630  if(iTkin == fTotBin) // Fermi plato, try from left
631  {
632      position = dNdxCut1*G4UniformRand() ;
633
634      for( iTransfer = 0;
635 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
636      {
637        if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
638      }
639      transfer = GetEnergyTransfer(iPlace,position,iTransfer);
640  }
641  else
642  {
643    dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ; 
644    // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
645    if(iTkin == 0) // Tkin is too small, trying from right only
646    {
647      position = dNdxCut2*G4UniformRand() ;
648
649      for( iTransfer = 0;
650  iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
651      {
652        if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
653      }
654      transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
655    } 
656    else // general case: Tkin between two vectors of the material
657    {
658      E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 
659      E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin)     ;
660      W  = 1.0/(E2 - E1) ;
661      W1 = (E2 - scaledTkin)*W ;
662      W2 = (scaledTkin - E1)*W ;
663
664      position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
665
666        // G4cout<<position<<"\t" ;
667
668      G4int iTrMax1, iTrMax2, iTrMax;
669
670      iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
671      iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength());
672
673      if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
674      else                    iTrMax = iTrMax1;
675
676
677      for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
678      {
679          if( position >=
680          ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
681            (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
682      }
683      transfer = GetEnergyTransfer(iPlace,position,iTransfer);
684    }
685  } 
686  // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
687  if(transfer < 0.0 ) transfer = 0.0 ;
688  // if(transfer < DBL_MIN ) transfer = DBL_MIN;
689
690  return transfer ;
691}
692
693///////////////////////////////////////////////////////////////////////
694//
695// Returns random PAI energy transfer according to passed
696// indexes of particle kinetic
697
698G4double
699G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
700{ 
701  G4int iTransferMax;
702  G4double x1, x2, y1, y2, energyTransfer;
703
704  if(iTransfer == 0)
705  {
706    energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
707  } 
708  else
709  {
710    iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
711
712    if ( iTransfer >= iTransferMax )  iTransfer = iTransferMax - 1;
713   
714    y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
715    y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
716
717    x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
718    x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
719
720    if ( x1 == x2 )    energyTransfer = x2;
721    else
722    {
723      if ( y1 == y2  ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
724      else
725      {
726        energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
727      }
728    }
729  }
730  return energyTransfer;
731}
732
733///////////////////////////////////////////////////////////////////////
734
735G4double G4PAIModel::SampleFluctuations( const G4Material* material,
736                                         const G4DynamicParticle* aParticle,
737                                               G4double&,
738                                               G4double& step,
739                                               G4double&)
740{
741  size_t jMat;
742  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
743  {
744    if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
745  }
746  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
747
748  fPAItransferTable = fPAIxscBank[jMat];
749  fdNdxCutVector   = fdNdxCutTable[jMat];
750
751  G4int iTkin, iTransfer, iPlace  ;
752  G4long numOfCollisions=0;
753
754  //  G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
755  //G4cout<<"in:  "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl ;
756
757  G4double loss = 0.0, charge2 ;
758  G4double stepSum = 0., stepDelta, lambda, omega; 
759  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
760  G4bool numb = true;
761  G4double Tkin       = aParticle->GetKineticEnergy() ;
762  G4double MassRatio  = fMass/aParticle->GetDefinition()->GetPDGMass() ;
763  G4double charge     = aParticle->GetDefinition()->GetPDGCharge() ;
764  charge2             = charge*charge ;
765  G4double TkinScaled = Tkin*MassRatio ;
766
767  for(iTkin=0;iTkin<fTotBin;iTkin++)
768  {
769    if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin))   break ;
770  }
771  iPlace = iTkin - 1 ; 
772  if(iPlace < 0) iPlace = 0;
773  //  G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
774  dNdxCut1 = (*fdNdxCutVector)(iPlace) ; 
775  //  G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
776
777
778  if(iTkin == fTotBin) // Fermi plato, try from left
779  {
780    meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2;
781    if(meanNumber < 0.) meanNumber = 0. ;
782    //  numOfCollisions = RandPoisson::shoot(meanNumber) ;
783    // numOfCollisions = G4Poisson(meanNumber) ;
784    if( meanNumber > 0.) lambda = step/meanNumber;
785    else                 lambda = DBL_MAX;
786    while(numb)
787    {
788     stepDelta = CLHEP::RandExponential::shoot(lambda);
789     stepSum += stepDelta;
790     if(stepSum >= step) break;
791     numOfCollisions++;
792    }   
793    //    G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
794
795    while(numOfCollisions)
796    {
797      position = dNdxCut1+
798                 ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
799
800      for( iTransfer = 0;
801   iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
802      {
803        if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
804      }
805      omega = GetEnergyTransfer(iPlace,position,iTransfer);
806      // G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
807      loss += omega;
808      numOfCollisions-- ;
809    }
810  }
811  else
812  {
813    dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ; 
814    //  G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
815 
816    if(iTkin == 0) // Tkin is too small, trying from right only
817    {
818      meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2;
819      if( meanNumber < 0. ) meanNumber = 0. ;
820      //  numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
821      //  numOfCollisions = G4Poisson(meanNumber) ;
822    if( meanNumber > 0.) lambda = step/meanNumber;
823    else                 lambda = DBL_MAX;
824    while(numb)
825    {
826     stepDelta = CLHEP::RandExponential::shoot(lambda);
827     stepSum += stepDelta;
828     if(stepSum >= step) break;
829     numOfCollisions++;
830    }   
831
832    //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
833
834      while(numOfCollisions)
835      {
836        position = dNdxCut2+
837                   ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
838   
839        for( iTransfer = 0;
840   iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
841        {
842          if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
843        }
844        omega = GetEnergyTransfer(iPlace,position,iTransfer);
845        // G4cout<<omega/keV<<"\t";
846        loss += omega;
847        numOfCollisions-- ;
848      }
849    } 
850    else // general case: Tkin between two vectors of the material
851    {
852      E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 
853      E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin)     ;
854       W = 1.0/(E2 - E1) ;
855      W1 = (E2 - TkinScaled)*W ;
856      W2 = (TkinScaled - E1)*W ;
857
858      // G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
859      //   (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
860      // G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
861      //     (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
862
863      meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 + 
864                   ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2;
865      if(meanNumber<0.0) meanNumber = 0.0;
866      //  numOfCollisions = RandPoisson::shoot(meanNumber) ;
867      // numOfCollisions = G4Poisson(meanNumber) ;
868    if( meanNumber > 0.) lambda = step/meanNumber;
869    else                 lambda = DBL_MAX;
870    while(numb)
871    {
872     stepDelta = CLHEP::RandExponential::shoot(lambda);
873     stepSum += stepDelta;
874     if(stepSum >= step) break;
875     numOfCollisions++;
876    }   
877
878    //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
879
880      while(numOfCollisions)
881      {
882        position = dNdxCut1*W1 + dNdxCut2*W2 +
883                 ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 + 
884                    dNdxCut2+
885                  ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand();
886
887        // G4cout<<position<<"\t" ;
888
889        for( iTransfer = 0;
890    iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
891        {
892          if( position >=
893          ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 + 
894            (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
895          {
896              break ;
897          }
898        }
899        omega =  GetEnergyTransfer(iPlace,position,iTransfer);
900        //  G4cout<<omega/keV<<"\t";
901        loss += omega;
902        numOfCollisions-- ;   
903      }
904    }
905  } 
906  //  G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
907  //  <<step/mm<<" mm"<<G4endl ;
908  if(loss > Tkin) loss=Tkin;
909  if(loss < 0.  ) loss = 0.;
910  return loss ;
911
912}
913
914//////////////////////////////////////////////////////////////////////
915//
916// Returns the statistical estimation of the energy loss distribution variance
917//
918
919
920G4double G4PAIModel::Dispersion( const G4Material* material, 
921                                 const G4DynamicParticle* aParticle,
922                                       G4double& tmax, 
923                                       G4double& step       )
924{
925  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
926  for(G4int i = 0 ; i < fMeanNumber; i++)
927  {
928    loss      = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
929    sumLoss  += loss;
930    sumLoss2 += loss*loss;
931  }
932  meanLoss = sumLoss/fMeanNumber;
933  sigma2   = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
934  return sigma2;
935}
936
937
938//
939//
940/////////////////////////////////////////////////
941
942
943
944
945
946
Note: See TracBrowser for help on using the repository browser.