G4ScreenedNuclearRecoil.hh

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ScreenedNuclearRecoil.hh

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
//
//
// G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp
// GEANT4 tag 
//
//
//
// Class Description
// Process for screened electromagnetic nuclear elastic scattering; 
// Physics comes from:
// Marcus H. Mendenhall and Robert A. Weller, 
// "Algorithms  for  the rapid  computation  of  classical  cross  sections  
// for  screened  Coulomb  collisions  "
// Nuclear  Instruments  and  Methods  in  Physics  Research  B58  (1991)  11-17
// The only input required is a screening function phi(r/a) which is the ratio
// of the actual interatomic potential for two atoms with atomic numbers Z1 and
// Z2,
// to the unscreened potential Z1*Z2*e^2/r where e^2 is elm_coupling in Geant4 
// units
// the actual screening tables are computed externally in a python module 
// "screened_scattering.py"
// to allow very specific screening functions to be added if desired, without 
// messing with the insides of this code.
//
// First version, April 2004, Marcus H. Mendenhall, Vanderbilt University
// May 1, 2008 -- Added code to allow process to have zero cross section above 
//  max energy, to coordinate with G4MSC.  -- mhm
//
// Class Description - End


#ifndef G4ScreenedNuclearRecoil_h
#define G4ScreenedNuclearRecoil_h 1

#include "globals.hh"
#include "G4VDiscreteProcess.hh"
#include "G4ParticleChange.hh"
#include "c2_function.hh"

#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"

#include <map>
#include <vector>

class G4VNIELPartition;

typedef c2_const_ptr<G4double> G4_c2_const_ptr;
typedef c2_ptr<G4double> G4_c2_ptr;
typedef c2_function<G4double> G4_c2_function;

typedef struct G4ScreeningTables { 
        G4double z1, z2, m1, m2, au, emin;
        G4_c2_const_ptr EMphiData; 
} G4ScreeningTables;

// A class for loading ScreenedCoulombCrossSections
class G4ScreenedCoulombCrossSectionInfo
{
public:
        G4ScreenedCoulombCrossSectionInfo() { }
        ~G4ScreenedCoulombCrossSectionInfo() { }
        
        static const char* CVSHeaderVers() { return 
 "G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp GEANT4 tag ";
        }
        static const char* CVSFileVers();
};

// A class for loading ScreenedCoulombCrossSections
class G4ScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSectionInfo
{
public:

      G4ScreenedCoulombCrossSection() : verbosity(1) { }
      G4ScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection &src) :
                G4ScreenedCoulombCrossSectionInfo(),verbosity(src.verbosity) { }
      virtual ~G4ScreenedCoulombCrossSection();
        
      typedef std::map<G4int, G4ScreeningTables> ScreeningMap;
        
      // a local, fast-access mapping of a particle's Z to its full definition
      typedef std::map<G4int, class G4ParticleDefinition *> ParticleCache;
        
      // LoadData is called by G4ScreenedNuclearRecoil::GetMeanFreePath
      // It loads the data tables, builds the elemental cross-section tables.
      virtual void LoadData(G4String screeningKey, G4int z1, G4double m1,
                            G4double recoilCutoff) = 0;
        
      // BuildMFPTables is called by G4ScreenedNuclearRecoil::GetMeanFreePath 
      //to build the MFP tables for each material
      void BuildMFPTables(void); // scan the MaterialsTable and construct MFP 
                                 //tables
        
      virtual G4ScreenedCoulombCrossSection *create() = 0; 
      // a 'virtual constructor' which clones the class
      const G4ScreeningTables *GetScreening(G4int Z)
                     { return &(screeningData[Z]); }
      void SetVerbosity(G4int v) { verbosity=v; }
        
      // this process needs element selection weighted only by number density
      G4ParticleDefinition* SelectRandomUnweightedTarget
                                    (const G4MaterialCutsCouple* couple);
        
      enum { nMassMapElements=116 };
                
      G4double standardmass(G4int z1) 
            { return z1 <= nMassMapElements ? massmap[z1] : 2.5*z1; }
        
      // get the mean-free-path table for the indexed material 
      const G4_c2_function * operator [] (G4int materialIndex) { 
              return MFPTables.find(materialIndex)!=MFPTables.end() ?
                     &(MFPTables[materialIndex].get()) : (G4_c2_function *)0;
      }
        
protected:
      ScreeningMap screeningData; // screening tables for each element
      ParticleCache targetMap;
      G4int verbosity;
      std::map<G4int, G4_c2_const_ptr > sigmaMap; 
      // total cross section for each element
      std::map<G4int, G4_c2_const_ptr > MFPTables; // MFP for each material
        
