physics/src/G4eLowEnergyLoss.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // -----------------------------------------------------------
0029 //      GEANT 4 class implementation file
0030 //
0031 //      History: based on object model of
0032 //      2nd December 1995, G.Cosmo
0033 //      ---------- G4eLowEnergyLoss physics process -----------
0034 //                by Laszlo Urban, 20 March 1997
0035 // **************************************************************
0036 // It calculates the energy loss of e+/e-.
0037 // --------------------------------------------------------------
0038 //
0039 // 08-05-97: small changes by L.Urban
0040 // 27-05-98: several bugs and inconsistencies are corrected,
0041 //           new table (the inverse of the range table) added ,
0042 //           AlongStepDoit uses now this new table. L.Urban
0043 // 08-09-98: cleanup
0044 // 24-09-98: rndmStepFlag false by default (no randomization of the step)
0045 // 14-10-98: messenger file added.
0046 // 16-10-98: public method SetStepFunction()
0047 // 20-01-99: important correction in AlongStepDoIt , L.Urban
0048 // 10/02/00  modifications , new e.m. structure, L.Urban
0049 // 11/04/00: Bug fix in dE/dx fluctuation simulation, Veronique Lefebure
0050 // 19-09-00  change of fluctuation sampling V.Ivanchenko
0051 // 20/09/00  update fluctuations V.Ivanchenko
0052 // 18/10/01  add fluorescence AlongStepDoIt V.Ivanchenko
0053 // 18/10/01  Revision to improve code quality and consistency with design, MGP
0054 // 19/10/01  update according to new design, V.Ivanchenko
0055 // 24/10/01  MGP - Protection against negative energy loss in AlongStepDoIt
0056 // 26/10/01  VI Clean up access to deexcitation
0057 // 23/11/01  VI Move static member-functions from header to source
0058 // 28/05/02  VI Remove flag fStopAndKill
0059 // 03/06/02  MGP - Restore fStopAndKill
0060 // 28/10/02  VI Optimal binning for dE/dx
0061 // 21/01/03  VI cut per region
0062 // 01/06/04  VI check if stopped particle has AtRest processes
0063 //
0064 // --------------------------------------------------------------
0065
0066 #include "G4eLowEnergyLoss.hh"
0067 #include "G4SystemOfUnits.hh"
0068 #include "G4EnergyLossMessenger.hh"
0069 #include "G4Poisson.hh"
0070 #include "G4ProductionCutsTable.hh"
0071
0072 //
0073
0074 // Initialisation of static data members
0075 // -------------------------------------
0076 // Contributing processes : ion.loss + soft brems->NbOfProcesses is initialized
0077 // to 2 . YOU DO NOT HAVE TO CHANGE this variable for a 'normal' run.
0078 //
0079 // You have to change NbOfProcesses if you invent a new process contributing
0080 // to the continuous energy loss.
0081 // The NbOfProcesses data member can be changed using the (public static)
0082 // functions Get/Set/Plus/MinusNbOfProcesses (see G4eLowEnergyLoss.hh)
0083
0084 G4int            G4eLowEnergyLoss::NbOfProcesses = 2;
0085
0086 G4int            G4eLowEnergyLoss::CounterOfElectronProcess = 0;
0087 G4int            G4eLowEnergyLoss::CounterOfPositronProcess = 0;
0088 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfElectronProcess =
0089                                            new G4PhysicsTable*[10];
0090 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfPositronProcess =
0091                                            new G4PhysicsTable*[10];
0092
0093
0094 G4PhysicsTable*  G4eLowEnergyLoss::theDEDXElectronTable         = 0;
0095 G4PhysicsTable*  G4eLowEnergyLoss::theDEDXPositronTable         = 0;
0096 G4PhysicsTable*  G4eLowEnergyLoss::theRangeElectronTable        = 0;
0097 G4PhysicsTable*  G4eLowEnergyLoss::theRangePositronTable        = 0;
0098 G4PhysicsTable*  G4eLowEnergyLoss::theInverseRangeElectronTable = 0;
0099 G4PhysicsTable*  G4eLowEnergyLoss::theInverseRangePositronTable = 0;
0100 G4PhysicsTable*  G4eLowEnergyLoss::theLabTimeElectronTable      = 0;
0101 G4PhysicsTable*  G4eLowEnergyLoss::theLabTimePositronTable      = 0;
0102 G4PhysicsTable*  G4eLowEnergyLoss::theProperTimeElectronTable   = 0;
0103 G4PhysicsTable*  G4eLowEnergyLoss::theProperTimePositronTable   = 0;
0104
0105 G4PhysicsTable*  G4eLowEnergyLoss::theeRangeCoeffATable         = 0;
0106 G4PhysicsTable*  G4eLowEnergyLoss::theeRangeCoeffBTable         = 0;
0107 G4PhysicsTable*  G4eLowEnergyLoss::theeRangeCoeffCTable         = 0;
0108 G4PhysicsTable*  G4eLowEnergyLoss::thepRangeCoeffATable         = 0;
0109 G4PhysicsTable*  G4eLowEnergyLoss::thepRangeCoeffBTable         = 0;
0110 G4PhysicsTable*  G4eLowEnergyLoss::thepRangeCoeffCTable         = 0;
0111
0112 G4double         G4eLowEnergyLoss::LowerBoundEloss = 10.*eV ;
0113 G4double         G4eLowEnergyLoss::UpperBoundEloss = 100.*GeV ;
0114 G4int            G4eLowEnergyLoss::NbinEloss = 360 ;
0115 G4double         G4eLowEnergyLoss::RTable ;
0116 G4double         G4eLowEnergyLoss::LOGRTable ;
0117
0118
0119 G4EnergyLossMessenger* G4eLowEnergyLoss::eLossMessenger         = 0;
0120
0121 //
0122
0123 // constructor and destructor
0124
0125 G4eLowEnergyLoss::G4eLowEnergyLoss(const G4String& processName)
0126    : G4RDVeLowEnergyLoss (processName),
0127      theLossTable(0),
0128      MinKineticEnergy(1.*eV),
0129      Charge(-1.),
0130      lastCharge(0.),
0131      theDEDXTable(0),
0132      CounterOfProcess(0),
0133      RecorderOfProcess(0),
0134      fdEdx(0),
0135      fRangeNow(0),
0136      linLossLimit(0.05),
0137      theFluo(false)
0138 {
0139
0140  //create (only once) EnergyLoss messenger
0141  if(!eLossMessenger) eLossMessenger = new G4EnergyLossMessenger();
0142 }
0143
0144 //
0145
0146 G4eLowEnergyLoss::~G4eLowEnergyLoss()
0147 {
0148      if (theLossTable)
0149        {
0150          theLossTable->clearAndDestroy();
0151          delete theLossTable;
0152        }
0153 }
0154
0155 void G4eLowEnergyLoss::SetNbOfProcesses(G4int nb)
0156 {
0157     NbOfProcesses=nb;
0158 }
0159
0160 void G4eLowEnergyLoss::PlusNbOfProcesses()
0161 {
0162     NbOfProcesses++;
0163 }
0164
0165 void G4eLowEnergyLoss::MinusNbOfProcesses()
0166 {
0167     NbOfProcesses--;
0168 }
0169
0170 G4int G4eLowEnergyLoss::GetNbOfProcesses()
0171 {
0172     return NbOfProcesses;
0173 }
0174
0175 void G4eLowEnergyLoss::SetLowerBoundEloss(G4double val)
0176 {
0177     LowerBoundEloss=val;
0178 }
0179
0180 void G4eLowEnergyLoss::SetUpperBoundEloss(G4double val)
0181 {
0182     UpperBoundEloss=val;
0183 }
0184
0185 void G4eLowEnergyLoss::SetNbinEloss(G4int nb)
0186 {
0187     NbinEloss=nb;
0188 }
0189
0190 G4double G4eLowEnergyLoss::GetLowerBoundEloss()
0191 {
0192     return LowerBoundEloss;
0193 }
0194
0195 G4double G4eLowEnergyLoss::GetUpperBoundEloss()
0196 {
0197     return UpperBoundEloss;
0198 }
0199
0200 G4int G4eLowEnergyLoss::GetNbinEloss()
0201 {
0202     return NbinEloss;
0203 }
0204 //
0205
0206 void G4eLowEnergyLoss::BuildDEDXTable(
0207                          const G4ParticleDefinition& aParticleType)
0208 {
0209   ParticleMass = aParticleType.GetPDGMass();
0210   Charge = aParticleType.GetPDGCharge()/eplus;
0211
0212   //  calculate data members LOGRTable,RTable first
0213
0214   G4double lrate = std::log(UpperBoundEloss/LowerBoundEloss);
0215   LOGRTable=lrate/NbinEloss;
0216   RTable   =std::exp(LOGRTable);
0217   // Build energy loss table as a sum of the energy loss due to the
0218   // different processes.
0219   //
0220
0221   const G4ProductionCutsTable* theCoupleTable=
0222         G4ProductionCutsTable::GetProductionCutsTable();
0223   size_t numOfCouples = theCoupleTable->GetTableSize();
0224
0225   // create table for the total energy loss
0226
0227   if (&aParticleType==G4Electron::Electron())
0228     {
0229       RecorderOfProcess=RecorderOfElectronProcess;
0230       CounterOfProcess=CounterOfElectronProcess;
0231       if (CounterOfProcess == NbOfProcesses)
0232         {
0233          if (theDEDXElectronTable)
0234            {
0235              theDEDXElectronTable->clearAndDestroy();
0236              delete theDEDXElectronTable;
0237            }
0238          theDEDXElectronTable = new G4PhysicsTable(numOfCouples);
0239          theDEDXTable = theDEDXElectronTable;
0240         }
0241     }
0242   if (&aParticleType==G4Positron::Positron())
0243     {
0244      RecorderOfProcess=RecorderOfPositronProcess;
0245      CounterOfProcess=CounterOfPositronProcess;
0246      if (CounterOfProcess == NbOfProcesses)
0247        {
0248         if (theDEDXPositronTable)
0249           {
0250             theDEDXPositronTable->clearAndDestroy();
0251             delete theDEDXPositronTable;
0252           }
0253         theDEDXPositronTable = new G4PhysicsTable(numOfCouples);
0254         theDEDXTable = theDEDXPositronTable;
0255        }
0256     }
0257
0258   if (CounterOfProcess == NbOfProcesses)
0259     {
0260      // fill the tables
0261      // loop for materials
0262      G4double LowEdgeEnergy , Value;
0263      G4bool isOutRange;
0264      G4PhysicsTable* pointer;
0265
0266      for (size_t J=0; J<numOfCouples; J++)
0267         {
0268          // create physics vector and fill it
0269
0270          G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
0271                     LowerBoundEloss, UpperBoundEloss, NbinEloss);
0272
0273          // loop for the kinetic energy
0274          for (G4int i=0; i<NbinEloss; i++)
0275             {
0276               LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
0277               //here comes the sum of the different tables created by the
0278               //processes (ionisation,bremsstrahlung,etc...)
0279               Value = 0.;
0280               for (G4int process=0; process < NbOfProcesses; process++)
0281                  {
0282                    pointer= RecorderOfProcess[process];
0283                    Value += (*pointer)[J]->GetValue(LowEdgeEnergy,isOutRange);
0284                  }
0285
0286               aVector->PutValue(i,Value) ;
0287             }
0288
0289          theDEDXTable->insert(aVector) ;
0290
0291         }
0292
0293
0294      //reset counter to zero
0295      if (&aParticleType==G4Electron::Electron()) CounterOfElectronProcess=0;
0296      if (&aParticleType==G4Positron::Positron()) CounterOfPositronProcess=0;
0297
0298      ParticleMass = aParticleType.GetPDGMass();
0299
0300      if (&aParticleType==G4Electron::Electron())
0301      {
0302        // Build range table
0303        theRangeElectronTable = BuildRangeTable(theDEDXElectronTable,
0304                                                theRangeElectronTable,
0305                               LowerBoundEloss,UpperBoundEloss,NbinEloss);
0306
0307        // Build lab/proper time tables
0308        theLabTimeElectronTable = BuildLabTimeTable(theDEDXElectronTable,
0309                          theLabTimeElectronTable,
0310                          LowerBoundEloss,UpperBoundEloss,NbinEloss);
0311        theProperTimeElectronTable = BuildProperTimeTable(theDEDXElectronTable,
0312                             theProperTimeElectronTable,
0313                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
0314
0315        // Build coeff tables for the energy loss calculation
0316        theeRangeCoeffATable = BuildRangeCoeffATable(theRangeElectronTable,
0317                              theeRangeCoeffATable,
0318                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
0319
0320        theeRangeCoeffBTable = BuildRangeCoeffBTable(theRangeElectronTable,
0321                              theeRangeCoeffBTable,
0322                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
0323
0324        theeRangeCoeffCTable = BuildRangeCoeffCTable(theRangeElectronTable,
0325                              theeRangeCoeffCTable,
0326                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
0327
0328        // invert the range table
0329        theInverseRangeElectronTable = BuildInverseRangeTable(theRangeElectronTable,
0330                               theeRangeCoeffATable,
0331                               theeRangeCoeffBTable,
0332                               theeRangeCoeffCTable,
0333                               theInverseRangeElectronTable,
0334                               LowerBoundEloss,UpperBoundEloss,NbinEloss);
0335      }
0336      if (&aParticleType==G4Positron::Positron())
0337      {
0338        // Build range table
0339        theRangePositronTable = BuildRangeTable(theDEDXPositronTable,
0340                                                theRangePositronTable,
0341                               LowerBoundEloss,UpperBoundEloss,NbinEloss);
0342
0343
0344        // Build lab/proper time tables
0345        theLabTimePositronTable = BuildLabTimeTable(theDEDXPositronTable,
0346                          theLabTimePositronTable,
0347                          LowerBoundEloss,UpperBoundEloss,NbinEloss);
0348        theProperTimePositronTable = BuildProperTimeTable(theDEDXPositronTable,
0349                             theProperTimePositronTable,
0350                             LowerBoundEloss,UpperBoundEloss,NbinEloss);
0351
0352        // Build coeff tables for the energy loss calculation
0353        thepRangeCoeffATable = BuildRangeCoeffATable(theRangePositronTable,
0354                              thepRangeCoeffATable,
0355                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
0356
0357        thepRangeCoeffBTable = BuildRangeCoeffBTable(theRangePositronTable,
0358                              thepRangeCoeffBTable,
0359                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
0360
0361        thepRangeCoeffCTable = BuildRangeCoeffCTable(theRangePositronTable,
0362                              thepRangeCoeffCTable,
0363                              LowerBoundEloss,UpperBoundEloss,NbinEloss);
0364
0365        // invert the range table
0366        theInverseRangePositronTable = BuildInverseRangeTable(theRangePositronTable,
0367                               thepRangeCoeffATable,
0368                               thepRangeCoeffBTable,
0369                               thepRangeCoeffCTable,
0370                               theInverseRangePositronTable,
0371                               LowerBoundEloss,UpperBoundEloss,NbinEloss);
0372      }
0373
0374      if(verboseLevel > 1) {
0375        G4cout << (*theDEDXElectronTable) << G4endl;
0376      }
0377
0378
0379      // make the energy loss and the range table available
0380      G4EnergyLossTables::Register(&aParticleType,
0381        (&aParticleType==G4Electron::Electron())?
0382        theDEDXElectronTable: theDEDXPositronTable,
0383        (&aParticleType==G4Electron::Electron())?
0384        theRangeElectronTable: theRangePositronTable,
0385        (&aParticleType==G4Electron::Electron())?
0386        theInverseRangeElectronTable: theInverseRangePositronTable,
0387        (&aParticleType==G4Electron::Electron())?
0388        theLabTimeElectronTable: theLabTimePositronTable,
0389        (&aParticleType==G4Electron::Electron())?
0390        theProperTimeElectronTable: theProperTimePositronTable,
0391        LowerBoundEloss, UpperBoundEloss, 1.,NbinEloss);
0392     }
0393 }
0394
0395 //
0396
0397 G4VParticleChange* G4eLowEnergyLoss::AlongStepDoIt( const G4Track& trackData,
0398                                                  const G4Step&  stepData)
0399 {
0400  // compute the energy loss after a Step
0401
0402   static const G4double faclow = 1.5 ;
0403
0404   // get particle and material pointers from trackData
0405   const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
0406   G4double E      = aParticle->GetKineticEnergy();
0407
0408   const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
0409
0410   G4double Step = stepData.GetStepLength();
0411
0412   aParticleChange.Initialize(trackData);
0413   //fParticleChange.Initialize(trackData);
0414
0415   G4double MeanLoss, finalT;
0416
0417   if (E < MinKineticEnergy)   finalT = 0.;
0418
0419   else if ( E< faclow*LowerBoundEloss)
0420   {
0421     if (Step >= fRangeNow)  finalT = 0.;
0422    //  else finalT = E*(1.-Step/fRangeNow) ;
0423     else finalT = E*(1.-std::sqrt(Step/fRangeNow)) ;
0424   }
0425
0426   else if (E>=UpperBoundEloss) finalT = E - Step*fdEdx;
0427
0428   else if (Step >= fRangeNow)  finalT = 0.;
0429
0430   else
0431   {
0432     if(Step/fRangeNow < linLossLimit) finalT = E-Step*fdEdx ;
0433     else
0434     {
0435       if (Charge<0.) finalT = G4EnergyLossTables::GetPreciseEnergyFromRange
0436                              (G4Electron::Electron(),fRangeNow-Step,couple);
0437       else           finalT = G4EnergyLossTables::GetPreciseEnergyFromRange
0438                              (G4Positron::Positron(),fRangeNow-Step,couple);
0439      }
0440   }
0441
0442   if(finalT < MinKineticEnergy) finalT = 0. ;
0443
0444   MeanLoss = E-finalT ;
0445
0446   //now the loss with fluctuation
0447   if ((EnlossFlucFlag) && (finalT > 0.) && (finalT < E)&&(E > LowerBoundEloss))
0448   {
0449     finalT = E-GetLossWithFluct(aParticle,couple,MeanLoss,Step);
0450     if (finalT < 0.) finalT = 0.;
0451   }
0452
0453   // kill the particle if the kinetic energy <= 0
0454   if (finalT <= 0. )
0455   {
0456     finalT = 0.;
0457     if(Charge > 0.0) aParticleChange.ProposeTrackStatus(fStopButAlive);
0458     else             aParticleChange.ProposeTrackStatus(fStopAndKill);
0459   }
0460
0461   G4double edep = E - finalT;
0462
0463   aParticleChange.ProposeEnergy(finalT);
0464
0465   // Deexcitation of ionised atoms
0466   std::vector<G4DynamicParticle*>* deexcitationProducts = 0;
0467   if (theFluo) deexcitationProducts = DeexciteAtom(couple,E,edep);
0468
0469   size_t nSecondaries = 0;
0470   if (deexcitationProducts != 0) nSecondaries = deexcitationProducts->size();
0471   aParticleChange.SetNumberOfSecondaries(nSecondaries);
0472
0473   if (nSecondaries > 0) {
0474
0475     const G4StepPoint* preStep = stepData.GetPreStepPoint();
0476     const G4StepPoint* postStep = stepData.GetPostStepPoint();
0477     G4ThreeVector r = preStep->GetPosition();
0478     G4ThreeVector deltaR = postStep->GetPosition();
0479     deltaR -= r;
0480     G4double t = preStep->GetGlobalTime();
0481     G4double deltaT = postStep->GetGlobalTime();
0482     deltaT -= t;
0483     G4double time, q;
0484     G4ThreeVector position;
0485
0486     for (size_t i=0; i<nSecondaries; i++) {
0487
0488       G4DynamicParticle* part = (*deexcitationProducts)[i];
0489       if (part != 0) {
0490         G4double eSecondary = part->GetKineticEnergy();
0491         edep -= eSecondary;
0492     if (edep > 0.)
0493       {
0494         q = G4UniformRand();
0495         time = deltaT*q + t;
0496         position  = deltaR*q;
0497         position += r;
0498         G4Track* newTrack = new G4Track(part, time, position);
0499         aParticleChange.AddSecondary(newTrack);
0500       }
0501     else
0502       {
0503         edep += eSecondary;
0504         delete part;
0505         part = 0;
0506       }
0507       }
0508     }
0509   }
0510   delete deexcitationProducts;
0511
0512   aParticleChange.ProposeLocalEnergyDeposit(edep);
0513
0514   return &aParticleChange;
0515 }
0516
0517 //
0518