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

Last change on this file since 1006 was 1005, checked in by garnier, 15 years ago

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