EnergyLossForExtrapolator.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <CLHEP/Units/SystemOfUnits.h>
 
#include <globals.hh>
#include <G4PhysicsTable.hh>
#include <G4DataVector.hh>
#include <vector>
 
class G4ParticleDefinition;
class G4Material;
class G4MaterialCutsCouple;
class G4ProductionCuts;
 
namespace Belle2 {
 
  namespace Simulation {

    class EnergyLossForExtrapolator {
 
    public:

      EnergyLossForExtrapolator();

      ~EnergyLossForExtrapolator();

      G4double ComputeDEDX(G4double kinEnergy, const G4ParticleDefinition*);

      G4double ComputeRange(G4double kinEnergy, const G4ParticleDefinition*);

      G4double ComputeEnergy(G4double range, const G4ParticleDefinition*);

      G4double EnergyAfterStep(G4double kinEnergy, G4double step,
                               const G4Material*, const G4ParticleDefinition*);

      G4double EnergyBeforeStep(G4double kinEnergy, G4double step,
                                const G4Material*, const G4ParticleDefinition*);

      G4double TrueStepLength(G4double kinEnergy, G4double step,
                              const G4Material*, const G4ParticleDefinition* part);

      //- inline G4double EnergyAfterStep(G4double kinEnergy, G4double step,
      //-                                const G4Material*, const G4String& particleName);

      //- inline G4double EnergyBeforeStep(G4double kinEnergy, G4double step,
      //-                                  const G4Material*, const G4String& particleName);

      inline G4double AverageScatteringAngle(G4double kinEnergy, G4double step,
                                             const G4Material*,
                                             const G4ParticleDefinition* part);

      //- inline G4double AverageScatteringAngle(G4double kinEnergy, G4double step,
      //-                                        const G4Material*,
      //-                                        const G4String& particleName);

      inline G4double ComputeTrueStep(const G4Material*, const G4ParticleDefinition* part,
                                      G4double kinEnergy, G4double stepLength);

      //- inline G4double EnergyDispersion(G4double kinEnergy, G4double step,
      //-                                  const G4Material*, const G4ParticleDefinition*);

      //- inline G4double EnergyDispersion(G4double kinEnergy, G4double step,
      //-                                  const G4Material*, const G4String& particleName);

      //- inline void SetVerbose(G4int val);

      inline void SetMinKinEnergy(G4double);

      inline void SetMaxKinEnergy(G4double);

      inline void SetMaxEnergyTransfer(G4double);
 
    private:

      void Initialisation();

      G4bool SetupKinematics(const G4ParticleDefinition*, const G4Material*,
                             G4double kinEnergy);

      G4PhysicsTable* PrepareTable();

      const G4ParticleDefinition* FindParticle(const G4String& name);

      void ComputeElectronDEDX(const G4ParticleDefinition* part, G4PhysicsTable* table);

      void ComputeMuonDEDX(const G4ParticleDefinition* part, G4PhysicsTable* table);

      void ComputeHadronDEDX(const G4ParticleDefinition* part, G4PhysicsTable* table);

      void ComputeTransportXS(const G4ParticleDefinition* part, G4PhysicsTable* table);

      inline G4double ComputeValue(G4double x, const G4PhysicsTable* table);

      EnergyLossForExtrapolator& operator=(const EnergyLossForExtrapolator& right);

      EnergyLossForExtrapolator(const EnergyLossForExtrapolator&);

      const G4ParticleDefinition* m_Particle;

      const G4ParticleDefinition* m_Electron;

      const G4ParticleDefinition* m_Positron;

      const G4ParticleDefinition* m_MuonPlus;

      const G4ParticleDefinition* m_MuonMinus;

      const G4ParticleDefinition* m_PionPlus;

      const G4ParticleDefinition* m_PionMinus;

      const G4ParticleDefinition* m_KaonPlus;

      const G4ParticleDefinition* m_KaonMinus;

      const G4ParticleDefinition* m_Proton;

      const G4ParticleDefinition* m_AntiProton;

      const G4ParticleDefinition* m_Deuteron;

      const G4ParticleDefinition* m_AntiDeuteron;

      G4DataVector m_Cuts;

      G4ProductionCuts* m_ProductionCuts;

      std::vector<const G4MaterialCutsCouple*> m_Couples;

      G4PhysicsTable*          m_DedxElectron;

      G4PhysicsTable*          m_DedxPositron;

      G4PhysicsTable*          m_DedxMuon;

      G4PhysicsTable*          m_DedxPion;

      G4PhysicsTable*          m_DedxKaon;

      G4PhysicsTable*          m_DedxProton;

      G4PhysicsTable*          m_DedxDeuteron;

      G4PhysicsTable*          m_RangeElectron;

      G4PhysicsTable*          m_RangePositron;

      G4PhysicsTable*          m_RangeMuon;

      G4PhysicsTable*          m_RangePion;

      G4PhysicsTable*          m_RangeKaon;