private:
      static const G4double massmap[nMassMapElements+1];

};

typedef struct G4CoulombKinematicsInfo {
        G4double impactParameter;
        G4ScreenedCoulombCrossSection *crossSection;
        G4double a1, a2, sinTheta, cosTheta, sinZeta, cosZeta, eRecoil;
        G4ParticleDefinition *recoilIon;        
        const G4Material *targetMaterial;
} G4CoulombKinematicsInfo;

class G4ScreenedCollisionStage {
public:
     virtual void DoCollisionStep(class G4ScreenedNuclearRecoil *master,
                    const class G4Track& aTrack, const class G4Step& aStep)=0;
     virtual ~G4ScreenedCollisionStage() {}
};

class G4ScreenedCoulombClassicalKinematics: 
  public G4ScreenedCoulombCrossSectionInfo, public G4ScreenedCollisionStage {

public:
     G4ScreenedCoulombClassicalKinematics();
     virtual void DoCollisionStep(class G4ScreenedNuclearRecoil *master,
             const class G4Track& aTrack, const class G4Step& aStep);
        
     G4bool DoScreeningComputation(class G4ScreenedNuclearRecoil *master,
             const G4ScreeningTables *screen,  G4double eps, G4double beta);
             
     virtual ~G4ScreenedCoulombClassicalKinematics() { }
     
protected:
     // the c2_functions we need to do the work.
     c2_const_plugin_function_p<G4double> &phifunc;
     c2_linear_p<G4double> &xovereps;
     G4_c2_ptr diff;        
};

class G4SingleScatter: public G4ScreenedCoulombCrossSectionInfo,
                       public G4ScreenedCollisionStage {

public:
     G4SingleScatter() { }
     virtual void DoCollisionStep(class G4ScreenedNuclearRecoil *master,
             const class G4Track& aTrack, const class G4Step& aStep);
     virtual ~G4SingleScatter() {}
};

