EnergyLossForExtrapolator.h

/**************************************************************************
 * BASF2 (Belle Analysis Framework 2)                                     *
 * Copyright(C) 2010-2011  Belle II Collaboration                         *
 *                                                                        *
 * Author: The Belle II Collaboration                                     *
 * Contributors: Leo Piilonen                                             *
 * Derived from: EnergyLossForExtrapolator.hh                             *
 *               (use geant4e-specific particles; include pions & kaons)  *
 *                                                                        *
 * This software is provided "as is" without any warranty.                *
 **************************************************************************/
 
#ifndef ENERGYLOSSFOREXTRAPOLATOR_H_
#define ENERGYLOSSFOREXTRAPOLATOR_H_
 
#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
 
#endif /* ENERGYLOSSFOREXTRAPOLATOR_H_ */