      G4PhysicsTable*          m_RangeProton;

      G4PhysicsTable*          m_RangeDeuteron;

      G4PhysicsTable*          m_InvRangeElectron;

      G4PhysicsTable*          m_InvRangePositron;

      G4PhysicsTable*          m_InvRangeMuon;

      G4PhysicsTable*          m_InvRangePion;

      G4PhysicsTable*          m_InvRangeKaon;

      G4PhysicsTable*          m_InvRangeProton;

      G4PhysicsTable*          m_InvRangeDeuteron;

      G4PhysicsTable*          m_MscatElectron;

      const G4Material* m_Material;

      G4int       m_MaterialIndex;

      G4double    m_ElectronDensity;

      G4double    m_RadLength;

      G4double    m_Mass;

      G4double    m_ChargeSq;

      G4double    m_KineticEnergy;

      G4double    m_Gamma;

      G4double    m_BetaGammaSq;

      G4double    m_BetaSq;

      G4double    m_Tmax;

      G4double    m_LinLossLimit;

      G4double    m_UserTmin;

      G4double    m_UserTmax;

      G4double    m_UserMaxEnergyTransfer;

      G4int       m_Nbins;

      G4int       m_NMaterials;

      G4bool      m_Initialised;
    };
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    //- inline G4double EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
    //-                                                            G4double step,
    //-                                                            const G4Material* mat,
    //-                                                            const G4String& name)
    //- {
    //-   return EnergyAfterStep(kinEnergy, step, mat, FindParticle(name));
    //- }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    //- inline G4double EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
    //-                                                             G4double step,
    //-                                                             const G4Material* mat,
    //-                                                             const G4String& name)
    //- {
    //-   return EnergyBeforeStep(kinEnergy, step, mat, FindParticle(name));
    //- }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    //- inline G4double EnergyLossForExtrapolator::AverageScatteringAngle(G4double kinEnergy,
    //-     G4double step,
    //-     const G4Material* mat,
    //-     const G4String& name)
    //- {
    //-   return AverageScatteringAngle(kinEnergy, step, mat, FindParticle(name));
    //- }
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    //- inline G4double
    //- EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
    //-                                             G4double step,
    //-                                             const G4Material* mat,
    //-                                             const G4String& name)
    //- {
    //-   return EnergyDispersion(kinEnergy, step, mat, FindParticle(name));
    //- }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    inline G4double
    EnergyLossForExtrapolator::AverageScatteringAngle(G4double kinEnergy,
                                                      G4double stepLength,
                                                      const G4Material* mat,
                                                      const G4ParticleDefinition* part)
    {
      G4double theta = 0.0;
      if (SetupKinematics(part, mat, kinEnergy)) {
        G4double t = stepLength / m_RadLength;
        G4double y = std::max(0.001, t);
        theta = 19.23 * CLHEP::MeV * std::sqrt(m_ChargeSq * t) * (1.0 + 0.038 * std::log(y)) / (m_BetaSq * m_Gamma * m_Mass);
      }
      return theta;
    }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    inline G4double
    EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat,
                                               const G4ParticleDefinition* part,
                                               G4double kinEnergy,
                                               G4double stepLength)
    {
      G4double theta = AverageScatteringAngle(kinEnergy, stepLength, mat, part);
      return stepLength * std::sqrt(1.0 + 0.625 * theta * theta);
    }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    //- inline
    //- G4double EnergyLossForExtrapolator::EnergyDispersion(G4double kinEnergy,
    //-                                                      G4double stepLength,
    //-                                                      const G4Material* mat,
    //-                                                      const G4ParticleDefinition* part)
    //- {
    //-   G4double sig2 = 0.0;
    //-   if (SetupKinematics(part, mat, kinEnergy)) {
    //-     G4double step = ComputeTrueStep(mat, part, kinEnergy, stepLength);
    //-     sig2 = (1.0 / beta2 - 0.5) * CLHEP::twopi_mc2_rcl2 * tmax * step * electronDensity * charge2;
    //-   }
    //-   return sig2;
    //- }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    inline G4double EnergyLossForExtrapolator::ComputeValue(G4double x,
                                                            const G4PhysicsTable* table)
    {
      G4double res = 0.0;
      G4bool b;
      if (table) res = ((*table)[m_MaterialIndex])->GetValue(x, b);
      return res;
    }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    //- inline void EnergyLossForExtrapolator::SetVerbose(G4int val)
    //- {
    //-   verbose = val;
    //- }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    inline void EnergyLossForExtrapolator::SetMinKinEnergy(G4double val)
    {
      m_UserTmin = val;
    }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    inline void EnergyLossForExtrapolator::SetMaxKinEnergy(G4double val)
    {
      m_UserTmax = val;
    }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
    inline void EnergyLossForExtrapolator::SetMaxEnergyTransfer(G4double val)
    {
      m_UserMaxEnergyTransfer = val;
    }
 
    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
  } //end namespace Simulation

} //end Belle2 namespace