class G4ScreenedNuclearRecoil : public G4ScreenedCoulombCrossSectionInfo,
                                public G4VDiscreteProcess
{
public:
        
     friend class G4ScreenedCollisionStage;

     
     G4ScreenedNuclearRecoil(const G4String& processName = "ScreenedElastic",
                const G4String &ScreeningKey="zbl", G4bool GenerateRecoils=1, 
                G4double RecoilCutoff=100.0*eV, G4double PhysicsCutoff=10.0*eV);
                
     virtual ~G4ScreenedNuclearRecoil();
     
     virtual G4double GetMeanFreePath(const G4Track&, G4double,
                                      G4ForceCondition* );
                                      
     virtual G4VParticleChange* PostStepDoIt(const G4Track& aTrack, 
                                             const G4Step& aStep);
     virtual G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
     virtual void BuildPhysicsTable(const G4ParticleDefinition& aParticleType);
     virtual void DumpPhysicsTable(const G4ParticleDefinition& aParticleType);
     
     virtual G4bool CheckNuclearCollision(G4double A, G4double A1,
                            G4double apsis); // return true if hard collision

     virtual G4ScreenedCoulombCrossSection *GetNewCrossSectionHandler(void);
        
     G4double GetNIEL() const { return NIEL; } 
        
     void ResetTables(); 
     // clear all data tables to allow changing energy cutoff, materials, etc.

     
     void SetMaxEnergyForScattering(G4double energy) {processMaxEnergy=energy;}
     
     std::string GetScreeningKey() const { return screeningKey; }
     
     
     void AllowEnergyDeposition(G4bool flag) {registerDepositedEnergy=flag;}
     
     G4bool GetAllowEnergyDeposition() const {return registerDepositedEnergy;}
     
     
     void EnableRecoils(G4bool flag) { generateRecoils=flag; }
     
     G4bool GetEnableRecoils() const { return generateRecoils; }
     
     
     void SetMFPScaling(G4double scale) { MFPScale=scale; }
     
     G4double GetMFPScaling() const { return MFPScale; }
     
     
     void AvoidNuclearReactions(G4bool flag) { avoidReactions=flag; }
     
     G4bool GetAvoidNuclearReactions() const { return avoidReactions; }
     
     
     void SetRecoilCutoff(G4double energy) { recoilCutoff=energy; }
     
     G4double GetRecoilCutoff() const { return recoilCutoff; }
     
     
     void SetPhysicsCutoff(G4double energy) { physicsCutoff=energy; 
                                              ResetTables(); }
                                              
     G4double GetPhysicsCutoff() const { return physicsCutoff; }
     
     
     void SetNIELPartitionFunction(const G4VNIELPartition *part);
     
     
     void SetCrossSectionHardening(G4double fraction, G4double HardeningFactor){
                hardeningFraction=fraction;
                hardeningFactor=HardeningFactor;
     }
     
     G4double GetHardeningFraction() const { return hardeningFraction; }
     
     G4double GetHardeningFactor() const { return hardeningFactor; }
     
     G4double GetCurrentInteractionLength() const {
                                    return currentInteractionLength; }
                                    
     
     void SetExternalCrossSectionHandler(G4ScreenedCoulombCrossSection *cs) { 
                            externalCrossSectionConstructor=cs; 
     }
     
     G4int GetVerboseLevel() const { return verboseLevel; }

     std::map<G4int, G4ScreenedCoulombCrossSection*> &GetCrossSectionHandlers() 
                { return crossSectionHandlers; }
                
     void ClearStages(void);
      
     void AddStage(G4ScreenedCollisionStage *stage) { 
                                     collisionStages.push_back(stage); }
                                     
     G4CoulombKinematicsInfo &GetKinematics() { return kinematics; }
     
     void SetValidCollision(G4bool flag) { validCollision=flag; }
     
     G4bool GetValidCollision() const { return validCollision; }

     class G4ParticleChange &GetParticleChange() 
              { return static_cast<G4ParticleChange &>(*pParticleChange); }
              
     
     void DepositEnergy(G4int z1, G4double a1, const G4Material *material,
                        G4double energy);
        
protected:
     G4double highEnergyLimit;
     
     G4double lowEnergyLimit;
     
     G4double processMaxEnergy;
     G4String screeningKey;
     G4bool generateRecoils, avoidReactions;
     G4double recoilCutoff, physicsCutoff;
     G4bool registerDepositedEnergy;
     G4double IonizingLoss, NIEL;
     G4double MFPScale;
     G4double hardeningFraction, hardeningFactor;
        
     G4ScreenedCoulombCrossSection *externalCrossSectionConstructor;
     std::vector<G4ScreenedCollisionStage *> collisionStages;
        
     std::map<G4int, G4ScreenedCoulombCrossSection*> crossSectionHandlers;
        
     G4bool validCollision;
     G4CoulombKinematicsInfo kinematics;
     const G4VNIELPartition *NIELPartitionFunction;
};

// A customized G4CrossSectionHandler which gets its data from 
// an external program

class G4NativeScreenedCoulombCrossSection: public G4ScreenedCoulombCrossSection 
{
public:
     G4NativeScreenedCoulombCrossSection();

     G4NativeScreenedCoulombCrossSection(
                              const G4NativeScreenedCoulombCrossSection &src) 
          : G4ScreenedCoulombCrossSection(src), phiMap(src.phiMap) { }
        
     G4NativeScreenedCoulombCrossSection(
                                     const G4ScreenedCoulombCrossSection &src) 
          : G4ScreenedCoulombCrossSection(src)  { }
          
     virtual ~G4NativeScreenedCoulombCrossSection();
        
     virtual void LoadData(G4String screeningKey, G4int z1, G4double m1,
                           G4double recoilCutoff);
                           
     virtual G4ScreenedCoulombCrossSection *create() 
                { return new G4NativeScreenedCoulombCrossSection(*this); } 
                
     // get a list of available keys
     std::vector<G4String> GetScreeningKeys() const;
        
     typedef G4_c2_function &(*ScreeningFunc)(G4int z1, G4int z2, 
                                 size_t nPoints, G4double rMax, G4double *au);
        
     void AddScreeningFunction(G4String name, ScreeningFunc fn) {
                           phiMap[name]=fn;
     }
        
private:
     // this is a map used to look up screening function generators
     std::map<std::string, ScreeningFunc> phiMap;
};

G4_c2_function &ZBLScreening(G4int z1, G4int z2, size_t npoints, 
                             G4double rMax, G4double *auval);
                             
G4_c2_function &MoliereScreening(G4int z1, G4int z2, size_t npoints, 
                                 G4double rMax, G4double *auval);
                                 
G4_c2_function &LJScreening(G4int z1, G4int z2, size_t npoints,
                            G4double rMax, G4double *auval);
                                                        
G4_c2_function &LJZBLScreening(G4int z1, G4int z2, size_t npoints,
                               G4double rMax, G4double *auval);

#endif