LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added to SVN Repository

Browse Source
This commit is contained in:
paraiso
2006-02-23 14:12:45 +00:00
parent 2ae8b89fce
commit 10aa47ae02
28 changed files with 5358 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

629
geant4/LEMuSR/G4Modified/G4ChordFinder.cc Normal file
View File

@ -0,0 +1,629 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4ChordFinder.cc,v 1.45 2003/11/13 19:46:56 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// 25.02.97 John Apostolakis, design and implimentation
// 05.03.97 V. Grichine , style modification
// -------------------------------------------------------------------
#include "G4ChordFinder.hh"
#include "G4MagneticField.hh"
#include "G4ElectricField.hh"
#include "G4El_UsualEqRhs.hh"
#include "G4Mag_UsualEqRhs.hh"
#include "G4ClassicalRK4.hh"
#include <iomanip>
#include "G4ios.hh"
#include "G4Transform3D.hh"
#include "G4UnitsTable.hh"
// ..........................................................................
G4ChordFinder::G4ChordFinder(G4MagInt_Driver* pIntegrationDriver)
: fDefaultDeltaChord( 0.25 * mm ),
fDeltaChord( fDefaultDeltaChord ),
fAllocatedStepper(false),
fEquation(0),
fDriversStepper(0),
fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
fMultipleRadius(15.0),
fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
fStatsVerbose(0)
{
// Simple constructor which does not create equation, ..
// fDeltaChord= fDefaultDeltaChord;
fIntgrDriver= pIntegrationDriver;
fAllocatedStepper= false;
fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
// check the values and set the other parameters
}
// ..........................................................................
G4ChordFinder::G4ChordFinder( G4MagneticField* theMagField,
G4double stepMinimum,
G4MagIntegratorStepper* pItsStepper )
: fDefaultDeltaChord( 0.25 * mm ),
fDeltaChord( fDefaultDeltaChord ),
fAllocatedStepper(false),
fEquation(0),
fDriversStepper(0),
fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
fMultipleRadius(15.0),
fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
fStatsVerbose(0)
{
// Construct the Chord Finder
// by creating in inverse order the Driver, the Stepper and EqRhs ...
G4Mag_EqRhs *pEquation = new G4Mag_UsualEqRhs(theMagField);
fEquation = pEquation;
fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
// G4FieldTrack ??
SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
// check the values and set the other parameters
// --->> Charge Q = 0
// --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
if( pItsStepper == 0 )
{
pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
fAllocatedStepper= true;
}
else
{
fAllocatedStepper= false;
}
fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
pItsStepper->GetNumberOfVariables() );
}
// .........................................................................
G4ChordFinder::G4ChordFinder( G4ElectricField* theElField,
G4double stepMinimum,
G4MagIntegratorStepper* pItsStepper )
: fDefaultDeltaChord( 0.25 * mm ),
fDeltaChord( fDefaultDeltaChord ),
fAllocatedStepper(false),
fEquation(0),
fDriversStepper(0),
fFirstFraction(0.999), fFractionLast(1.00), fFractionNextEstimate(0.98),
fMultipleRadius(15.0),
fTotalNoTrials_FNC(0), fNoCalls_FNC(0), fmaxTrials_FNC(0),
fStatsVerbose(0)
{ // Construct the Chord Finder
// by creating in inverse order the Driver, the Stepper and EqRhs ...
G4El_UsualEqRhs *pEquation = new G4El_UsualEqRhs(theElField);
fEquation = pEquation;
fLastStepEstimate_Unconstrained = DBL_MAX; // Should move q, p to
// G4FieldTrack ??
SetFractions_Last_Next( fFractionLast, fFractionNextEstimate);
// check the values and set the other parameters
// --->> Charge Q = 0
// --->> Momentum P = 1 NOMINAL VALUES !!!!!!!!!!!!!!!!!!
if( pItsStepper == 0 )
{
pItsStepper = fDriversStepper = new G4ClassicalRK4(pEquation);
fAllocatedStepper= true;
}
else
{
fAllocatedStepper= false;
}
fIntgrDriver = new G4MagInt_Driver(stepMinimum, pItsStepper,
pItsStepper->GetNumberOfVariables() );
G4cout << " CHORD FINDER :: init..."
<< G4endl;
// tao
mfield=theElField;
}
// ......................................................................
void
G4ChordFinder::SetFractions_Last_Next( G4double fractLast, G4double fractNext )
{
// Use -1.0 as request for Default.
if( fractLast == -1.0 ) fractLast = 1.0; // 0.9;
if( fractNext == -1.0 ) fractNext = 0.98; // 0.9;
// fFirstFraction = 0.999; // Orig 0.999 A safe value, range: ~ 0.95 - 0.999
// fMultipleRadius = 15.0; // For later use, range: ~ 2 - 20
if( fStatsVerbose ) {
G4cout << " ChordFnd> Trying to set fractions: "
<< " first " << fFirstFraction
<< " last " << fractLast
<< " next " << fractNext
<< " and multiple " << fMultipleRadius
<< G4endl;
}
if( (fractLast > 0.0) && (fractLast <=1.0) )
{ fFractionLast= fractLast; }
else
G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
<< " fraction Last = " << fractLast
<< " must be 0 < fractionLast <= 1 " << G4endl;
if( (fractNext > 0.0) && (fractNext <1.0) )
{ fFractionNextEstimate = fractNext; }
else
G4cerr << "G4ChordFinder:: SetFractions_Last_Next: Invalid "
<< " fraction Next = " << fractNext
<< " must be 0 < fractionNext < 1 " << G4endl;
}
// ......................................................................
G4ChordFinder::~G4ChordFinder()
{
delete fEquation; // fIntgrDriver->pIntStepper->theEquation_Rhs;
if( fAllocatedStepper)
{
delete fDriversStepper;
} // fIntgrDriver->pIntStepper;}
delete fIntgrDriver;
if( fStatsVerbose ) { PrintStatistics(); }
}
void
G4ChordFinder::PrintStatistics()
{
// Print Statistics
G4cout << "G4ChordFinder statistics report: " << G4endl;
G4cout
<< " No trials: " << fTotalNoTrials_FNC
<< " No Calls: " << fNoCalls_FNC
<< " Max-trial: " << fmaxTrials_FNC
<< G4endl;
G4cout
<< " Parameters: "
<< " fFirstFraction " << fFirstFraction
<< " fFractionLast " << fFractionLast
<< " fFractionNextEstimate " << fFractionNextEstimate
<< G4endl;
}
// ......................................................................
G4double
G4ChordFinder::AdvanceChordLimited( G4FieldTrack& yCurrent,
G4double stepMax,
G4double epsStep,
const G4ThreeVector latestSafetyOrigin,
G4double latestSafetyRadius
)
{
G4double stepPossible;
G4double dyErr;
G4FieldTrack yEnd( yCurrent);
G4double startCurveLen= yCurrent.GetCurveLength();
G4double nextStep;
// *************
stepPossible= FindNextChord(yCurrent, stepMax, yEnd, dyErr, epsStep, &nextStep
, latestSafetyOrigin, latestSafetyRadius
);
// *************
G4bool good_advance;
if ( dyErr < epsStep * stepPossible )
{
// Accept this accuracy.
yCurrent = yEnd;
good_advance = true;
}
else
{
// Advance more accurately to "end of chord"
// ***************
good_advance = fIntgrDriver->AccurateAdvance(yCurrent, stepPossible, epsStep, nextStep);
// ***************
if ( ! good_advance ){
// In this case the driver could not do the full distance
stepPossible= yCurrent.GetCurveLength()-startCurveLen;
}
}
#ifdef G4DEBUG_FIELD
G4cout << "Exiting FindNextChord Limited with:" << G4endl
<< " yCurrent: " << yCurrent<< G4endl
<<" step possible: " << stepPossible <<G4endl;
G4cout << "G4ChordFinder report: final fieldtrack \n" << yCurrent <<"\n"
<< "spin" << yCurrent.GetSpin() <<"\n"
<<"time of flight" <<yCurrent.GetProperTimeOfFlight()/ns <<"\n"
<<G4endl;
#endif
return stepPossible;
}
// #define TEST_CHORD_PRINT 1
// ............................................................................
G4double
G4ChordFinder::FindNextChord( const G4FieldTrack yStart,
G4double stepMax,
G4FieldTrack& yEnd, // Endpoint
G4double& dyErrPos, // Error of endpoint
G4double epsStep,
G4double* pStepForAccuracy,
const G4ThreeVector, // latestSafetyOrigin,
G4double // latestSafetyRadius
)
// Returns Length of Step taken
{
// G4int stepRKnumber=0;
G4FieldTrack yCurrent= yStart;
G4double stepTrial, stepForAccuracy;
G4double dydx[G4FieldTrack::ncompSVEC];
// 1.) Try to "leap" to end of interval
// 2.) Evaluate if resulting chord gives d_chord that is good enough.
// 2a.) If d_chord is not good enough, find one that is.
G4bool validEndPoint= false;
G4double dChordStep, lastStepLength; // stepOfLastGoodChord;
#ifdef DEBUG_FIELD
G4cout <<"\n---------- G4ChordFinder :: init FieldTrack : " << yCurrent << "\n------------";
#endif
fIntgrDriver-> GetDerivatives( yCurrent, dydx ) ;
#ifdef DEBUG_FIELD
G4cout <<"\n---------- G4ChordFinder :: derivative got, entering loop : dydx "
<< dydx[0] <<" " << dydx[1] <<" "<< dydx[2] <<" "<<"\n"
<< dydx[3] <<" " << dydx[4] <<" "<< dydx[5] <<" "<<"\n"
<< dydx[6] <<" " << dydx[7] <<" "<< dydx[8] <<" "<<"\n"
<< dydx[9] <<" " << dydx[10] <<" "<< dydx[11] <<" "<<"\n"
<<G4endl;
#endif
G4int noTrials=0;
const G4double safetyFactor= fFirstFraction; // 0.975 or 0.99 ? was 0.999
stepTrial = std::min( stepMax,
safetyFactor * fLastStepEstimate_Unconstrained );
G4double newStepEst_Uncons= 0.0;
do
{
#ifdef DEBUG_FIELD
G4cout <<"\n---------- G4ChordFinder :: quick Advance : " << yCurrent << "\n------------";
#endif
G4double stepForChord;
yCurrent = yStart; // Always start from initial point
// ************
fIntgrDriver->QuickAdvance( yCurrent, dydx, stepTrial,
dChordStep, dyErrPos);
// ************
#ifdef DEBUG_FIELD
G4cout <<"\n---------- G4ChordFinder :: current FieldTrack : " << yCurrent << "\n------------";
#endif
// We check whether the criterion is met here.
validEndPoint = AcceptableMissDist(dChordStep);
// && (dyErrPos < eps) ;
lastStepLength = stepTrial;
// This method estimates to step size for a good chord.
stepForChord = NewStep(stepTrial, dChordStep, newStepEst_Uncons );
if( ! validEndPoint ) {
if( stepTrial<=0.0 )
stepTrial = stepForChord;
else if (stepForChord <= stepTrial)
// Reduce by a fraction, possibly up to 20%
stepTrial = std::min( stepForChord,
fFractionLast * stepTrial);
else
stepTrial *= 0.1;
// if(dbg) G4cerr<<"Dchord too big. Try new hstep="<<stepTrial<<G4endl;
}
// #ifdef TEST_CHORD_PRINT
// TestChordPrint( noTrials, lastStepLength, dChordStep, stepTrial );
// #endif
noTrials++;
}
while( ! validEndPoint ); // End of do-while RKD
if( newStepEst_Uncons > 0.0 ){
fLastStepEstimate_Unconstrained= newStepEst_Uncons;
}
AccumulateStatistics( noTrials );
// stepOfLastGoodChord = stepTrial;
if( pStepForAccuracy ){
// Calculate the step size required for accuracy, if it is needed
G4double dyErr_relative = dyErrPos/(epsStep*lastStepLength);
if( dyErr_relative > 1.0 ) {
stepForAccuracy =
fIntgrDriver->ComputeNewStepSize( dyErr_relative,
lastStepLength );
}else{
stepForAccuracy = 0.0; // Convention to show step was ok
}
*pStepForAccuracy = stepForAccuracy;
}
#ifdef TEST_CHORD_PRINT
// static int dbg=0;
// if( dbg )
G4cout << "ChordF/FindNextChord: NoTrials= " << noTrials
<< " StepForGoodChord=" << std::setw(10) << stepTrial << G4endl;
#endif
yEnd= yCurrent;
return stepTrial;
}
// ----------------------------------------------------------------------------
#if 0
// #ifdef G4VERBOSE
if( dbg ) {
G4cerr << "Returned from QuickAdvance with: yCur=" << yCurrent <<G4endl;
G4cerr << " dChordStep= "<< dChordStep <<" dyErr=" << dyErr << G4endl;
}
#endif
// ----------------------------------------------------------------------------
// ...........................................................................
G4double G4ChordFinder::NewStep(G4double stepTrialOld,
G4double dChordStep, // Curr. dchord achieved
G4double& stepEstimate_Unconstrained )
//
// Is called to estimate the next step size, even for successful steps,
// in order to predict an accurate 'chord-sensitive' first step
// which is likely to assist in more performant 'stepping'.
//
{
G4double stepTrial;
static G4double lastStepTrial = 1., lastDchordStep= 1.;
#if 1
// const G4double threshold = 1.21, multiplier = 0.9;
// 0.9 < 1 / sqrt(1.21)
if (dChordStep > 0.0)
{
stepEstimate_Unconstrained = stepTrialOld*sqrt( fDeltaChord / dChordStep );
// stepTrial = 0.98 * stepEstimate_Unconstrained;
stepTrial = fFractionNextEstimate * stepEstimate_Unconstrained;
}
else
{
// Should not update the Unconstrained Step estimate: incorrect!
stepTrial = stepTrialOld * 2.;
}
// if ( dChordStep < threshold * fDeltaChord ){
// stepTrial= stepTrialOld * multiplier;
// }
if( stepTrial <= 0.001 * stepTrialOld)
{
if ( dChordStep > 1000.0 * fDeltaChord ){
stepTrial= stepTrialOld * 0.03;
}else{
if ( dChordStep > 100. * fDeltaChord ){
stepTrial= stepTrialOld * 0.1;
}else{
// Try halving the length until dChordStep OK
stepTrial= stepTrialOld * 0.5;
}
}
}else if (stepTrial > 1000.0 * stepTrialOld)
{
stepTrial= 1000.0 * stepTrialOld;
}
if( stepTrial == 0.0 ){
stepTrial= 0.000001;
}
lastStepTrial = stepTrialOld;
lastDchordStep= dChordStep;
#else
if ( dChordStep > 1000. * fDeltaChord ){
stepTrial= stepTrialOld * 0.03;
}else{
if ( dChordStep > 100. * fDeltaChord ){
stepTrial= stepTrialOld * 0.1;
}else{
// Keep halving the length until dChordStep OK
stepTrial= stepTrialOld * 0.5;
}
}
#endif
// A more sophisticated chord-finder could figure out a better
// stepTrial, from dChordStep and the required d_geometry
// eg
// Calculate R, r_helix (eg at orig point)
// if( stepTrial < 2 pi R )
// stepTrial = R arc_cos( 1 - fDeltaChord / r_helix )
// else
// ??
return stepTrial;
}
//
// Given a starting curve point A (CurveA_PointVelocity), a later
// curve point B (CurveB_PointVelocity) and a point E which is (generally)
// not on the curve, find and return a point F which is on the curve and
// which is close to E. While advancing towards F utilise eps_step
// as a measure of the relative accuracy of each Step.
G4FieldTrack
G4ChordFinder::ApproxCurvePointV( const G4FieldTrack& CurveA_PointVelocity,
const G4FieldTrack& CurveB_PointVelocity,
const G4ThreeVector& CurrentE_Point,
G4double eps_step)
{
// 1st implementation:
// if r=|AE|/|AB|, and s=true path lenght (AB)
// return the point that is r*s along the curve!
G4FieldTrack Current_PointVelocity= CurveA_PointVelocity;
G4ThreeVector CurveA_Point= CurveA_PointVelocity.GetPosition();
G4ThreeVector CurveB_Point= CurveB_PointVelocity.GetPosition();
G4ThreeVector ChordAB_Vector= CurveB_Point - CurveA_Point;
G4ThreeVector ChordAE_Vector= CurrentE_Point - CurveA_Point;
G4double ABdist= ChordAB_Vector.mag();
G4double curve_length; // A curve length of AB
G4double AE_fraction;
curve_length= CurveB_PointVelocity.GetCurveLength()
- CurveA_PointVelocity.GetCurveLength();
// const
G4double integrationInaccuracyLimit= std::max( perMillion, 0.5*eps_step );
if( curve_length < ABdist * (1. - integrationInaccuracyLimit) ){
#ifdef G4DEBUG_FIELD
G4cerr << " Warning in G4ChordFinder::ApproxCurvePoint: "
<< G4endl
<< " The two points are further apart than the curve length "
<< G4endl
<< " Dist = " << ABdist
<< " curve length = " << curve_length
<< " relativeDiff = " << (curve_length-ABdist)/ABdist
<< G4endl;
if( curve_length < ABdist * (1. - 10*eps_step) ) {
G4cerr << " ERROR: the size of the above difference"
<< " exceeds allowed limits. Aborting." << G4endl;
G4Exception("G4ChordFinder::ApproxCurvePointV()", "PrecisionError",
FatalException, "Unphysical curve length.");
}
#endif
// Take default corrective action:
// --> adjust the maximum curve length.
// NOTE: this case only happens for relatively straight paths.
curve_length = ABdist;
}
G4double new_st_length;
if ( ABdist > 0.0 ){
AE_fraction = ChordAE_Vector.mag() / ABdist;
}else{
AE_fraction = 0.5; // Guess .. ?;
#ifdef G4DEBUG_FIELD
G4cout << "Warning in G4ChordFinder::ApproxCurvePoint:"
<< " A and B are the same point!" << G4endl
<< " Chord AB length = " << ChordAE_Vector.mag() << G4endl
<< G4endl;
#endif
}
if( (AE_fraction> 1.0 + perMillion) || (AE_fraction< 0.) ){
#ifdef G4DEBUG_FIELD
G4cerr << " G4ChordFinder::ApproxCurvePointV - Warning:"
<< " Anomalous condition:AE > AB or AE/AB <= 0 " << G4endl
<< " AE_fraction = " << AE_fraction << G4endl
<< " Chord AE length = " << ChordAE_Vector.mag() << G4endl
<< " Chord AB length = " << ABdist << G4endl << G4endl;
G4cerr << " OK if this condition occurs after a recalculation of 'B'"
<< G4endl << " Otherwise it is an error. " << G4endl ;
#endif
// This course can now result if B has been re-evaluated,
// without E being recomputed (1 July 99)
// In this case this is not a "real error" - but it undesired
// and we cope with it by a default corrective action ...
AE_fraction = 0.5; // Default value
}
new_st_length= AE_fraction * curve_length;
G4bool good_advance;
if ( AE_fraction > 0.0 ) {
good_advance =
fIntgrDriver->AccurateAdvance(Current_PointVelocity,
new_st_length,
eps_step ); // Relative accuracy
// In this case it does not matter if it cannot advance the full distance
}
// If there was a memory of the step_length actually require at the start
// of the integration Step, this could be re-used ...
return Current_PointVelocity;
}
void
G4ChordFinder::TestChordPrint( G4int noTrials,
G4int lastStepTrial,
G4double dChordStep,
G4double nextStepTrial )
{
G4int oldprec= G4cout.precision(5);
G4cout << " ChF/fnc: notrial " << std::setw( 3) << noTrials
<< " this_step= " << std::setw(10) << lastStepTrial;
if( fabs( (dChordStep / fDeltaChord) - 1.0 ) < 0.001 ){
G4cout.precision(8);
}else{ G4cout.precision(6); }
G4cout << " dChordStep= " << std::setw(12) << dChordStep;
if( dChordStep > fDeltaChord ) { G4cout << " d+"; }
else { G4cout << " d-"; }
G4cout.precision(5);
G4cout << " new_step= " << std::setw(10)
<< fLastStepEstimate_Unconstrained
<< " new_step_constr= " << std::setw(10)
<< lastStepTrial << G4endl;
G4cout << " nextStepTrial = " << std::setw(10) << nextStepTrial << G4endl;
G4cout.precision(oldprec);
}

198
geant4/LEMuSR/G4Modified/G4ChordFinder.hh Normal file
View File

@ -0,0 +1,198 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4ChordFinder.hh,v 1.16 2003/11/13 17:53:47 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// class G4ChordFinder
//
// Class description:
//
// A class that provides RK integration of motion ODE (as does g4magtr)
// and also has a method that returns an Approximate point on the curve
// near to a (chord) point.
// History:
// - 25.02.97 John Apostolakis, design and implementation
// - 05.03.97 V. Grichine , makeup to G4 'standard'
// -------------------------------------------------------------------
#ifndef G4CHORDFINDER_HH
#define G4CHORDFINDER_HH
#include "G4MagIntegratorDriver.hh"
#include "G4FieldTrack.hh"
class G4MagneticField;
class G4ElectricField;
class G4ChordFinder
{
public: // with description
G4ChordFinder( G4MagInt_Driver* pIntegrationDriver );
G4ChordFinder( G4MagneticField* itsMagField,
G4double stepMinimum = 1.0e-2 * mm,
G4MagIntegratorStepper* pItsStepper = 0 );
// A constructor that creates defaults for all "children" classes.
G4ChordFinder( G4ElectricField* itsElField,
G4double stepMinimum = 1.0e-2 * mm,
G4MagIntegratorStepper* pItsStepper = 0 );
virtual ~G4ChordFinder();
G4double AdvanceChordLimited( G4FieldTrack& yCurrent,
G4double stepInitial,
G4double epsStep_Relative,
const G4ThreeVector latestSafetyOrigin,
G4double lasestSafetyRadius);
// Uses ODE solver's driver to find the endpoint that satisfies
// the chord criterion: that d_chord < delta_chord
// -> Returns Length of Step taken.
G4FieldTrack ApproxCurvePointV(const G4FieldTrack& curveAPointVelocity,
const G4FieldTrack& curveBPointVelocity,
const G4ThreeVector& currentEPoint,
G4double epsStep);
inline G4double GetDeltaChord() const;
inline void SetDeltaChord(G4double newval);
inline void SetChargeMomentumMass(G4double pCharge, // in e+ units
G4double pMomentum,
G4double pMass );
// Function to inform integration driver of charge, speed.
inline void SetIntegrationDriver(G4MagInt_Driver* IntegrationDriver);
inline G4MagInt_Driver* GetIntegrationDriver();
// Access and set Driver.
inline void ResetStepEstimate();
// Clear internal state (last step estimate)
inline G4int GetNoCalls();
inline G4int GetNoTrials(); // Total number of trials
inline G4int GetNoMaxTrials(); // Maximum # of trials for one call
// Get statistics about number of calls & trials in FindNextChord
virtual void PrintStatistics();
// A report with the above -- and possibly other stats
inline G4int SetVerbose( G4int newvalue=1);
// Set verbosity and return old value
protected: // .........................................................
inline void AccumulateStatistics( G4int noTrials );
// Accumulate the basic statistics
// - other specialised ones must be kept by derived classes
inline G4bool AcceptableMissDist(G4double dChordStep) const;
G4double NewStep( G4double stepTrialOld,
G4double dChordStep, // Current dchord estimate
G4double& stepEstimate_Unconstrained ) ;
virtual G4double FindNextChord( const G4FieldTrack yStart,
G4double stepMax,
G4FieldTrack& yEnd,
G4double& dyErr, // Error of endpoint
G4double epsStep,
G4double* pNextStepForAccuracy, // = 0,
const G4ThreeVector latestSafetyOrigin,
G4double latestSafetyRadius
);
void PrintDchordTrial(G4int noTrials,
G4double stepTrial,
G4double oldStepTrial,
G4double dChordStep);
public: // no description
void TestChordPrint( G4int noTrials,
G4int lastStepTrial,
G4double dChordStep,
G4double nextStepTrial );
// Printing for monitoring ...
inline G4double GetFirstFraction(); // Originally 0.999
inline G4double GetFractionLast(); // Originally 1.000
inline G4double GetFractionNextEstimate(); // Originally 0.980
inline G4double GetMultipleRadius(); // No original value
// Parameters for adapting performance ... use with great care
public: // with description
void SetFractions_Last_Next( G4double fractLast= 0.90,
G4double fractNext= 0.95 );
// Parameters for performance ... change with great care
inline void SetFirstFraction(G4double fractFirst);
// Parameter for performance ... change with great care
protected:
inline G4double GetLastStepEstimateUnc();
inline void SetLastStepEstimateUnc( G4double stepEst );
private: // ............................................................
G4ChordFinder(const G4ChordFinder&);
G4ChordFinder& operator=(const G4ChordFinder&);
// Private copy constructor and assignment operator.
private: // ............................................................
// G4int nOK, nBAD;
G4MagInt_Driver* fIntgrDriver;
const G4double fDefaultDeltaChord; // SET in G4ChordFinder.cc = 0.25 mm
G4double fDeltaChord; // Maximum miss distance
G4double fLastStepEstimate_Unconstrained; // State information for efficiency
// Variables used in construction/destruction
G4bool fAllocatedStepper;
G4EquationOfMotion* fEquation;
G4MagIntegratorStepper* fDriversStepper;
// Parameters
G4double fFirstFraction, fFractionLast, fFractionNextEstimate;
G4double fMultipleRadius;
// For Statistics
// -- G4int fNoTrials, fNoCalls;
G4int fTotalNoTrials_FNC, fNoCalls_FNC, fmaxTrials_FNC; // fnoTimesMaxTrFNC;
G4int fStatsVerbose; // if > 0, print Statistics in destructor
G4double itime, ftime;
G4ElectricField* mfield;
G4int firstChord;
};
// Inline function implementation:
#include "G4ChordFinder.icc"
#endif // G4CHORDFINDER_HH

80
geant4/LEMuSR/G4Modified/G4El_EqRhs.cc Normal file
View File

@ -0,0 +1,80 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4El_EqRhs.cc,v 1.10 2003/11/05 16:33:55 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
// This is the standard right-hand side for equation of motion
// in a pure Electric Field .
//
// Other that might be required are:
// i) is when using a moving reference frame ... or
// ii) extending for other forces, eg an electric field
//
// J. Apostolakis, January 13th, 1997
//
// --------------------------------------------------------------------
#include "G4ElectricField.hh"
#include "G4El_EqRhs.hh"
#include "globals.hh"
//const G4double G4El_EqRhs::fUnitConstant = 0.299792458 * (GeV/(tesla*m));
// Constructor Implementation
//
G4El_EqRhs::G4El_EqRhs( G4ElectricField *elField )
: G4EquationOfMotion(elField)
{
}
G4El_EqRhs::~G4El_EqRhs() { }
void G4El_EqRhs::RightHandSide( const G4double y[],
G4double dydx[] ) const
{
G4double Field[3];
G4double PositionAndTime[4];
// Position
PositionAndTime[0] = y[0];
PositionAndTime[1] = y[1];
PositionAndTime[2] = y[2];
// Global Time
PositionAndTime[3] = y[7]; // See G4FieldTrack::LoadFromArray
G4cout <<"EL_EQ RIGHT HAND SIDE!"<<G4endl;
GetFieldValue(PositionAndTime, Field) ;
EvaluateRhsGivenB( y, Field, dydx );
}
void
G4El_EqRhs::SetChargeMomentumMass( G4double particleCharge, // e+ units
G4double , // MomentumXc
G4double ) // particleMass
{
fCof_val = particleCharge*eplus*c_light ;
}

87
geant4/LEMuSR/G4Modified/G4El_EqRhs.hh Normal file
View File

@ -0,0 +1,87 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4El_EqRhs.hh,v 1.8 2003/11/05 12:54:13 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// class G4El_EqRhs
//
// Class description:
//
// The "standard" equation of motion of a particle in a pure electric field.
// History:
// - Created. J.Apostolakis, January 13th 1997
// --------------------------------------------------------------------
#ifndef G4_EL_EQRHS_DEF
#define G4_EL_EQRHS_DEF
#include "G4Types.hh"
#include "G4EquationOfMotion.hh"
class G4ElectricField;
class G4El_EqRhs : public G4EquationOfMotion
{
public: // with description
G4El_EqRhs( G4ElectricField *elField );
virtual ~G4El_EqRhs();
// Constructor and destructor. No actions.
virtual void EvaluateRhsGivenB( const G4double y[],
const G4double E[3],
G4double dydx[] ) const = 0;
// Given the value of the field "B", this function
// calculates the value of the derivative dydx.
// This is the _only_ function a subclass must define.
// The other two functions use Rhs_givenB.
inline G4double FCof() const;
void RightHandSide( const G4double y[], G4double dydx[] ) const;
virtual void SetChargeMomentumMass( G4double particleCharge, // in e+ units
G4double MomentumXc,
G4double mass);
private:
G4double fCof_val;
static const G4double fUnitConstant; // Set in G4El_EqRhs.cc
// to 0.299792458
// Coefficient in the Lorentz motion equation (Lorentz force), if the
// electric field B is in Tesla, the particle charge in units of the
// elementary (positron?) charge, the momentum P in MeV/c, and the
// space coordinates and path along the trajectory in mm .
};
inline
G4double G4El_EqRhs::FCof() const
{
return fCof_val;
}
#endif /* G4_EL_EQRHS_DEF */

76
geant4/LEMuSR/G4Modified/G4El_MagEqRhs.cc Normal file
View File

@ -0,0 +1,76 @@
#include "G4El_MagEqRhs.hh"
#include "G4Mag_SpinEqRhs.hh"
#include "G4MagneticField.hh"
#include "G4ThreeVector.hh"
#include "G4El_UsualEqRhs.hh"
#include "G4ElectricField.hh"
#include "G4ios.hh"
G4El_MagEqRhs::G4El_MagEqRhs( G4Mag_EqRhs *Meq, G4El_EqRhs *Eeq,G4Field* field)
:G4EquationOfMotion(field)
{
fMagEq=Meq;
fElEq=Eeq;
}
G4El_MagEqRhs::~G4El_MagEqRhs()
{;}
void G4El_MagEqRhs::RightHandSide( const G4double y[],
G4double dydx[] )const
{
G4double MagField[3];
G4double ElField[3];
G4double dydx1[12];
G4double dydx2[12];
G4double PositionAndTime[4];
// Position
PositionAndTime[0] = y[0];
PositionAndTime[1] = y[1];
PositionAndTime[2] = y[2];
// Global Time
PositionAndTime[3] = y[7];
// Get Respective Field Values
fMagEq->GetFieldValue(PositionAndTime, MagField) ;
fElEq->GetFieldValue(PositionAndTime, ElField) ;
fMagEq->EvaluateRhsGivenB( y, MagField, &dydx1[0] );
fElEq->EvaluateRhsGivenB( y, ElField, &dydx2[0] );
G4int i;
i=0;
for(i=0;i==18;i++)
{
dydx[i] = dydx1[i]+dydx2[i];
}
}
void G4El_MagEqRhs::SetChargeMomentumMass( G4double particleCharge, // e+ units
G4double MomentumXc, // MomentumXc
G4double mass) // particleMass
{
fMagEq->SetChargeMomentumMass( particleCharge, MomentumXc , mass);
fElEq->SetChargeMomentumMass( particleCharge,MomentumXc , mass);
}
void G4El_MagEqRhs::EvaluateRhsGivenB( const G4double y[],
const G4double B[3],
G4double dydx[] ) const
{}

42
geant4/LEMuSR/G4Modified/G4El_MagEqRhs.hh Normal file
View File

@ -0,0 +1,42 @@
#ifndef G4EL_MAGEQRHS
#define G4EL_MAGEQRHS
#include "G4Types.hh"
#include "G4Mag_EqRhs.hh"
#include "G4El_UsualEqRhs.hh"
#include "G4ios.hh"
#include "G4EquationOfMotion.hh"
class G4MagneticField;
class G4ElectricField;
class G4Mag_EqRhs;
class G4El_UsualEqRhs;
class G4El_MagEqRhs : public G4EquationOfMotion
{
public: // with description
G4El_MagEqRhs( G4Mag_EqRhs* , G4El_EqRhs*,G4Field* field );
~G4El_MagEqRhs();
// Constructor and destructor. No actions.
void RightHandSide( const G4double y[], G4double dydx[] ) const;
virtual void SetChargeMomentumMass( G4double particleCharge, // in e+ units
G4double MomentumXc,
G4double mass);
void EvaluateRhsGivenB( const G4double y[],
const G4double E[3],
G4double dydx[] ) const;
private:
G4Mag_EqRhs *fMagEq;
G4El_EqRhs *fElEq;
};
#endif /* G4EL_MAGEQRHS */

110
geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.cc Normal file
View File

@ -0,0 +1,110 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4El_UsualEqRhs.cc,v 1.10 2003/11/05 17:31:31 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// This is the 'standard' right-hand side for the equation of motion
// of a charged particle in a magnetic field.
//
// Initial version: J. Apostolakis, January 13th, 1997
//
// --------------------------------------------------------------------
#include "G4UnitsTable.hh"
#include "G4El_UsualEqRhs.hh"
#include "G4ElectricField.hh"
#include "G4ios.hh"
G4El_UsualEqRhs::G4El_UsualEqRhs( G4ElectricField* ElField )
: G4El_EqRhs( ElField ) {}
G4El_UsualEqRhs::~G4El_UsualEqRhs() {}
void
G4El_UsualEqRhs::EvaluateRhsGivenB( const G4double y[],
const G4double E[3],
G4double dydx[] ) const
{
G4double momentum_square = y[3]*y[3] + y[4]*y[4] + y[5]*y[5];
G4double inv_momentum_magnitude = 1.0 / sqrt( momentum_square );
G4double cof = cst*inv_momentum_magnitude;
dydx[0] = y[3]*inv_momentum_magnitude; // (d/ds)x = Vx/V
dydx[1] = y[4]*inv_momentum_magnitude; // (d/ds)y = Vy/V
dydx[2] = y[5]*inv_momentum_magnitude; // (d/ds)z = Vz/V
dydx[3] = cof*(E[0]) ; // Ax = a*(Ex)
dydx[4] = cof*(E[1]) ; // Ay = a*(Ey)
dydx[5] = cof*(E[2]) ; // Az = a*(Ez)
#ifdef DEBUG_FIELD
G4cout<<"LEMuSREl_UsualEqRhs :: posmomE \n"
<< y[0]/100 <<" " << y[1]/100 <<" "<< y[2]/100+5.67 <<" "<<"\n"
<< y[3] <<" " << y[4] <<" "<< y[5] <<" "<<"\n"
<< E[0]/volt*meter <<" " << E[1]/volt*meter <<" "<< E[2]/volt*meter <<" "<<"\n"
<<G4endl;
G4cout<<"LEMuSREl_UsualEqRhs :: dydx \n"
<< dydx[0] <<" " << dydx[1] <<" "<< dydx[2] <<" "<<"\n"
<< dydx[3] <<" " << dydx[4] <<" "<< dydx[5] <<" "<<"\n"
<< dydx[6] <<" " << dydx[7] <<" "<< dydx[8] <<" "<<"\n"
<< dydx[9] <<" " << dydx[10] <<" "<< dydx[11] <<" "<<"\n"
<<G4endl;
// getchar();
#endif
return ;
}
void
G4El_UsualEqRhs:: SetChargeMomentumMass( G4double particleCharge, // in e+ units
G4double MomentumXc,
G4double mass )//mass
{
fInvCurrentMomentumXc= 1.0 / MomentumXc;
cst = particleCharge*eplus*mass;//*c_light;
}
void G4El_UsualEqRhs::RightHandSide( const G4double y[],
G4double dydx[] ) const
{
G4double Field[3];
G4double PositionAndTime[4];
// Position
PositionAndTime[0] = y[0];
PositionAndTime[1] = y[1];
PositionAndTime[2] = y[2];
// Global Time
PositionAndTime[3] = y[7]; // See G4FieldTrack::LoadFromArray
#ifdef DEBUG_FIELD
G4cout <<"EL_USUALEQ RIGHT HAND SIDE!"<<G4endl;
#endif
GetFieldValue(PositionAndTime, Field) ;
EvaluateRhsGivenB( y, Field, dydx );
}

76
geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.hh Normal file
View File

@ -0,0 +1,76 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4El_UsualEqRhs.hh,v 1.6 2003/10/31 14:35:52 gcosmo Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// class G4El_UsualEqRhs
//
// Class description:
//
// This is the standard right-hand side for equation of motion.
// The only case another is required is when using a moving reference
// frame ... or extending the class to include additional Forces,
// eg an electric field
// History:
// - Created: J. Apostolakis, January 13th 1997.
// --------------------------------------------------------------------
#ifndef G4EL_USUAL_EQRHS
#define G4EL_USUAL_EQRHS
#include "G4El_EqRhs.hh"
class G4ElectricField;
class G4El_UsualEqRhs : public G4El_EqRhs
{
public: // with description
G4El_UsualEqRhs( G4ElectricField* ElField );
~G4El_UsualEqRhs();
// Constructor and destructor. No actions.
void EvaluateRhsGivenB( const G4double y[],
const G4double E[3],
G4double dydx[] ) const;
// Given the value of the electric field E, this function
// calculates the value of the derivative dydx.
virtual void SetChargeMomentumMass( G4double particleCharge, // in e+ units
G4double MomentumXc,
G4double mass);
void RightHandSide( const G4double y[], G4double dydx[] ) const;
private:
G4double cst;
G4double fInvCurrentMomentumXc; // This extra state enables us
// to save a square root in a
// critical method.
};
#endif /* G4EL_USUAL_EQRHS */

119
geant4/LEMuSR/G4Modified/G4FieldManager.cc Normal file
View File

@ -0,0 +1,119 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4FieldManager.cc,v 1.13 2003/11/08 04:08:13 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
// - Augsut 05 T.K.Paraïso add magcomponent l.138
// -------------------------------------------------------------------
#include "G4FieldManager.hh"
#include "G4Field.hh"
#include "G4MagneticField.hh"
#include "G4ChordFinder.hh"
G4FieldManager::G4FieldManager(G4Field *detectorField,
G4ChordFinder *pChordFinder,
G4bool fieldChangesEnergy
)
: fDetectorField(detectorField),
fChordFinder(pChordFinder),
fAllocatedChordFinder(false),
fDefault_Delta_One_Step_Value(0.01*mm),
fDefault_Delta_Intersection_Val(0.001*mm),
fEpsilonMinDefault(5.0e-5),
fEpsilonMaxDefault(0.001),
fEpsilonMin( fEpsilonMinDefault ),
fEpsilonMax( fEpsilonMaxDefault),
fMagComponent(false)
{
fDelta_One_Step_Value= fDefault_Delta_One_Step_Value;
fDelta_Intersection_Val= fDefault_Delta_Intersection_Val;
if ( detectorField )
fFieldChangesEnergy= detectorField->DoesFieldChangeEnergy();
else
fFieldChangesEnergy= fieldChangesEnergy;
}
G4FieldManager::G4FieldManager(G4MagneticField *detectorField)
: fDetectorField(detectorField), fAllocatedChordFinder(true),
fFieldChangesEnergy(false),
fDefault_Delta_One_Step_Value(0.01*mm),
fDefault_Delta_Intersection_Val(0.001*mm),
fEpsilonMinDefault(5.0e-5),
fEpsilonMaxDefault(0.001),
fEpsilonMin( fEpsilonMinDefault ),
fEpsilonMax( fEpsilonMaxDefault),
fMagComponent(true)
{
fChordFinder= new G4ChordFinder( detectorField );
fDelta_One_Step_Value= fDefault_Delta_One_Step_Value;
fDelta_Intersection_Val= fDefault_Delta_Intersection_Val;
}
void G4FieldManager::ConfigureForTrack( const G4Track * )
{
// Default is to do nothing!
;
}
G4FieldManager::~G4FieldManager()
{
if( fAllocatedChordFinder ){
delete fChordFinder;
}
}
void
G4FieldManager::CreateChordFinder(G4MagneticField *detectorMagField)
{
if ( fAllocatedChordFinder )
delete fChordFinder;
fChordFinder= new G4ChordFinder( detectorMagField );
fAllocatedChordFinder= true;
}
G4bool G4FieldManager::SetDetectorField(G4Field *pDetectorField)
{
fDetectorField= pDetectorField;
if ( pDetectorField )
fFieldChangesEnergy= pDetectorField->DoesFieldChangeEnergy();
else
fFieldChangesEnergy= false; // No field
return false;
}
G4bool G4FieldManager::SetDetectorField(G4Field *pDetectorField,G4bool mag)
{
fDetectorField= pDetectorField;
fMagComponent = mag;
if ( pDetectorField )
fFieldChangesEnergy= pDetectorField->DoesFieldChangeEnergy();
else
fFieldChangesEnergy= false; // No field
return false;
}

191
geant4/LEMuSR/G4Modified/G4FieldManager.hh Normal file
View File

@ -0,0 +1,191 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4FieldManager.hh,v 1.13 2003/11/08 03:55:16 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// class G4FieldManager
//
// Class description:
//
// A class to manage (Store) a pointer to the Field subclass that
// describes the field of a detector (magnetic, electric or other).
// Also stores a reference to the chord finder.
//
// The G4FieldManager class exists to allow the user program to specify
// the electric, magnetic and/or other field(s) of the detector.
// (OR, in the future, of a part of it - planned to be a logical volume).
// It also stores a pointer to the ChordFinder object that can do the
// propagation in this field. All geometrical track "advancement"
// in the field is handled by this ChordFinder object.
//
// G4FieldManager allows the other classes/object (of the MagneticField
// & other class categories) to find out whether a detector field object
// exists and what that object is.
//
// The Chord Finder must be created either by calling CreateChordFinder
// for a Magnetic Field or by the user creating a a Chord Finder object
// "manually" and setting this pointer.
//
// A default FieldManager is created by the singleton class
// G4NavigatorForTracking and exists before main is called.
// However a new one can be created and given to G4NavigatorForTracking.
//
// Our current design envisions that one Field manager is
// valid for each region detector.
// History:
// - 10.03.97 John Apostolakis, design and implementation.
// - Augsut 05 T.K.Paraïso add magcomponent l.138
// -------------------------------------------------------------------
#ifndef G4FIELDMANAGER_HH
#define G4FIELDMANAGER_HH 1
#include "globals.hh"
class G4Field;
class G4MagneticField;
class G4ElectricField;
class G4ChordFinder;
class G4Track; // Forward reference for parameter configuration
class G4FieldManager
{
public: // with description
G4FieldManager(G4Field *detectorField=0,
G4ChordFinder *pChordFinder=0,
G4bool b=true ); // fieldChangesEnergy is taken from field
// General constructor for any field.
// -> Must be set with field and chordfinder for use.
G4FieldManager(G4MagneticField *detectorMagneticField);
// Creates ChordFinder
// - assumes pure magnetic field (so Energy constant)
virtual ~G4FieldManager();
G4bool SetDetectorField(G4Field *detectorField);
// TAO
G4bool SetDetectorField(G4Field *detectorField, G4bool magcomponent);
inline const G4Field* GetDetectorField() const;
inline G4bool DoesFieldExist() const;
// Set, get and check the field object
void CreateChordFinder(G4MagneticField *detectorMagField);
inline void SetChordFinder(G4ChordFinder *aChordFinder);
inline G4ChordFinder* GetChordFinder();
inline const G4ChordFinder* GetChordFinder() const;
// Create, set or get the associated Chord Finder
virtual void ConfigureForTrack( const G4Track * );
// Setup the choice of the configurable parameters
// relying on the current track's energy, particle identity, ..
// Note: In addition to the values of member variables,
// a user can use this to change the ChordFinder, the field, ...
public: // with description
inline G4double GetDeltaIntersection() const; // virtual ?
// Accuracy for boundary intersection.
inline G4double GetDeltaOneStep() const; // virtual ?
// Accuracy for one tracking/physics step.
inline void SetAccuraciesWithDeltaOneStep(G4double valDeltaOneStep);
// Sets both accuracies, maintaining a fixed ratio for accuracties
// of volume Intersection and Integration (in One Step)
inline void SetDeltaOneStep(G4double valueD1step);
// Set accuracy for integration of one step. (only)
inline void SetDeltaIntersection(G4double valueDintersection);
// Set accuracy of intersection of a volume. (only)
inline G4double GetMinimumEpsilonStep() const;
inline void SetMinimumEpsilonStep( G4double newEpsMin );
// Minimum for Relative accuracy of a Step
inline G4double GetMaximumEpsilonStep() const;
inline void SetMaximumEpsilonStep( G4double newEpsMax );
// Maximum for Relative accuracy of a Step
inline G4bool DoesFieldChangeEnergy() const;
inline void SetFieldChangesEnergy(G4bool value);
// For electric field this should be true
// For electromagnetic field this should be true
// For magnetic field this should be false
inline G4bool FieldHasMagComponent() const;
inline void SetFieldMagComponent(G4bool value);
// For electric field this should be true
// For magnetic field this should be false
private:
G4FieldManager(const G4FieldManager&);
G4FieldManager& operator=(const G4FieldManager&);
// Private copy constructor and assignment operator.
private:
G4Field* fDetectorField;
G4ChordFinder* fChordFinder;
G4bool fAllocatedChordFinder; // Did we used "new" to
// create fChordFinder ?
G4bool fFieldChangesEnergy;
// Values for the required accuracies
//
G4double fDelta_One_Step_Value; // for one tracking/physics step
G4double fDelta_Intersection_Val; // for boundary intersection
G4double fDefault_Delta_One_Step_Value; // = 0.25 * mm;
G4double fDefault_Delta_Intersection_Val; // = 0.1 * mm;
// Values for the small possible relative accuracy of a step
// (corresponding to the greatest possible integration accuracy)
G4double fEpsilonMinDefault; // Can be 1.0e-5 to 1.0e-10 ...
G4double fEpsilonMaxDefault; // Can be 1.0e-3 to 1.0e-8 ...
G4double fEpsilonMin;
G4double fEpsilonMax;
//TAO (nb if we add a value which is initialized in constructor: respect the order of declaration - initialization. therefor fMagComponent should be declared after fEpsilonMax .-cf constructors)
G4bool fMagComponent;
};
// Our current design envisions that one Field manager is valid for a region of the detector.
// (eg a detector with electric E and magnetic B field will now treat
// them as one field - and could treat any other field of importance
// as additional components of a single field.)
// Does it make sense to have several instead ?
// Is the lack of elegance of the design (of G4Field) made up
// for by the simplification it allows ?
// Implementation of inline functions
#include "G4FieldManager.icc"
#endif /* G4FIELDMANAGER_HH */

133
geant4/LEMuSR/G4Modified/G4FieldManager.icc Normal file
View File

@ -0,0 +1,133 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4FieldManager.icc,v 1.9 2003/11/08 03:39:39 japost Exp $
// GEANT4 tag $Name: geant4-06-00-patch-01 $
//
//
// G4FieldManager inline implementation
//
// -------------------------------------------------------------------
inline
const G4Field* G4FieldManager::GetDetectorField() const
{
// If pointer is null, should this raise an exception ??
return fDetectorField;
}
inline
G4bool G4FieldManager::DoesFieldExist() const
{
return (fDetectorField != 0);
}
inline
void G4FieldManager::SetChordFinder(G4ChordFinder *aChordFinder)
{
fChordFinder= aChordFinder;
}
inline
G4ChordFinder* G4FieldManager::GetChordFinder()
{
return fChordFinder;
}
inline
G4double G4FieldManager::GetDeltaIntersection() const
{
return fDelta_Intersection_Val;
}
inline
G4double G4FieldManager::GetDeltaOneStep() const
{
return fDelta_One_Step_Value;
}
inline
void G4FieldManager::SetDeltaOneStep(G4double valDeltaOneStep)
{
fDelta_One_Step_Value= valDeltaOneStep;
}
inline
void G4FieldManager::SetDeltaIntersection(G4double valDeltaIntersection)
{
fDelta_Intersection_Val = valDeltaIntersection;
}
inline
void G4FieldManager::SetAccuraciesWithDeltaOneStep(G4double valDeltaOneStep)
{
fDelta_One_Step_Value= valDeltaOneStep;
fDelta_Intersection_Val = 0.4 * fDelta_One_Step_Value;
}
inline G4bool G4FieldManager::DoesFieldChangeEnergy() const
{ return fFieldChangesEnergy;}
inline void G4FieldManager::SetFieldChangesEnergy(G4bool value)
{ fFieldChangesEnergy = value; }
inline G4bool G4FieldManager:: FieldHasMagComponent() const
{ return fMagComponent; }
inline void G4FieldManager::SetFieldMagComponent(G4bool value)
{ fMagComponent = value;}
// Minimum for Relative accuracy of any Step
inline
G4double G4FieldManager::GetMinimumEpsilonStep() const
{
return fEpsilonMin;
}
inline
void G4FieldManager::SetMinimumEpsilonStep( G4double newEpsMin )
{
if( (newEpsMin > 0.0) && (fabs(1.0+newEpsMin) > 1.0) )
{
fEpsilonMin = newEpsMin;
}
}
// Maximum for Relative accuracy of any Step
inline
G4double G4FieldManager::GetMaximumEpsilonStep() const
{
return fEpsilonMax;
}
inline
void G4FieldManager::SetMaximumEpsilonStep( G4double newEpsMax )
{
if( (newEpsMax > 0.0)
&& (newEpsMax >= fEpsilonMin )
&& (fabs(1.0+newEpsMax)>1.0) )
{
fEpsilonMax = newEpsMax;
}
}

951
geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc Normal file
View File

@ -0,0 +1,951 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4MultipleScattering52.cc,v 1.2 2004/12/01 19:37:14 vnivanch Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
//
// -----------------------------------------------------------------------------
// 16/05/01 value of cparm changed , L.Urban
// 18/05/01 V.Ivanchenko Clean up against Linux ANSI compilation
// 07/08/01 new methods Store/Retrieve PhysicsTable (mma)
// 23-08-01 new angle and z distribution,energy dependence reduced,
// Store,Retrieve methods commented out temporarily, L.Urban
// 27-08-01 in BuildPhysicsTable:aParticleType.GetParticleName()=="mu+" (mma)
// 28-08-01 GetContinuousStepLimit and AlongStepDoIt moved from .icc file (mma)
// 03-09-01 value of data member factlim changed, L.Urban
// 10-09-01 small change in GetContinuousStepLimit, L.Urban
// 11-09-01 G4MultipleScatteringx put as default G4MultipleScattering
// store/retrieve physics table reactivated (mma)
// 13-09-01 corr. in ComputeTransportCrossSection, L.Urban
// 14-09-01 protection in GetContinuousStepLimit, L.Urban
// 17-09-01 migration of Materials to pure STL (mma)
// 27-09-01 value of data member factlim changed, L.Urban
// 31-10-01 big fixed in PostStepDoIt,L.Urban
// 24-04-02 some minor changes in boundary algorithm, L.Urban
// 06-05-02 bug fixed in GetContinuousStepLimit, L.Urban
// 24-05-02 changes in angle distribution and boundary algorithm, L.Urban
// 11-06-02 bug fixed in ComputeTransportCrossSection, L.Urban
// 12-08-02 bug fixed in PostStepDoIt (lateral displacement), L.Urban
// 15-08-02 new angle distribution, L.Urban
// 26-09-02 angle distribution + boundary algorithm modified, L.Urban
// 15-10-02 temporary fix for proton scattering
// 30-10-02 modified angle distribution,mods in boundary algorithm,
// changes in data members, L.Urban
// 30-10-02 rename variable cm - Ecm, V.Ivanchenko
// 11-12-02 precision problem in ComputeTransportCrossSection
// for small Tkin/for heavy particles cured, L.Urban
// 05-02-03 changes in data members, new sampling for geom.
// path length, step dependence reduced with new
// method
// 17-03-03 cut per region, V.Ivanchenko
// 13-04-03 add initialisation in GetContinuesStepLimit
// + change table size (V.Ivanchenko)
// 26-04-03 fix problems of retrieve tables (M.Asai)
// 23-05-03 important change in angle distribution for muons/hadrons
// the central part now is similar to the Highland parametrization +
// minor correction in angle sampling algorithm (for all particles)
// (L.Urban)
// 24-05-03 bug in nuclear size corr.computation fixed thanks to Vladimir(L.Urban)
// 30-05-03 misprint in PostStepDoIt corrected(L.Urban)
// 08-08-03 This class is frozen at the release 5.2 (V.Ivanchenko)
// 08-11-04 Remove Store/Retrieve tables (V.Ivantchenko)
// -----------------------------------------------------------------------------
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4MultipleScattering52.hh"
#include "G4StepStatus.hh"
#include "G4Navigator.hh"
#include "G4TransportationManager.hh"
#include "Randomize.hh"
#include "G4ProductionCutsTable.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
using namespace std;
G4MultipleScattering52::G4MultipleScattering52(const G4String& processName)
: G4VContinuousDiscreteProcess(processName),
theTransportMeanFreePathTable(0),
taubig(8.0),tausmall(1.e-14),taulim(1.e-5),
LowestKineticEnergy(0.1*keV),
HighestKineticEnergy(100.*TeV),
TotBin(100),
materialIndex(0),
tLast (0.0),
zLast (0.0),
boundary(true),
facrange(0.199),tlimit(1.e10*mm),tlimitmin(1.e-7*mm),
cf(1.001),
stepno(0),stepnolastmsc(-1000000),nsmallstep(5),
laststep(0.),
valueGPILSelectionMSC(NotCandidateForSelection),
zmean(0.),samplez(true),
range(1.),T0(1.),T1(1.),lambda0(1.),lambda1(-1.),
Tlow(0.),alam(1.),blam(1.),dtrl(0.15),
lambdam(-1.),clam(1.),zm(1.),cthm(1.),
fLatDisplFlag(true),
NuclCorrPar (0.0615),
FactPar(0.40),
facxsi(1.)
{ }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4MultipleScattering52::~G4MultipleScattering52()
{
if(theTransportMeanFreePathTable)
{
theTransportMeanFreePathTable->clearAndDestroy();
delete theTransportMeanFreePathTable;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void G4MultipleScattering52::BuildPhysicsTable(
const G4ParticleDefinition& aParticleType)
{
// set values of some data members
if((aParticleType.GetParticleName() == "e-") ||
(aParticleType.GetParticleName() == "e+"))
{
// parameters for e+/e-
alfa1 = 1.45 ;
alfa2 = 0.60 ;
alfa3 = 0.30 ;
b = 1. ;
xsi = facxsi*2.22 ;
c0 = 2.30 ;
}
else
{
// parameters for heavy particles
alfa1 = 1.10 ;
alfa2 = 0.14 ;
alfa3 = 0.07 ;
b = 1. ;
xsi = facxsi*2.70 ;
c0 = 1.40 ;
}
// ..............................
Tlow = aParticleType.GetPDGMass();
// tables are built for MATERIALS
const G4double sigmafactor = twopi*classic_electr_radius*
classic_electr_radius;
G4double KineticEnergy,AtomicNumber,AtomicWeight,sigma,lambda;
G4double density;
// destroy old tables if any
if (theTransportMeanFreePathTable)
{
theTransportMeanFreePathTable->clearAndDestroy();
delete theTransportMeanFreePathTable;
}
// create table
const G4ProductionCutsTable* theCoupleTable=
G4ProductionCutsTable::GetProductionCutsTable();
size_t numOfCouples = theCoupleTable->GetTableSize();
theTransportMeanFreePathTable = new G4PhysicsTable(numOfCouples);
// loop for materials
for (size_t i=0; i<numOfCouples; i++)
{
// create physics vector and fill it
G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
LowestKineticEnergy,HighestKineticEnergy,TotBin);
// get elements in the material
const G4MaterialCutsCouple* couple = theCoupleTable->
GetMaterialCutsCouple(i);
const G4Material* material = couple->GetMaterial();
const G4ElementVector* theElementVector = material->GetElementVector();
const G4double* NbOfAtomsPerVolume =
material->GetVecNbOfAtomsPerVolume();
const G4int NumberOfElements = material->GetNumberOfElements();
density = material->GetDensity();
// loop for kinetic energy values
for (G4int i=0; i<TotBin; i++)
{
KineticEnergy = aVector->GetLowEdgeEnergy(i);
sigma = 0.;
// loop for element in the material
for (G4int iel=0; iel<NumberOfElements; iel++)
{
AtomicNumber = (*theElementVector)[iel]->GetZ();
AtomicWeight = (*theElementVector)[iel]->GetA();
sigma += NbOfAtomsPerVolume[iel]*
ComputeTransportCrossSection(aParticleType,KineticEnergy,
AtomicNumber,AtomicWeight);
}
sigma *= sigmafactor;
lambda = 1./sigma;
aVector->PutValue(i,lambda);
}
theTransportMeanFreePathTable->insert(aVector);
}
if((aParticleType.GetParticleName() == "e-" ) ||
(aParticleType.GetParticleName() == "mu+" ) ||
(aParticleType.GetParticleName() == "Mu" ) ||
(aParticleType.GetParticleName() == "proton") ) PrintInfoDefinition();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double G4MultipleScattering52::ComputeTransportCrossSection(
const G4ParticleDefinition& aParticleType,
G4double KineticEnergy,
G4double AtomicNumber,G4double AtomicWeight)
{
const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
Bohr_radius*Bohr_radius/(hbarc*hbarc);
const G4double epsmin = 1.e-4 , epsmax = 1.e10;
const G4double Zdat[15] = { 4., 6.,13.,20.,26.,29.,32.,38.,47.,
50.,56.,64.,74.,79.,82. };
const G4double Tdat[23] = {0.0001*MeV,0.0002*MeV,0.0004*MeV,0.0007*MeV,
0.001*MeV,0.002*MeV,0.004*MeV,0.007*MeV,
0.01*MeV,0.02*MeV,0.04*MeV,0.07*MeV,
0.1*MeV,0.2*MeV,0.4*MeV,0.7*MeV,
1.*MeV,2.*MeV,4.*MeV,7.*MeV,10.*MeV,20.*MeV,
10000.0*MeV};
// corr. factors for e-/e+ lambda
G4double celectron[15][23] =
{{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
1.112,1.108,1.100,1.093,1.089,1.087,0.7235 },
{1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
1.109,1.105,1.097,1.090,1.086,1.082,0.7925 },
{2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
1.131,1.124,1.113,1.104,1.099,1.098,0.9147 },
{3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
1.112,1.105,1.096,1.089,1.085,1.098,0.9700 },
{6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
1.073,1.070,1.064,1.059,1.056,1.056,1.0022 },
{9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
1.074,1.070,1.063,1.059,1.056,1.052,1.0158 },
{11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
1.068,1.064,1.059,1.054,1.051,1.050,1.0284 },
{18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
1.039,1.037,1.034,1.031,1.030,1.036,1.0515 },
{18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
1.031,1.028,1.024,1.022,1.021,1.024,1.0834 },
{14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
1.020,1.017,1.015,1.013,1.013,1.020,1.0937 },
{14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
0.995,0.993,0.993,0.993,0.993,1.011,1.1140 },
{22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
0.974,0.972,0.973,0.974,0.975,0.987,1.1410 },
{50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
0.950,0.947,0.949,0.952,0.954,0.963,1.1750 },
{65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
0.941,0.938,0.940,0.944,0.946,0.954,1.1922 },
// {45.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239, // paper.....
{55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
0.933,0.930,0.933,0.936,0.939,0.949,1.2026 }};
G4double cpositron[15][23] = {
{2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
1.131,1.126,1.117,1.108,1.103,1.100,0.7235 },
{3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
1.138,1.132,1.122,1.113,1.108,1.102,0.7925 },
{7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
1.203,1.190,1.173,1.159,1.151,1.145,0.9147 },
{9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
1.225,1.210,1.191,1.175,1.166,1.174,0.9700 },
{17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
1.217,1.203,1.184,1.169,1.160,1.151,1.0022 },
{24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
1.237,1.222,1.201,1.184,1.174,1.159,1.0158 },
{23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
1.252,1.234,1.212,1.194,1.183,1.170,1.0284 },
{22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
1.254,1.237,1.214,1.195,1.185,1.179,1.0515 },
{33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
1.312,1.288,1.258,1.235,1.221,1.205,1.0834 },
{32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
1.320,1.294,1.264,1.240,1.226,1.214,1.0937 },
{29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
1.328,1.302,1.270,1.245,1.231,1.233,1.1140 },
{38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
1.361,1.330,1.294,1.267,1.251,1.239,1.1410 },
{49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
1.409,1.372,1.330,1.298,1.280,1.258,1.1750 },
{59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
1.442,1.400,1.354,1.319,1.299,1.272,1.1922 },
{56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
1.456,1.412,1.364,1.328,1.307,1.282,1.2026 }};
G4double sigma;
G4double Z23 = 2.*log(AtomicNumber)/3.; Z23 = exp(Z23);
G4double ParticleMass = aParticleType.GetPDGMass();
G4double ParticleKineticEnergy = KineticEnergy ;
// correction if particle .ne. e-/e+
// compute equivalent kinetic energy
// lambda depends on p*beta ....
G4double Mass = ParticleMass ;
if((aParticleType.GetParticleName() != "e-") &&
(aParticleType.GetParticleName() != "e+") )
{
G4double TAU = KineticEnergy/Mass ;
G4double c = Mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
G4double w = c-2. ;
G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
KineticEnergy = electron_mass_c2*tau ;
Mass = electron_mass_c2 ;
}
G4double Charge = aParticleType.GetPDGCharge();
G4double ChargeSquare = Charge*Charge/(eplus*eplus);
G4double TotalEnergy = KineticEnergy + Mass ;
G4double beta2 = KineticEnergy*(TotalEnergy+Mass)
/(TotalEnergy*TotalEnergy);
G4double bg2 = KineticEnergy*(TotalEnergy+Mass)
/(Mass*Mass);
G4double eps = epsfactor*bg2/Z23;
if (eps<epsmin) sigma = 2.*eps*eps;
else if(eps<epsmax) sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
else sigma = log(2.*eps)-1.+1./eps;
sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
// nuclear size effect correction for high energy
// ( a simple approximation at present)
G4double corrnuclsize,a,x0,w1,w2,w;
x0 = 1. - NuclCorrPar*ParticleMass/(ParticleKineticEnergy*
exp(log(AtomicWeight/(g/mole))/3.));
if ( (x0 < -1.) || (ParticleKineticEnergy <= 10.*MeV))
{ x0 = -1.; corrnuclsize = 1.;}
else
{ a = 1.+1./eps;
if (eps > epsmax) w1=log(2.*eps)+1./eps-3./(8.*eps*eps);
else w1=log((a+1.)/(a-1.))-2./(a+1.);
w = 1./((1.-x0)*eps);
if (w < epsmin) w2=-log(w)-1.+2.*w-1.5*w*w;
else w2 = log((a-x0)/(a-1.))-(1.-x0)/(a-x0);
corrnuclsize = w1/w2;
corrnuclsize = exp(-FactPar*ParticleMass/ParticleKineticEnergy)*
(corrnuclsize-1.)+1.;
}
// interpolate in AtomicNumber and beta2
// get bin number in Z
G4int iZ = 14;
while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
if (iZ==14) iZ = 13;
if (iZ==-1) iZ = 0 ;
G4double Z1 = Zdat[iZ];
G4double Z2 = Zdat[iZ+1];
G4double ratZ = (AtomicNumber-Z1)/(Z2-Z1);
// get bin number in T (beta2)
G4int iT = 22;
while ((iT>=0)&&(Tdat[iT]>=KineticEnergy)) iT -= 1;
if(iT==22) iT = 21;
if(iT==-1) iT = 0 ;
// calculate betasquare values
G4double T = Tdat[iT], E = T + electron_mass_c2;
G4double b2small = T*(E+electron_mass_c2)/(E*E);
T = Tdat[iT+1]; E = T + electron_mass_c2;
G4double b2big = T*(E+electron_mass_c2)/(E*E);
G4double ratb2 = (beta2-b2small)/(b2big-b2small);
G4double c1,c2,cc1,cc2,corr;
if (Charge < 0.)
{
c1 = celectron[iZ][iT];
c2 = celectron[iZ+1][iT];
cc1 = c1+ratZ*(c2-c1);
c1 = celectron[iZ][iT+1];
c2 = celectron[iZ+1][iT+1];
cc2 = c1+ratZ*(c2-c1);
corr = cc1+ratb2*(cc2-cc1);
sigma /= corr;
}
if (Charge > 0.)
{
c1 = cpositron[iZ][iT];
c2 = cpositron[iZ+1][iT];
cc1 = c1+ratZ*(c2-c1);
c1 = cpositron[iZ][iT+1];
c2 = cpositron[iZ+1][iT+1];
cc2 = c1+ratZ*(c2-c1);
corr = cc1+ratb2*(cc2-cc1);
sigma /= corr;
}
// nucl. size correction for particles other than e+/e- only at present !!!!
if((aParticleType.GetParticleName() != "e-") &&
(aParticleType.GetParticleName() != "e+") )
sigma /= corrnuclsize;
return sigma;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double G4MultipleScattering52::GetContinuousStepLimit(
const G4Track& track,
G4double,
G4double currentMinimumStep,
G4double&)
{
G4double zPathLength,tPathLength;
const G4DynamicParticle* aParticle;
G4double tau,zt,cz,cz1,grej,grej0;
const G4double expmax = 100., ztmax = (2.*expmax+1.)/(2.*expmax+3.) ;
const G4double tmax = 1.e20*mm ;
G4bool isOut;
// this process is not a candidate for selection by default
valueGPILSelectionMSC = NotCandidateForSelection;
tPathLength = currentMinimumStep;
const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
materialIndex = couple->GetIndex();
aParticle = track.GetDynamicParticle();
T0 = aParticle->GetKineticEnergy();
lambda0 = (*theTransportMeanFreePathTable)
(materialIndex)->GetValue(T0,isOut);
range = G4EnergyLossTables::GetRange(aParticle->GetDefinition(),
T0,couple);
//VI Initialisation at the beginning of the step
cthm = 1.;
lambda1 = -1.;
lambdam = -1.;
alam = range;
blam = 1.+alam/lambda0 ;
zm = 1.;
// special treatment near boundaries ?
if (boundary && range >= currentMinimumStep)
{
// step limitation at boundary ?
stepno = track.GetCurrentStepNumber() ;
if(stepno == 1)
{
stepnolastmsc = -1000000 ;
tlimit = 1.e10 ;
}
if(stepno > 1)
{
if(track.GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary)
{
stepnolastmsc = stepno ;
// if : diff.treatment for small/not small Z
if(range > lambda0)
tlimit = facrange*range ;
else
tlimit = facrange*lambda0 ;
if(tlimit < tlimitmin) tlimit = tlimitmin ;
laststep = tlimit ;
if(tPathLength > tlimit)
{
tPathLength = tlimit ;
valueGPILSelectionMSC = CandidateForSelection;
}
}
else if(stepno > stepnolastmsc)
{
if((stepno - stepnolastmsc) < nsmallstep)
{
if(tPathLength > tlimit)
{
laststep *= cf ;
tPathLength = laststep ;
valueGPILSelectionMSC = CandidateForSelection;
}
}
}
}
}
// do the true -> geom transformation
zmean = tPathLength;
tau = tPathLength/lambda0 ;
if (tau < tausmall || range < currentMinimumStep) zPathLength = tPathLength;
else
{
if(tPathLength/range < dtrl) zmean = lambda0*(1.-exp(-tau));
else
{
T1 = G4EnergyLossTables::GetPreciseEnergyFromRange(
aParticle->GetDefinition(),range-tPathLength,couple);
lambda1 = (*theTransportMeanFreePathTable)
(materialIndex)->GetValue(T1,isOut);
if(T0 < Tlow)
alam = range ;
else
alam = lambda0*tPathLength/(lambda0-lambda1) ;
blam = 1.+alam/lambda0 ;
if(tPathLength/range < 2.*dtrl)
{
zmean = alam*(1.-exp(blam*log(1.-tPathLength/alam)))/blam ;
lambdam = -1. ;
}
else
{
G4double w = 1.-0.5*tPathLength/alam ;
lambdam = lambda0*w ;
clam = 1.+alam/lambdam ;
cthm = exp(alam*log(w)/lambda0) ;
zm = alam*(1.-exp(blam*log(w)))/blam ;
zmean = zm + alam*(1.-exp(clam*log(w)))*cthm/clam ;
}
}
// sample z
zt = zmean/tPathLength ;
if (samplez && (zt < ztmax) && (zt > 0.5))
{
cz = 0.5*(3.*zt-1.)/(1.-zt) ;
if(tPathLength < exp(log(tmax)/(2.*cz)))
{
cz1 = 1.+cz ;
grej0 = exp(cz1*log(cz*tPathLength/cz1))/cz ;
do
{
zPathLength = tPathLength*exp(log(G4UniformRand())/cz1) ;
grej = exp(cz*log(zPathLength))*(tPathLength-zPathLength)/grej0 ;
} while (grej < G4UniformRand()) ;
}
else zPathLength = zmean;
}
else zPathLength = zmean;
}
// protection against z > lambda
if(zPathLength > lambda0)
zPathLength = lambda0 ;
tLast = tPathLength;
zLast = zPathLength;
return zPathLength;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4VParticleChange* G4MultipleScattering52::AlongStepDoIt(
const G4Track& track,const G4Step& step)
{
// only a geom path->true path transformation is performed
fParticleChange.Initialize(track);
G4double geomPathLength = step.GetStepLength();
G4double truePathLength = 0. ;
//VI change order of if operators
if(geomPathLength == zLast) truePathLength = tLast;
else if(geomPathLength/lambda0 < tausmall) truePathLength = geomPathLength;
else
{
if(lambda1 < 0.) truePathLength = -lambda0*log(1.-geomPathLength/lambda0) ;
else if(lambdam < 0.)
{
if(blam*geomPathLength/alam < 1.)
truePathLength = alam*(1.-exp(log(1.-blam*geomPathLength/alam)/
blam)) ;
else
truePathLength = tLast;
}
else
{
if(geomPathLength <= zm)
{
if(blam*geomPathLength/alam < 1.)
truePathLength = alam*(1.-exp(log(1.-blam*geomPathLength/alam)/
blam)) ;
else
truePathLength = 0.5*tLast;
lambdam = -1. ;
}
else
{
if(clam*(geomPathLength-zm)/(alam*cthm) < 1.)
truePathLength = 0.5*tLast + alam*(1.-
exp(log(1.-clam*(geomPathLength-zm)/(alam*cthm)))/clam) ;
else
truePathLength = tLast ;
}
}
// protection ....
if(truePathLength > tLast)
truePathLength = tLast ;
}
//VI truePath length cannot be smaller than geomPathLength
if (truePathLength < geomPathLength) truePathLength = geomPathLength;
fParticleChange.ProposeTrueStepLength(truePathLength);
return &fParticleChange;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4VParticleChange* G4MultipleScattering52::PostStepDoIt(
const G4Track& trackData,
const G4Step& stepData)
{
// angle distribution parameters
const G4double kappa = 2.5, kappapl1 = kappa+1., kappami1 = kappa-1. ;
fParticleChange.Initialize(trackData);
G4double truestep = stepData.GetStepLength();
const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
G4double KineticEnergy = aParticle->GetKineticEnergy();
G4double Mass = aParticle->GetDefinition()->GetPDGMass() ;
// do nothing for stopped particles !
if(KineticEnergy > 0.)
{
// change direction first ( scattering )
G4double cth = 1.0 ;
G4double tau = truestep/lambda0 ;
if (tau < tausmall) cth = 1.;
else if(tau > taubig) cth = -1.+2.*G4UniformRand();
else
{
if(lambda1 > 0.)
{
if(lambdam < 0.)
tau = -alam*log(1.-truestep/alam)/lambda0 ;
else
tau = -log(cthm)-alam*log(1.-(truestep-0.5*tLast)/alam)/lambdam ;
}
if(tau > taubig) cth = -1.+2.*G4UniformRand();
else
{
const G4double amax=25. ;
const G4double tau0 = 0.02 ;
const G4double c_highland = 13.6*MeV, corr_highland=0.038 ;
const G4double x1fac1 = exp(-xsi) ;
const G4double x1fac2 = (1.-(1.+xsi)*x1fac1)/(1.-x1fac1) ;
const G4double x1fac3 = 1.3 ; // x1fac3 >= 1. !!!!!!!!!
G4double a,x0,c,xmean1,xmean2,
xmeanth,prob,qprob ;
G4double ea,eaa,b1,bx,eb1,ebx,cnorm1,cnorm2,f1x0,f2x0,w ;
// for heavy particles take the width of the cetral part
// from the Highland formula
// (Particle Physics Booklet, July 2002, eq. 26.10)
if(Mass > electron_mass_c2) // + other conditions (beta, x/X0,...?)
{
G4double Q = fabs(aParticle->GetDefinition()->GetPDGCharge()) ;
G4double X0 = trackData.GetMaterialCutsCouple()->
GetMaterial()->GetRadlen() ;
G4double xx0 = truestep/X0 ;
G4double betacp = KineticEnergy*(KineticEnergy+2.*Mass)/
(KineticEnergy+Mass) ;
G4double theta0=c_highland*Q*sqrt(xx0)*
(1.+corr_highland*log(xx0))/betacp ;
if(theta0 > tausmall)
a = 0.5/(1.-cos(theta0)) ;
else
a = 1./(theta0*theta0) ;
}
else
{
w = log(tau/tau0) ;
if(tau < tau0)
a = (alfa1-alfa2*w)/tau ;
else
a = (alfa1+alfa3*w)/tau ;
}
xmeanth = exp(-tau) ;
x0 = 1.-xsi/a ;
if(x0 < -1.) x0 = -1. ;
if(x0 == -1.)
{
// 1 model fuction only
// in order to have xmean1 > xmeanth -> qprob < 1
if((1.-1./a) < xmeanth)
a = 1./(1.-xmeanth) ;
if(a*(1.-x0) < amax)
ea = exp(-a*(1.-x0)) ;
else
ea = 0. ;
eaa = 1.-ea ;
xmean1 = 1.-1./a+(1.-x0)*ea/eaa ;
c = 2. ;
b1 = b+1. ;
bx = b1 ;
eb1 = b1 ;
ebx = b1 ;
xmean2 = 0. ;
prob = 1. ;
qprob = xmeanth/xmean1 ;
}
else
{
// 2 model fuctions
// in order to have xmean1 > xmeanth
if((1.-x1fac2/a) < xmeanth)
{
a = x1fac3*x1fac2/(1.-xmeanth) ;
if(a*(1.-x0) < amax)
ea = exp(-a*(1.-x0)) ;
else
ea = 0. ;
eaa = 1.-ea ;
xmean1 = 1.-1./a+(1.-x0)*ea/eaa ;
}
else
{
ea = x1fac1 ;
eaa = 1.-x1fac1 ;
xmean1 = 1.-x1fac2/a ;
}
// from continuity of the 1st derivatives
c = a*(b-x0) ;
if(a*tau < c0)
c = c0*(b-x0)/tau ;
if(c == 1.) c=1.000001 ;
if(c == 2.) c=2.000001 ;
if(c == 3.) c=3.000001 ;
b1 = b+1. ;
bx=b-x0 ;
eb1=exp((c-1.)*log(b1)) ;
ebx=exp((c-1.)*log(bx)) ;
xmean2 = (x0*eb1+ebx+(eb1*bx-b1*ebx)/(2.-c))/(eb1-ebx) ;
cnorm1 = a/eaa ;
f1x0 = cnorm1*exp(-a*(1.-x0)) ;
cnorm2 = (c-1.)*eb1*ebx/(eb1-ebx) ;
f2x0 = cnorm2/exp(c*log(b-x0)) ;
// from continuity at x=x0
prob = f2x0/(f1x0+f2x0) ;
// from xmean = xmeanth
qprob = (f1x0+f2x0)*xmeanth/(f2x0*xmean1+f1x0*xmean2) ;
}
// protection against prob or qprob > 1 and
// prob or qprob < 0
// ***************************************************************
if((qprob > 1.) || (qprob < 0.) || (prob > 1.) || (prob < 0.))
{
// this print possibility has been left intentionally
// for debugging purposes ..........................
G4bool pr = false ;
// pr = true ;
if(pr)
{
const G4double prlim = 0.10 ;
if((fabs((xmeanth-xmean2)/(xmean1-xmean2)-prob)/prob > prlim) ||
((xmeanth-xmean2)/(xmean1-xmean2) > 1.) ||
((xmeanth-xmean2)/(xmean1-xmean2) < 0.) )
{
G4cout.precision(5) ;
G4cout << "\nparticle=" << aParticle->GetDefinition()->
GetParticleName() << " in material "
<< trackData.GetMaterialCutsCouple()->
GetMaterial()->GetName() << " with kinetic energy "
<< KineticEnergy << " MeV," << G4endl ;
G4cout << " step length="
<< truestep << " mm" << G4endl ;
G4cout << "p=" << prob << " q=" << qprob << " -----> "
<< "p=" << (xmeanth-xmean2)/(xmean1-xmean2)
<< " q=" << 1. << G4endl ;
}
}
qprob = 1. ;
prob = (xmeanth-xmean2)/(xmean1-xmean2) ;
}
// **************************************************************
// sampling of costheta
if(G4UniformRand() < qprob)
{
if(G4UniformRand() < prob)
cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
else
cth = b-b1*bx/exp(log(ebx-G4UniformRand()*(ebx-eb1))/(c-1.)) ;
}
else
cth = -1.+2.*G4UniformRand() ;
}
}
G4double sth = sqrt(1.-cth*cth);
G4double phi = twopi*G4UniformRand();
G4double dirx = sth*cos(phi), diry = sth*sin(phi), dirz = cth;
G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection();
G4ThreeVector newDirection(dirx,diry,dirz);
newDirection.rotateUz(ParticleDirection);
fParticleChange.ProposeMomentumDirection(newDirection.x(),
newDirection.y(),
newDirection.z());
if (fLatDisplFlag)
{
// compute mean lateral displacement, only for safety > tolerance !
G4double safetyminustolerance = stepData.GetPostStepPoint()->GetSafety();
G4double rmean, etau;
if (safetyminustolerance > 0.)
{
if (tau < tausmall) rmean = 0.;
else if(tau < taulim) rmean = kappa*tau*tau*tau*(1.-kappapl1*tau/4.)/6.;
else
{
if(tau<taubig) etau = exp(-tau);
else etau = 0.;
rmean = -kappa*tau;
rmean = -exp(rmean)/(kappa*kappami1);
rmean += tau-kappapl1/kappa+kappa*etau/kappami1;
}
if (rmean>0.) rmean = 2.*lambda0*sqrt(rmean/3.);
else rmean = 0.;
// for rmean > 0) only
if (rmean > 0.)
{
if (rmean>safetyminustolerance) rmean = safetyminustolerance;
// sample direction of lateral displacement
phi = twopi*G4UniformRand();
dirx = cos(phi); diry = sin(phi); dirz = 0.;
G4ThreeVector latDirection(dirx,diry,dirz);
latDirection.rotateUz(ParticleDirection);
// compute new endpoint of the Step
G4ThreeVector newPosition = stepData.GetPostStepPoint()->GetPosition()
+ rmean*latDirection;
G4Navigator* navigator =
G4TransportationManager::GetTransportationManager()
->GetNavigatorForTracking();
navigator->LocateGlobalPointWithinVolume(newPosition);
fParticleChange.ProposePosition(newPosition);
}
}
}
}
return &fParticleChange;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void G4MultipleScattering52::PrintInfoDefinition()
{
G4String comments = " Tables of transport mean free paths.";
comments += "\n New model of MSC , computes the lateral \n";
comments += " displacement of the particle , too.";
G4cout << G4endl << GetProcessName() << ": " << comments
<< "\n PhysicsTables from "
<< G4BestUnit(LowestKineticEnergy ,"Energy")
<< " to " << G4BestUnit(HighestKineticEnergy,"Energy")
<< " in " << TotBin << " bins. \n";
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

190
geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc Normal file
View File

@ -0,0 +1,190 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4MuonDecayChannel.cc,v 1.11 2004/12/10 18:02:04 gcosmo Exp $
// GEANT4 tag $Name: geant4-07-00-cand-05 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 30 May 1997 H.Kurashige
//
// Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige
// ------------------------------------------------------------
#include "G4ParticleDefinition.hh"
#include "G4DecayProducts.hh"
#include "G4VDecayChannel.hh"
#include "G4MuonDecayChannel.hh"
#include "Randomize.hh"
#include "G4LorentzVector.hh"
#include "G4LorentzRotation.hh"
G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName,
G4double theBR)
:G4VDecayChannel("Muon Decay",1)
{
// set names for daughter particles
if (theParentName == "mu+") {
SetBR(theBR);
SetParent("mu+");
SetNumberOfDaughters(3);
SetDaughter(0, "e+");
SetDaughter(1, "nu_e");
SetDaughter(2, "anti_nu_mu");
} else if (theParentName == "Mu") {
SetBR(theBR);
SetParent("Mu");
SetNumberOfDaughters(3);
SetDaughter(0, "e+");
SetDaughter(1, "nu_e");
SetDaughter(2, "anti_nu_mu");
} else if (theParentName == "mu-") {
SetBR(theBR);
SetParent("mu-");
SetNumberOfDaughters(3);
SetDaughter(0, "e-");
SetDaughter(1, "anti_nu_e");
SetDaughter(2, "nu_mu");
} else {
#ifdef G4VERBOSE
if (GetVerboseLevel()>0) {
G4cout << "G4MuonDecayChannel:: constructor :";
G4cout << " parent particle is not muon but ";
G4cout << theParentName << G4endl;
}
#endif
}
}
G4MuonDecayChannel::~G4MuonDecayChannel()
{
}
G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double)
{
// this version neglects muon polarization
// assumes the pure V-A coupling
// gives incorrect energy spectrum for neutrinos
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt ";
#endif
if (parent == 0) FillParent();
if (daughters == 0) FillDaughters();
// parent mass
G4double parentmass = parent->GetPDGMass();
//daughters'mass
G4double daughtermass[3];
G4double sumofdaughtermass = 0.0;
for (G4int index=0; index<3; index++){
daughtermass[index] = daughters[index]->GetPDGMass();
sumofdaughtermass += daughtermass[index];
}
//create parent G4DynamicParticle at rest
G4ThreeVector dummy;
G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
//create G4Decayproducts
G4DecayProducts *products = new G4DecayProducts(*parentparticle);
delete parentparticle;
// calculate daughter momentum
G4double daughtermomentum[3];
G4double energy;
// calcurate electron energy
G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
G4double x;
G4double r;
do {
do {
r = G4UniformRand();
x = xmax*G4UniformRand();
} while (r > (3.0 - 2.0*x)*x*x);
energy = x*parentmass/2.0 - daughtermass[0];
} while (energy <0.0);
//create daughter G4DynamicParticle
// daughter 0 (electron)
daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
G4double costheta, sintheta, phi, sinphi, cosphi;
costheta = 2.*G4UniformRand()-1.0;
sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
phi = twopi*G4UniformRand()*rad;
sinphi = std::sin(phi);
cosphi = std::cos(phi);
G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
G4DynamicParticle * daughterparticle
= new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
products->PushProducts(daughterparticle);
// daughter 1 ,2 (nutrinos)
// create neutrinos in the C.M frame of two neutrinos
G4double energy2 = parentmass*(1.0 - x/2.0);
G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
G4double beta = -1.0*daughtermomentum[0]/energy2;
G4double costhetan = 2.*G4UniformRand()-1.0;
G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
G4double phin = twopi*G4UniformRand()*rad;
G4double sinphin = std::sin(phin);
G4double cosphin = std::cos(phin);
G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
G4DynamicParticle * daughterparticle1
= new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
G4DynamicParticle * daughterparticle2
= new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
// boost to the muon rest frame
G4LorentzVector p4;
p4 = daughterparticle1->Get4Momentum();
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
daughterparticle1->Set4Momentum(p4);
p4 = daughterparticle2->Get4Momentum();
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
daughterparticle2->Set4Momentum(p4);
products->PushProducts(daughterparticle1);
products->PushProducts(daughterparticle2);
daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
// output message
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) {
G4cout << "G4MuonDecayChannel::DecayIt ";
G4cout << " create decay products in rest frame " <<G4endl;
products->DumpInfo();
}
#endif
return products;
}

63
geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh Normal file
View File

@ -0,0 +1,63 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4MuonDecayChannel.hh,v 1.5 2001/07/11 10:01:56 gunter Exp $
// GEANT4 tag $Name: geant4-07-00-cand-01 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 30 May 1997 H.Kurashige
// ------------------------------------------------------------
#ifndef G4MuonDecayChannel_h
#define G4MuonDecayChannel_h 1
#include "G4ios.hh"
#include "globals.hh"
#include "G4VDecayChannel.hh"
class G4MuonDecayChannel :public G4VDecayChannel
{
// Class Decription
// This class describes muon decay kinemtics.
// This version neglects muon polarization
// assumes the pure V-A coupling
// gives incorrect energy spectrum for neutrinos
//
public: // With Description
//Constructors
G4MuonDecayChannel(const G4String& theParentName,
G4double theBR);
// Destructor
virtual ~G4MuonDecayChannel();
public: // With Description
virtual G4DecayProducts *DecayIt(G4double);
};
#endif

264
geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc Normal file
View File

@ -0,0 +1,264 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History:
// 17 August 2004 P.Gumplinger and T.MacPhail
// samples Michel spectrum including 1st order
// radiative corrections
// Reference: Florian Scheck "Muon Physics", in Physics Reports
// (Review Section of Physics Letters) 44, No. 4 (1978)
// 187-248. North-Holland Publishing Company, Amsterdam
// at page 210 cc.
//
// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
//
// ------------------------------------------------------------
//
#include "G4MuonDecayChannelWithSpin.hh"
#include "Randomize.hh"
#include "G4DecayProducts.hh"
#include "G4LorentzVector.hh"
G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4String& theParentName,
G4double theBR)
: G4MuonDecayChannel(theParentName,theBR)
{
}
G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin()
{
}
G4DecayProducts *G4MuonDecayChannelWithSpin::DecayIt(G4double)
{
// This version assumes V-A coupling with 1st order radiative correctons,
// the standard model Michel parameter values, but
// gives incorrect energy spectrum for neutrinos
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
#endif
if (parent == 0) FillParent();
if (daughters == 0) FillDaughters();
// parent mass
G4double parentmass = parent->GetPDGMass();
EMMU = parentmass;
//daughters'mass
G4double daughtermass[3];
G4double sumofdaughtermass = 0.0;
for (G4int index=0; index<3; index++){
daughtermass[index] = daughters[index]->GetPDGMass();
sumofdaughtermass += daughtermass[index];
}
EMASS = daughtermass[0];
//create parent G4DynamicParticle at rest
G4ThreeVector dummy;
G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
//create G4Decayproducts
G4DecayProducts *products = new G4DecayProducts(*parentparticle);
delete parentparticle;
// calcurate electron energy
G4double michel_rho = 0.75; //Standard Model Michel rho
G4double michel_delta = 0.75; //Standard Model Michel delta
G4double michel_xsi = 1.00; //Standard Model Michel xsi
G4double michel_eta = 0.00; //Standard Model eta
G4double rndm, x, ctheta;
G4double FG;
G4double FG_max = 2.00;
G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
G4double x0 = EMASS/W_mue;
G4double x0_squared = x0*x0;
// ***************************************************
// x0 <= x <= 1. and -1 <= y <= 1
//
// F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y
// ***************************************************
// ***** sampling F(x,y) directly (brute force) *****
do{
// Sample the positron energy by sampling from F
rndm = G4UniformRand();
x = x0 + rndm*(1.-x0);
G4double x_squared = x*x;
G4double F_IS, F_AS, G_IS, G_AS;
F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
G_IS = G_IS + michel_eta*(1.-x)*x0;
G_AS = 3.*(michel_xsi-1.)*(1.-x);
G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
F_IS = F_IS + G_IS;
F_AS = F_AS + G_AS;
// *** Radiative Corrections ***
G4double R_IS = F_c(x,x0);
G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
// *** Radiative Corrections ***
G4double R_AS = F_theta(x,x0);
rndm = G4UniformRand();
ctheta = 2.*rndm-1.;
G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
if(FG>FG_max){
G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
FG_max = FG;
}
rndm = G4UniformRand();
}while(FG<rndm*FG_max);
G4double energy = x * W_mue;
rndm = G4UniformRand();
G4double phi = twopi * rndm;
if(energy < EMASS) energy = EMASS;
// calculate daughter momentum
G4double daughtermomentum[3];
daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
G4double stheta = std::sqrt(1.-ctheta*ctheta);
G4double cphi = std::cos(phi);
G4double sphi = std::sin(phi);
//Coordinates of the decay positron with respect to the muon spin
G4double px = stheta*cphi;
G4double py = stheta*sphi;
G4double pz = ctheta;
G4ThreeVector direction0(px,py,pz);
direction0.rotateUz(parent_polarization);
G4DynamicParticle * daughterparticle0
= new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
products->PushProducts(daughterparticle0);
// daughter 1 ,2 (neutrinos)
// create neutrinos in the C.M frame of two neutrinos
G4double energy2 = parentmass*(1.0 - x/2.0);
G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
G4double beta = -1.0*daughtermomentum[0]/energy2;
G4double costhetan = 2.*G4UniformRand()-1.0;
G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
G4double phin = twopi*G4UniformRand()*rad;
G4double sinphin = std::sin(phin);
G4double cosphin = std::cos(phin);
G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
G4DynamicParticle * daughterparticle1
= new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
G4DynamicParticle * daughterparticle2
= new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
// boost to the muon rest frame
G4LorentzVector p4;
p4 = daughterparticle1->Get4Momentum();
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
daughterparticle1->Set4Momentum(p4);
p4 = daughterparticle2->Get4Momentum();
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
daughterparticle2->Set4Momentum(p4);
products->PushProducts(daughterparticle1);
products->PushProducts(daughterparticle2);
daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
// output message
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) {
G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
G4cout << " create decay products in rest frame " <<G4endl;
products->DumpInfo();
}
#endif
return products;
}
G4double G4MuonDecayChannelWithSpin::R_c(G4double x){
G4int n_max = (int)(100.*x);
if(n_max<10)n_max=10;
G4double L2 = 0.0;
for(G4int n=1; n<=n_max; n++){
L2 += std::pow(x,n)/(n*n);
}
G4double omega = std::log(EMMU/EMASS);
G4double r_c;
r_c = 2.*L2-(pi*pi/3.)-2.;
r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
return r_c;
}

122
geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.hh Normal file
View File

@ -0,0 +1,122 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History:
// 17 August 2004 P.Gumplinger and T.MacPhail
// samples Michel spectrum including 1st order
// radiative corrections
// Reference: Florian Scheck "Muon Physics", in Physics Reports
// (Review Section of Physics Letters) 44, No. 4 (1978)
// 187-248. North-Holland Publishing Company, Amsterdam
// at page 210 cc.
//
// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
//
// ------------------------------------------------------------
#ifndef G4MuonDecayChannelWithSpin_hh
#define G4MuonDecayChannelWithSpin_hh 1
#include "globals.hh"
#include "G4ThreeVector.hh"
#include "G4MuonDecayChannel.hh"
class G4MuonDecayChannelWithSpin : public G4MuonDecayChannel
{
// Class Decription
// This class describes muon decay kinemtics.
// This version assumes V-A coupling with 1st order radiative correctons,
// the standard model Michel parameter values, but
// gives incorrect energy spectrum for neutrinos
public: // With Description
//Constructors
G4MuonDecayChannelWithSpin(const G4String& theParentName,
G4double theBR);
// Destructor
virtual ~G4MuonDecayChannelWithSpin();
public: // With Description
virtual G4DecayProducts *DecayIt(G4double);
void SetPolarization(G4ThreeVector);
const G4ThreeVector& GetPolarization() const;
private:
G4ThreeVector parent_polarization;
// Radiative Correction Factors
G4double F_c(G4double x, G4double x0);
G4double F_theta(G4double x, G4double x0);
G4double R_c(G4double x);
G4double EMMU;
G4double EMASS;
};
inline void G4MuonDecayChannelWithSpin::SetPolarization(G4ThreeVector polar)
{
parent_polarization = polar;
}
inline const G4ThreeVector& G4MuonDecayChannelWithSpin::GetPolarization() const
{
return parent_polarization;
}
inline G4double G4MuonDecayChannelWithSpin::F_c(G4double x, G4double x0)
{
G4double omega = std::log(EMMU/EMASS);
G4double f_c;
f_c = (5.+17.*x-34.*x*x)*(omega+std::log(x))-22.*x+34.*x*x;
f_c = (1.-x)/(3.*x*x)*f_c;
f_c = (6.-4.*x)*R_c(x)+(6.-6.*x)*std::log(x) + f_c;
f_c = (fine_structure_const/twopi) * (x*x-x0*x0) * f_c;
return f_c;
}
inline G4double G4MuonDecayChannelWithSpin::F_theta(G4double x, G4double x0)
{
G4double omega = std::log(EMMU/EMASS);
G4double f_theta;
f_theta = (1.+x+34*x*x)*(omega+std::log(x))+3.-7.*x-32.*x*x;
f_theta = f_theta + ((4.*(1.-x)*(1.-x))/x)*std::log(1.-x);
f_theta = (1.-x)/(3.*x*x) * f_theta;
f_theta = (2.-4.*x)*R_c(x)+(2.-6.*x)*std::log(x)-f_theta;
f_theta = (fine_structure_const/twopi) * (x*x-x0*x0) * f_theta;
return f_theta;
}
#endif

119
geant4/LEMuSR/G4Modified/G4Muonium.cc Normal file
View File

@ -0,0 +1,119 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4Muonium.cc,v 1.8 2003/06/16 16:57:55 gunter Exp $
// GEANT4 tag $Name: geant4-07-00-cand-01 $
//
//
// ----------------------------------------------------------------------
// GEANT 4 class implementation file
//
// History: first implementation, based on object model of
// 4th April 1996, G.Cosmo
// **********************************************************************
// Added particle definitions, H.Kurashige, 19 April 1996
// Added SetCuts implementation, L.Urban, 10 May 1996
// Operators (+=, *=, ++, -> etc.) correctly used, P. Urban, 26/6/96
// Add MuonPlusDefinition(), H.Kurashige 4 July 1996
// Add MuoniumDefinition(),T.K.Paraiso 7 April 2005
// ----------------------------------------------------------------------
#include <fstream>
#include <iomanip>
#include "G4Muonium.hh"
#include "G4MuonDecayChannel.hh"
#include "G4DecayTable.hh"
#include "Randomize.hh"
// ######################################################################
// ### MUONPLUS ###
// ######################################################################
G4Muonium::G4Muonium(
const G4String& aName, G4double mass,
G4double width, G4double charge,
G4int iSpin, G4int iParity,
G4int iConjugation, G4int iIsospin,
G4int iIsospin3, G4int gParity,
const G4String& pType, G4int lepton,
G4int baryon, G4int encoding,
G4bool stable, G4double lifetime,
G4DecayTable *decaytable )
: G4ParticleDefinition( aName,mass,width,charge,iSpin,iParity,
iConjugation,iIsospin,iIsospin3,gParity,pType,
lepton,baryon,encoding,stable,lifetime,decaytable )
{
SetParticleSubType("mu");
SetPDGStable(false);
SetMuoniumGamma();
//create Decay Table
G4DecayTable* table = GetDecayTable();
if (table!=NULL) delete table;
table = new G4DecayTable();
SetDecayTable(table);
}
// ......................................................................
// ... static member definitions ...
// ......................................................................
//
// Arguments for constructor are as follows
// name mass width charge
// 2*spin parity C-conjugation
// 2*Isospin 2*Isospin3 G-parity
// type lepton number baryon number PDG encoding
// stable lifetime decay table
G4Muonium G4Muonium::theMuonium(
"Mu", 0.1056584*GeV, 2.99591e-16*MeV, 0.*eplus,
1, 0, 0,
0, 0, 0,
"lepton", -1, 0, -1313,
true, 2197.03*ns, NULL
);
G4Muonium* G4Muonium::MuoniumDefinition() {return &theMuonium;}
G4Muonium* G4Muonium::Muonium()
{
return &theMuonium;
}
void G4Muonium::SetMuoniumGamma()
{
G4double gamma;
gamma = 0.5*((1.*eplus)/(0.1056584*GeV/(c_light*c_light))-(1.*eplus)/(0.51099906*MeV/(c_light*c_light)));
// SetGyromagneticRatio(gamma);
}

119
geant4/LEMuSR/G4Modified/G4Muonium.cc~ Normal file
View File

@ -0,0 +1,119 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4Muonium.cc,v 1.8 2003/06/16 16:57:55 gunter Exp $
// GEANT4 tag $Name: geant4-07-00-cand-01 $
//
//
// ----------------------------------------------------------------------
// GEANT 4 class implementation file
//
// History: first implementation, based on object model of
// 4th April 1996, G.Cosmo
// **********************************************************************
// Added particle definitions, H.Kurashige, 19 April 1996
// Added SetCuts implementation, L.Urban, 10 May 1996
// Operators (+=, *=, ++, -> etc.) correctly used, P. Urban, 26/6/96
// Add MuonPlusDefinition(), H.Kurashige 4 July 1996
// Add MuoniumDefinition(),T.K.Paraiso 7 April 2005
// ----------------------------------------------------------------------
#include <fstream>
#include <iomanip>
#include "G4Muonium.hh"
#include "G4MuonDecayChannel.hh"
#include "G4DecayTable.hh"
#include "Randomize.hh"
// ######################################################################
// ### MUONPLUS ###
// ######################################################################
G4Muonium::G4Muonium(
const G4String& aName, G4double mass,
G4double width, G4double charge,
G4int iSpin, G4int iParity,
G4int iConjugation, G4int iIsospin,
G4int iIsospin3, G4int gParity,
const G4String& pType, G4int lepton,
G4int baryon, G4int encoding,
G4bool stable, G4double lifetime,
G4DecayTable *decaytable )
: G4ParticleDefinition( aName,mass,width,charge,iSpin,iParity,
iConjugation,iIsospin,iIsospin3,gParity,pType,
lepton,baryon,encoding,stable,lifetime,decaytable )
{
SetParticleSubType("mu");
SetPDGStable(false);
SetMuoniumGamma();
//create Decay Table
G4DecayTable* table = GetDecayTable();
if (table!=NULL) delete table;
table = new G4DecayTable();
SetDecayTable(table);
}
// ......................................................................
// ... static member definitions ...
// ......................................................................
//
// Arguments for constructor are as follows
// name mass width charge
// 2*spin parity C-conjugation
// 2*Isospin 2*Isospin3 G-parity
// type lepton number baryon number PDG encoding
// stable lifetime decay table
G4Muonium G4Muonium::theMuonium(
"Mu", 0.1056584*GeV, 2.99591e-16*MeV, 0.*eplus,
1, 0, 0,
0, 0, 0,
"lepton", -1, 0, -1313,
true, 2197.03*ns, NULL
);
G4Muonium* G4Muonium::MuoniumDefinition() {return &theMuonium;}
G4Muonium* G4Muonium::Muonium()
{
return &theMuonium;
}
void G4Muonium::SetMuoniumGamma()
{
G4double gamma;
gamma = 0.5*((1.*eplus)/(0.1056584*GeV/(c_light*c_light))-(1.*eplus)/(0.51099906*MeV/(c_light*c_light)));
SetGyromagneticRatio(gamma);
}

90
geant4/LEMuSR/G4Modified/G4Muonium.hh Normal file
View File

@ -0,0 +1,90 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4Muonium.hh,v 1.7 2001/10/16 08:16:15 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-cand-01 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 4-th April 1996, G.Cosmo
// ****************************************************************
// Added particle definitions, H.Kurashige, 19 April 1996
// Added SetCuts, L.Urban, 12 June 1996
// Added not static GetEnergyCuts() and GetLengthCuts(), G.Cosmo, 11 July 1996
// Muonium defined as a Mu, T.K.Paraiso 07-04-05
// ----------------------------------------------------------------
// Each class inheriting from G4VLepton
// corresponds to a particle type; one and only one
// instance for each class is guaranteed.
#ifndef G4Muonium_h
#define G4Muonium_h 1
#include "globals.hh"
#include "G4ios.hh"
#include "G4ParticleDefinition.hh"
// ######################################################################
// ### MUONPLUS ###
// ######################################################################
class G4Muonium : public G4ParticleDefinition
{
private: // constructors are hide as private
G4Muonium(
const G4String& aName, G4double mass,
G4double width, G4double charge,
G4int iSpin, G4int iParity,
G4int iConjugation, G4int iIsospin,
G4int iIsospin3, G4int gParity,
const G4String& pType, G4int lepton,
G4int baryon, G4int encoding,
G4bool stable, G4double lifetime,
G4DecayTable *decaytable
);
public:
virtual ~G4Muonium(){;};
void SetMuoniumGamma();
static G4Muonium* MuoniumDefinition();
static G4Muonium* Muonium();
private:
static G4Muonium theMuonium;
};
#endif

90
geant4/LEMuSR/G4Modified/G4Muonium.hh~ Normal file
View File

@ -0,0 +1,90 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4Muonium.hh,v 1.7 2001/10/16 08:16:15 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-cand-01 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 4-th April 1996, G.Cosmo
// ****************************************************************
// Added particle definitions, H.Kurashige, 19 April 1996
// Added SetCuts, L.Urban, 12 June 1996
// Added not static GetEnergyCuts() and GetLengthCuts(), G.Cosmo, 11 July 1996
// Muonium defined as a Mu0, T.K.Paraiso 07-04-05
// ----------------------------------------------------------------
// Each class inheriting from G4VLepton
// corresponds to a particle type; one and only one
// instance for each class is guaranteed.
#ifndef G4Muonium_h
#define G4Muonium_h 1
#include "globals.hh"
#include "G4ios.hh"
#include "G4ParticleDefinition.hh"
// ######################################################################
// ### MUONPLUS ###
// ######################################################################
class G4Muonium : public G4ParticleDefinition
{
private:
static G4Muonium theMuonium;
private: // constructors are hide as private
G4Muonium(
const G4String& aName, G4double mass,
G4double width, G4double charge,
G4int iSpin, G4int iParity,
G4int iConjugation, G4int iIsospin,
G4int iIsospin3, G4int gParity,
const G4String& pType, G4int lepton,
G4int baryon, G4int encoding,
G4bool stable, G4double lifetime,
G4DecayTable *decaytable
);
public:
virtual ~G4Muonium(){}
void SetMuoniumGamma();
static G4Muonium* MuoniumDefinition();
static G4Muonium* Muonium();
};
#endif

290
geant4/LEMuSR/G4Modified/G4ParticleDefinition.cc Normal file
View File

@ -0,0 +1,290 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4ParticleDefinition.cc,v 1.23 2005/01/14 03:00:39 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-patch-01 $
//
//
// --------------------------------------------------------------
// GEANT 4 class implementation file
//
// History: first implementation, based on object model of
// 2nd December 1995, G.Cosmo
// ---------------- G4ParticleDefinition -----------------
// first implementation by Makoto Asai, 29 January 1996
// revised by G.Cosmo, 29 February 1996
// revised by H.Kurashige, 19 April 1996
// Code uses operators (+=, *=, ++, -> etc.) correctly, P. Urban, 26/6/96
// revised by H.Kurashige, 4 July 1996
// revised by H.Kurashige, 16 Feb 1997
// revised by H.Kurashige, 10 Nov 1997
// remove new/delete G4ProcessManager by H.Kurashige 06 June 1998
// added Resonance flag and ApplyCuts flag H.Kurashige 27 June 1998
// modify FillQuarkContents() for quarks/diquarks H.Kurashige 30 June 1998
// modify encoding rule H.Kurashige 23 Oct. 98
// modify FillQuarkContents() for deltas 25 Nov.,98 H.Kurashige
//
// modify FillQuarkContents() to use G4PDGCodeChecker 17 Aug. 99 H.Kurashige
// added Gyromagnetic Factors T.K.Paraïso Okt. 05
// --------------------------------------------------------------
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
#include "G4DecayTable.hh"
#include "G4PDGCodeChecker.hh"
G4ParticleDefinition::G4ParticleDefinition(
const G4String& aName,
G4double mass,
G4double width,
G4double charge,
G4int iSpin,
G4int iParity,
G4int iConjugation,
G4int iIsospin,
G4int iIsospin3,
G4int gParity,
const G4String& pType,
G4int lepton,
G4int baryon,
G4int encoding,
G4bool stable,
G4double lifetime,
G4DecayTable *decaytable,
G4bool shortlived,
const G4String& subType,
G4int anti_encoding,
G4double gFactor)
: theParticleName(aName),
thePDGMass(mass),
thePDGWidth(width),
thePDGCharge(charge),
thePDGiSpin(iSpin),
thePDGSpin(iSpin*0.5),
thePDGiParity(iParity),
thePDGiConjugation(iConjugation),
thePDGiGParity(gParity),
thePDGiIsospin(iIsospin),
thePDGiIsospin3(iIsospin3),
thePDGIsospin(iIsospin*0.5),
thePDGIsospin3(iIsospin3*0.5),
theLeptonNumber(lepton),
theBaryonNumber(baryon),
theParticleType(pType),
theParticleSubType(subType),
thePDGEncoding(encoding),
theAntiPDGEncoding(-1*encoding),
fShortLivedFlag(shortlived),
thePDGStable(stable),
thePDGLifeTime(lifetime),
theDecayTable(decaytable),
theProcessManager(0),
theAtomicNumber(0),
theAtomicMass(0),
theGyromagneticFactor(gFactor),//absolute value tao
theGyromagneticRatio(gFactor*charge/(2*mass/(c_light*c_light))),
verboseLevel(1),
fApplyCutsFlag(false)
{
// tao
// check name and register this particle into ParticleTable
theParticleTable = G4ParticleTable::GetParticleTable();
theParticleTable->Insert(this);
if (anti_encoding !=0) theAntiPDGEncoding = anti_encoding;
// check quark contents
#ifdef G4VERBOSE
if (this->FillQuarkContents() != thePDGEncoding) {
if (verboseLevel>0) {
// cerr bnot G4cerr is used intentionally
// because G4ParticleDefinition constructor may be called
// before G4cerr object is instantiated !!
G4cerr << "Particle " << aName << " has a strange PDGEncoding " <<G4endl;
}
}
#endif
}
G4ParticleDefinition::G4ParticleDefinition(const G4ParticleDefinition &)
{
G4Exception("You call Copy Constructor of G4ParticleDefinition ");
}
G4ParticleDefinition::G4ParticleDefinition()
{
G4Exception("You call Default Constructor of G4ParticleDefinition ");
}
G4ParticleDefinition::~G4ParticleDefinition()
{
if (theDecayTable!= 0) delete theDecayTable;
}
const G4ParticleDefinition & G4ParticleDefinition::operator=(const G4ParticleDefinition &right)
{
if (this != &right) {
}
return *this;
}
G4int G4ParticleDefinition::operator==(const G4ParticleDefinition &right) const
{
return (this->theParticleName == right.theParticleName);
}
G4int G4ParticleDefinition::operator!=(const G4ParticleDefinition &right) const
{
return (this->theParticleName != right.theParticleName);
}
G4int G4ParticleDefinition::FillQuarkContents()
// calculate quark and anti-quark contents
// return value is PDG encoding for this particle.
// It means error if the return value is differnt from
// this->thePDGEncoding.
{
G4int flavor;
for (flavor= 0; flavor<NumberOfQuarkFlavor; flavor++){
theQuarkContent[flavor] = 0;
theAntiQuarkContent[flavor] = 0;
}
G4PDGCodeChecker checker;
G4int temp = checker.CheckPDGCode(thePDGEncoding, theParticleType);
if ( temp != 0) {
for (flavor= 0; flavor<NumberOfQuarkFlavor; flavor++){
theQuarkContent[flavor] = checker.GetQuarkContent(flavor);
theAntiQuarkContent[flavor] = checker.GetAntiQuarkContent(flavor);
}
if ((theParticleType == "meson")||(theParticleType == "baryon")) {
// check charge
if (!checker.CheckCharge(thePDGCharge) ){
temp = 0;
#ifdef G4VERBOSE
if (verboseLevel>1) {
G4cout << " illegal charge (" << thePDGCharge/eplus;
G4cout << " PDG code=" << thePDGEncoding <<G4endl;
}
#endif
}
// check spin
if (checker.GetSpin() != thePDGiSpin) {
temp=0;
#ifdef G4VERBOSE
if (verboseLevel>1) {
G4cout << " illegal SPIN (" << thePDGiSpin << "/2";
G4cout << " PDG code=" << thePDGEncoding <<G4endl;
}
#endif
}
}
}
return temp;
}
void G4ParticleDefinition::DumpTable() const
{
G4cout << G4endl;
G4cout << "--- G4ParticleDefinition ---" << G4endl;
G4cout << " Particle Name : " << theParticleName << G4endl;
G4cout << " PDG particle code : " << thePDGEncoding;
G4cout << " [PDG anti-particle code: " << this->GetAntiPDGEncoding() << "]"<< G4endl;
G4cout << " Mass [GeV/c2] : " << thePDGMass/GeV ;
G4cout << " Width : " << thePDGWidth/GeV << G4endl;
G4cout << " Lifetime [nsec] : " << thePDGLifeTime/ns << G4endl;
G4cout << " Charge [e]: " << thePDGCharge/eplus << G4endl;
G4cout << " Spin : " << thePDGiSpin << "/2" << G4endl;
G4cout << " Parity : " << thePDGiParity << G4endl;
G4cout << " Charge conjugation : " << thePDGiConjugation << G4endl;
G4cout << " Isospin : (I,Iz): (" << thePDGiIsospin <<"/2";
G4cout << " , " << thePDGiIsospin3 << "/2 ) " << G4endl;
G4cout << " GParity : " << thePDGiGParity << G4endl;
G4cout << " Quark contents (d,u,s,c,b,t) : " << theQuarkContent[0];
G4cout << ", " << theQuarkContent[1];
G4cout << ", " << theQuarkContent[2];
G4cout << ", " << theQuarkContent[3];
G4cout << ", " << theQuarkContent[4];
G4cout << ", " << theQuarkContent[5] << G4endl;
G4cout << " AntiQuark contents : " << theAntiQuarkContent[0];
G4cout << ", " << theAntiQuarkContent[1];
G4cout << ", " << theAntiQuarkContent[2];
G4cout << ", " << theAntiQuarkContent[3];
G4cout << ", " << theAntiQuarkContent[4];
G4cout << ", " << theAntiQuarkContent[5] << G4endl;
G4cout << " Lepton number : " << theLeptonNumber;
G4cout << " Baryon number : " << theBaryonNumber << G4endl;
G4cout << " Particle type : " << theParticleType ;
G4cout << " [" << theParticleSubType << "]" << G4endl;
if ( fShortLivedFlag ){
G4cout << " ShortLived : ON" << G4endl;
}
if ( thePDGStable ){
G4cout << " Stable : stable" << G4endl;
} else {
if( theDecayTable != 0 ){
theDecayTable->DumpInfo();
} else {
G4cout << "Decay Table is not defined !!" <<G4endl;
}
}
}
void G4ParticleDefinition::SetApplyCutsFlag(G4bool flg)
{
if(theParticleName=="gamma"
|| theParticleName=="e-"
|| theParticleName=="e+")
{ fApplyCutsFlag = flg; }
else
{
G4cerr
<< "G4ParticleDefinition::SetApplyCutsFlag() for " << theParticleName
<< G4endl;
G4cerr
<< "becomes obsolete. Production threshold is applied only for "
<< "gamma, e- and e+." << G4endl;
}
}

334
geant4/LEMuSR/G4Modified/G4ParticleDefinition.hh Normal file
View File

@ -0,0 +1,334 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4ParticleDefinition.hh,v 1.27 2005/01/14 03:00:38 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-patch-01 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 2nd December 1995, G.Cosmo
// ---------------- G4ParticleDefinition ----------------
// first implementation by Makoto Asai, 29 January 1996
// revised by G.Cosmo, 29 February 1996
// revised by H.Kurashige, 19 April 1996
// revised by H.Kurashige, 4 July 1996
// added GetEnergyCuts() and GetLengthCuts() by G.Cosmo, 11 July 1996
// added Set/GetVerboseLevel() H.Kurashige 11 Nov. 1997
// added SetCuts() and ResetCuts H.Kurashige 15 Nov.1996
// change SetProcessManager as public H.Kurashige 06 June 1998
// added GetEnergyThreshold H.Kurashige 08 June 1998
// added ShortLived flag and ApplyCuts flag H.Kurashige 27 June 1998
// fixed some improper codings H.Kurashige 08 Apr. 1999
// added sub-type H.Kurashige 15 Feb. 2000
// added RestoreCuts H.Kurashige 09 Mar. 2001
// restructuring for Cuts per Region by Hisaya 11 MAr.2003
// added Gyromagnetic Factors T.K.Paraïso Okt. 05
// ------------------------------------------------------------
#ifndef G4ParticleDefinition_h
#define G4ParticleDefinition_h 1
#include "globals.hh"
#include "G4ios.hh"
#include <vector>
class G4ProcessManager;
class G4DecayTable;
class G4ParticleTable;
class G4ParticlePropertyTable;
class G4ParticleDefinition
{
// Class Description
// This class containes all the static data of a particle.
// It also has uses a process manager in order to collect
// all the processes this kind of particle can undertake.
//
friend class G4ParticlePropertyTable;
public: // With Description
// Only one type of constructor can be used for G4ParticleDefinition.
// If you want to create new particle, you must set name of the particle
// at construction. Most of members seen as arguments of the constructor
// (except last 3 arguments concerning with decay ) are "constant"
// and can not be changed later. (No "SET" methods are available)
// Each type of particle must be constructed as a unique static object
// of special class derived from G4ParticleDefinition.
// see G4ParticleTypes for detail
G4ParticleDefinition(const G4String& aName,
G4double mass,
G4double width,
G4double charge,
G4int iSpin,
G4int iParity,
G4int iConjugation,
G4int iIsospin,
G4int iIsospinZ,
G4int gParity,
const G4String& pType,
G4int lepton,
G4int baryon,
G4int encoding,
G4bool stable,
G4double lifetime,
G4DecayTable *decaytable,
G4bool shortlived = false,
const G4String& subType ="",
G4int anti_encoding =0,
G4double gFactor=2.0023193);
virtual ~G4ParticleDefinition();
public: // With Description
// By these following Getxxxx methods, you can get values
// for members which can not be changed
//
const G4String& GetParticleName() const { return theParticleName; }
G4double GetPDGMass() const { return thePDGMass; }
G4double GetPDGWidth() const { return thePDGWidth; }
G4double GetPDGCharge() const { return thePDGCharge; }
G4double GetPDGSpin() const { return thePDGSpin; }
G4int GetPDGiSpin() const { return thePDGiSpin; }
G4int GetPDGiParity() const { return thePDGiParity; }
G4int GetPDGiConjugation() const { return thePDGiConjugation; }
G4double GetPDGIsospin() const { return thePDGIsospin; }
G4double GetPDGIsospin3() const { return thePDGIsospin3; }
G4int GetPDGiIsospin() const { return thePDGiIsospin; }
G4int GetPDGiIsospin3() const { return thePDGiIsospin3; }
G4int GetPDGiGParity() const { return thePDGiGParity; }
const G4String& GetParticleType() const { return theParticleType; }
const G4String& GetParticleSubType() const { return theParticleSubType; }
G4int GetLeptonNumber() const { return theLeptonNumber; }
G4int GetBaryonNumber() const { return theBaryonNumber; }
G4int GetPDGEncoding() const { return thePDGEncoding; }
G4int GetAntiPDGEncoding() const { return theAntiPDGEncoding; }
void SetAntiPDGEncoding(G4int aEncoding);
G4int GetQuarkContent(G4int flavor) const;
G4int GetAntiQuarkContent(G4int flavor) const;
// return the number of quark with flavor contained in this particle.
// The value of flavor is assigned as follows
// 1:d, 2:u, 3:s, 4:c, 5:b, 6:t
public: // With Description
// ShortLived flag
G4bool IsShortLived() const { return fShortLivedFlag; }
G4bool GetPDGStable() const { return thePDGStable; }
void SetPDGStable(const G4bool aFlag) { thePDGStable=aFlag; }
G4double GetPDGLifeTime() const { return thePDGLifeTime; }
void SetPDGLifeTime(G4double aLifeTime) { thePDGLifeTime = aLifeTime; }
public:// With Description
G4DecayTable* GetDecayTable();
void SetDecayTable(G4DecayTable* aDecayTable);
// Set/Get Decay Table
// !! Decay Table can be modified !!
public: // With Description
G4ProcessManager* GetProcessManager() const;
void SetProcessManager(G4ProcessManager* aProcessManager);
// Set/Get Process Manager
// !! Process Manager can be modified !!
G4ParticleTable* GetParticleTable();
// get pointer to the particle table
void DumpTable() const;
// Prints information of data members.
protected:
G4int FillQuarkContents();
// calculate quark and anti-quark contents
// return value is PDG encoding for this particle.
// It means error if the return value is deffernt from
// this->thePDGEncoding.
void SetParticleSubType(const G4String& subtype);
public:
void SetAtomicNumber(G4int);
G4int GetAtomicNumber() const;
void SetAtomicMass(G4int);
G4int GetAtomicMass() const;
G4double GetGyromagneticFactor() const { return theGyromagneticFactor; }
void SetGyromagneticFactor(G4double ge) { theGyromagneticFactor = ge; }
G4double GetGyromagneticRatio() const { return theGyromagneticRatio; }
void SetGyromagneticRatio(G4double gamma) { theGyromagneticRatio = gamma; }
public:
void SetVerboseLevel(G4int value);
G4int GetVerboseLevel() const;
// controle flag for output message
// 0: Silent
// 1: Warning message
// 2: More
protected:
// !!! can not use "copy constructor" nor "default constructor" !!!!
G4ParticleDefinition(const G4ParticleDefinition &right);
G4ParticleDefinition();
private:
// !!! Assignment operation is forbidden !!!
const G4ParticleDefinition & operator=(const G4ParticleDefinition &right);
public:
G4int operator==(const G4ParticleDefinition &right) const;
G4int operator!=(const G4ParticleDefinition &right) const;
private:
// Values following can not be changed
// i.e. No Setxxxx Methods for them
G4String theParticleName;
// The name of the particle.
// Each object must have its specific name!!
// --- following member values must be defined with Units
G4double thePDGMass;
// The mass of the particle, in units of equivalent energy.
G4double thePDGWidth;
// The decay width of the particle, usually the width of a
// Breit-Wigner function, assuming that you are near the
// mass center anyway. (in units of equivalent energy)
G4double thePDGCharge;
// The charge of the particle.(in units of Coulomb)
// ---- following members are quantum number
// i.e. discrete numbers can be allowded
// So, you can defined only by using integer in constructor
G4int thePDGiSpin;
// The total spin of the particle, also often denoted as
// capital J, in units of 1/2.
G4double thePDGSpin;
// The total spin of the particle, in units of 1.
G4int thePDGiParity;
// The parity quantum number, in units of 1. If the parity
// is not defined for this particle, we will set this to 0.
G4int thePDGiConjugation;
// This charge conjugation quantum number in units of 1.
G4int thePDGiGParity;
// The value of the G-parity quantum number.
G4int thePDGiIsospin;
G4int thePDGiIsospin3;
// The isospin and its 3rd-component in units of 1/2.
G4double thePDGIsospin;
G4double thePDGIsospin3;
// The isospin quantum number in units of 1.
G4int theLeptonNumber;
// The lepton quantum number.
G4int theBaryonNumber;
// The baryon quantum number.
G4String theParticleType;
// More general textual type description of the particle.
G4String theParticleSubType;
// Textual type description of the particle
// eg. pion, lamda etc.
G4int thePDGEncoding;
// The Particle Data Group integer identifier of this particle
G4int theAntiPDGEncoding;
// The Particle Data Group integer identifier of the anti-particle
protected:
enum {NumberOfQuarkFlavor = 6};
G4int theQuarkContent[NumberOfQuarkFlavor];
G4int theAntiQuarkContent[NumberOfQuarkFlavor];
// the number of quark (minus Sign means anti-quark) contents
// The value of flavor is assigned as follows
// 0:d, 1:u, 2:s, 3:c, 4:b, 5:t
private:
// Following members can be changed after construction
G4bool fShortLivedFlag;
// Particles which have true value of this flag
// will not be tracked by TrackingManager
G4bool thePDGStable;
// Is an indicator that this particle is stable. It must
// not decay. If the user tries to assign a kind of decay
// object to it, it will refuse to take it.
G4double thePDGLifeTime;
// Is related to the decay width of the particle. The mean
// life time is given in seconds.
class G4DecayTable *theDecayTable;
// Points DecayTable
private:
class G4ProcessManager *theProcessManager;
// Points to G4ProcessManager
G4ParticleTable* theParticleTable;
private:
G4int theAtomicNumber;
G4int theAtomicMass;
G4double theGyromagneticFactor;
// The gyromanetic factor of the particle
G4double theGyromagneticRatio;
// The gyromanetic ratio of the particle
private:
G4int verboseLevel;
private:
G4bool fApplyCutsFlag;
public:
void SetApplyCutsFlag(G4bool);
G4bool GetApplyCutsFlag() const;
};
#include "G4ParticleDefinition.icc"
#endif

122
geant4/LEMuSR/G4Modified/G4ParticleDefinition.icc Normal file
View File

@ -0,0 +1,122 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4ParticleDefinition.icc,v 1.8 2004/12/02 08:08:58 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
//
inline
G4ParticleTable* G4ParticleDefinition::GetParticleTable()
{
return theParticleTable;
}
inline
G4DecayTable* G4ParticleDefinition::GetDecayTable()
{
return theDecayTable;
}
inline
void G4ParticleDefinition::SetDecayTable(G4DecayTable* aDecayTable)
{
theDecayTable = aDecayTable;
}
inline
void G4ParticleDefinition::SetVerboseLevel(G4int value)
{
verboseLevel = value;
}
inline
G4int G4ParticleDefinition::GetVerboseLevel() const
{
return verboseLevel;
}
inline
G4ProcessManager* G4ParticleDefinition::GetProcessManager() const
{
return theProcessManager;
}
inline
void G4ParticleDefinition::SetProcessManager(G4ProcessManager *aProcessManager)
{
theProcessManager = aProcessManager;
}
inline
G4int G4ParticleDefinition::GetQuarkContent(G4int flavor) const
{
G4int content = 0;
if ((flavor>0) && (flavor<=NumberOfQuarkFlavor)){
content = theQuarkContent[flavor-1];
}else {
#ifdef G4VERBOSE
if (verboseLevel >0) {
G4cout << "Invalid Quark Flavor for G4ParticleDefinition::GetQuarkContent";
G4cout << ": flavor=" << flavor <<G4endl;
}
#endif
}
return content;
}
inline
G4int G4ParticleDefinition::GetAntiQuarkContent(G4int flavor) const
{
G4int content = 0;
if ((flavor>0) && (flavor<=NumberOfQuarkFlavor)){
content = theAntiQuarkContent[flavor-1];
}else {
#ifdef G4VERBOSE
if (verboseLevel >0) {
G4cout <<"Invalid Quark Flavor for G4ParticleDefinition::GetAntiQuarkContent";
G4cout << ": flavor=" << flavor <<G4endl;
}
#endif
}
return content;
}
inline
void G4ParticleDefinition::SetParticleSubType(const G4String& subtype)
{
theParticleSubType = subtype;
}
inline
void G4ParticleDefinition::SetAntiPDGEncoding(G4int aEncoding)
{
theAntiPDGEncoding = aEncoding;
}
inline
G4bool G4ParticleDefinition::GetApplyCutsFlag() const
{
return fApplyCutsFlag;
}

460
geant4/LEMuSR/G4Modified/G4VDecayChannel.cc Normal file
View File

@ -0,0 +1,460 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4VDecayChannel.cc,v 1.16 2004/12/02 08:09:00 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 27 July 1996 H.Kurashige
// 30 May 1997 H.Kurashige
// 23 Mar. 2000 H.Weber : add GetAngularMomentum
// ------------------------------------------------------------
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
#include "G4DecayTable.hh"
#include "G4DecayProducts.hh"
#include "G4VDecayChannel.hh"
const G4String G4VDecayChannel::noName = " ";
G4VDecayChannel::G4VDecayChannel(const G4String &aName, G4int Verbose)
:kinematics_name(aName),
rbranch(0.0),
numberOfDaughters(0),
parent_name(0), daughters_name(0),
particletable(0),
parent(0), daughters(0),
parent_mass(0.0), daughters_mass(0),
verboseLevel(Verbose)
{
// set pointer to G4ParticleTable (static and singleton object)
particletable = G4ParticleTable::GetParticleTable();
}
G4VDecayChannel::G4VDecayChannel(const G4String &aName,
const G4String& theParentName,
G4double theBR,
G4int theNumberOfDaughters,
const G4String& theDaughterName1,
const G4String& theDaughterName2,
const G4String& theDaughterName3,
const G4String& theDaughterName4 )
:kinematics_name(aName),
rbranch(theBR),
numberOfDaughters(theNumberOfDaughters),
parent_name(0), daughters_name(0),
particletable(0),
parent(0), daughters(0),
parent_mass(0.0), daughters_mass(0),
verboseLevel(1)
{
// set pointer to G4ParticleTable (static and singleton object)
particletable = G4ParticleTable::GetParticleTable();
// parent name
parent_name = new G4String(theParentName);
// cleate array
daughters_name = new G4String*[numberOfDaughters];
for (G4int index=0;index<numberOfDaughters;index++) daughters_name[index]=0;
// daughters' name
if (numberOfDaughters>0) daughters_name[0] = new G4String(theDaughterName1);
if (numberOfDaughters>1) daughters_name[1] = new G4String(theDaughterName2);
if (numberOfDaughters>2) daughters_name[2] = new G4String(theDaughterName3);
if (numberOfDaughters>3) daughters_name[3] = new G4String(theDaughterName4);
}
G4VDecayChannel::G4VDecayChannel(const G4VDecayChannel &right)
{
kinematics_name = right.kinematics_name;
verboseLevel = right.verboseLevel;
rbranch = right.rbranch;
// copy parent name
parent_name = new G4String(*right.parent_name);
parent = 0;
parent_mass = 0.0;
//create array
numberOfDaughters = right.numberOfDaughters;
if ( numberOfDaughters >0 ) {
daughters_name = new G4String*[numberOfDaughters];
//copy daughters name
for (G4int index=0; index < numberOfDaughters; index++)
{
daughters_name[index] = new G4String(*right.daughters_name[index]);
}
}
//
daughters_mass = 0;
daughters = 0;
// particle table
particletable = G4ParticleTable::GetParticleTable();
}
G4VDecayChannel & G4VDecayChannel::operator=(const G4VDecayChannel &right)
{
if (this != &right) {
kinematics_name = right.kinematics_name;
verboseLevel = right.verboseLevel;
rbranch = right.rbranch;
// copy parent name
parent_name = new G4String(*right.parent_name);
// clear daughters_name array
ClearDaughtersName();
// recreate array
numberOfDaughters = right.numberOfDaughters;
if ( numberOfDaughters >0 ) {
daughters_name = new G4String*[numberOfDaughters];
//copy daughters name
for (G4int index=0; index < numberOfDaughters; index++) {
daughters_name[index] = new G4String(*right.daughters_name[index]);
}
}
}
//
parent = 0;
daughters = 0;
parent_mass = 0.0;
daughters_mass = 0;
// particle table
particletable = G4ParticleTable::GetParticleTable();
return *this;
}
G4VDecayChannel::~G4VDecayChannel()
{
if (parent_name != 0) delete parent_name;
ClearDaughtersName();
if (daughters_mass != 0) delete [] daughters_mass;
}
void G4VDecayChannel::ClearDaughtersName()
{
if ( daughters_name != 0) {
if (numberOfDaughters>0) {
#ifdef G4VERBOSE
if (verboseLevel>1) {
G4cout << "G4VDecayChannel::ClearDaughtersName ";
G4cout << "clear all daughters " << G4endl;
}
#endif
for (G4int index=0; index < numberOfDaughters; index++) {
if (daughters_name[index] != 0) delete daughters_name[index];
}
}
delete [] daughters_name;
daughters_name = 0;
}
//
if (daughters != 0) delete [] daughters;
if (daughters_mass != 0) delete [] daughters_mass;
daughters = 0;
daughters_mass = 0;
numberOfDaughters = 0;
}
void G4VDecayChannel::SetNumberOfDaughters(G4int size)
{
if (size >0) {
// remove old contents
ClearDaughtersName();
// cleate array
daughters_name = new G4String*[size];
for (G4int index=0;index<size;index++) daughters_name[index]=0;
numberOfDaughters = size;
}
}
void G4VDecayChannel::SetDaughter(G4int anIndex,
const G4String &particle_name)
{
// check numberOfDaughters is positive
if (numberOfDaughters<=0) {
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::SetDaughter: ";
G4cout << "Number of daughters is not defined" << G4endl;
}
#endif
return;
}
// check existence of daughters_name array
if (daughters_name == 0) {
// cleate array
daughters_name = new G4String*[numberOfDaughters];
for (G4int index=0;index<numberOfDaughters;index++) {
daughters_name[index]=0;
}
}
// check an index
if ( (anIndex<0) || (anIndex>=numberOfDaughters) ) {
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::SetDaughter";
G4cout << "index out of range " << anIndex << G4endl;
}
#endif
} else {
// delete the old name if it exists
if (daughters_name[anIndex]!=0) delete daughters_name[anIndex];
// fill the name
daughters_name[anIndex] = new G4String(particle_name);
// refill the array of daughters[] if it exists
if (daughters != 0) FillDaughters();
#ifdef G4VERBOSE
if (verboseLevel>1) {
G4cout << "G4VDecayChannel::SetDaughter[" << anIndex <<"] :";
G4cout << daughters_name[anIndex] << ":" << *daughters_name[anIndex]<<G4endl;
}
#endif
}
}
void G4VDecayChannel::SetDaughter(G4int anIndex, const G4ParticleDefinition * parent_type)
{
if (parent_type != 0) SetDaughter(anIndex, parent_type->GetParticleName());
}
void G4VDecayChannel::FillDaughters()
{
G4int index;
#ifdef G4VERBOSE
if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <<G4endl;
#endif
if (daughters != 0) delete [] daughters;
// parent mass
if (parent == 0) FillParent();
G4double parentmass = parent->GetPDGMass();
//
G4double sumofdaughtermass = 0.0;
if ((numberOfDaughters <=0) || (daughters_name == 0) ){
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::FillDaughters ";
G4cout << "[ " << parent->GetParticleName() << " ]";
G4cout << "numberOfDaughters is not defined yet";
}
#endif
daughters = 0;
G4Exception("G4VDecayChannel::FillDaughters");
}
//create and set the array of pointers to daughter particles
daughters = new G4ParticleDefinition*[numberOfDaughters];
if (daughters_mass != 0) delete [] daughters_mass;
daughters_mass = new G4double[numberOfDaughters];
// loop over all daughters
for (index=0; index < numberOfDaughters; index++) {
if (daughters_name[index] == 0) {
// daughter name is not defined
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::FillDaughters ";
G4cout << "[ " << parent->GetParticleName() << " ]";
G4cout << index << "-th daughter is not defined yet" << G4endl;
}
#endif
daughters[index] = 0;
G4Exception("G4VDecayChannel::FillDaughters");
}
//search daughter particles in the particle table
daughters[index] = particletable->FindParticle(*daughters_name[index]);
if (daughters[index] == 0) {
// can not find the daughter particle
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::FillDaughters ";
G4cout << "[ " << parent->GetParticleName() << " ]";
G4cout << index << ":" << *daughters_name[index];
G4cout << " is not defined !!" << G4endl;
G4cout << " The BR of this decay mode is set to zero " << G4endl;
}
#endif
SetBR(0.0);
return;
}
#ifdef G4VERBOSE
if (verboseLevel>1) {
G4cout << index << ":" << *daughters_name[index];
G4cout << ":" << daughters[index] << G4endl;
}
#endif
daughters_mass[index] = daughters[index]->GetPDGMass();
sumofdaughtermass += daughters[index]->GetPDGMass();
} // end loop over all daughters
// check sum of daghter mass
G4double widthMass = parent->GetPDGWidth();
if ( (parent->GetParticleType() != "nucleus") &&
(sumofdaughtermass > parentmass + 5*widthMass) ){
// !!! illegal mass !!!
#ifdef G4VERBOSE
if (GetVerboseLevel()>0) {
G4cout << "G4VDecayChannel::FillDaughters ";
G4cout << "[ " << parent->GetParticleName() << " ]";
G4cout << " Energy/Momentum conserevation breaks " <<G4endl;
if (GetVerboseLevel()>1) {
G4cout << " parent:" << *parent_name;
G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
for (index=0; index < numberOfDaughters; index++){
G4cout << " daughter " << index << ":" << *daughters_name[index];
G4cout << " mass:" << daughters[index]->GetPDGMass()/GeV;
G4cout << "[GeV/c/c]" <<G4endl;
}
}
}
#endif
}
}
void G4VDecayChannel::FillParent()
{
if (parent_name == 0) {
// parent name is not defined
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::FillParent ";
G4cout << ": parent name is not defined !!" << G4endl;
}
#endif
parent = 0;
G4Exception("G4VDecayChannel::FillParent");
}
// search parent particle in the particle table
parent = particletable->FindParticle(*parent_name);
if (parent == 0) {
// parent particle does not exist
#ifdef G4VERBOSE
if (verboseLevel>0) {
G4cout << "G4VDecayChannel::FillParent ";
G4cout << *parent_name << " does not exist !!" << G4endl;
}
#endif
G4Exception("G4VDecayChannel::FillParent");
}
parent_mass = parent->GetPDGMass();
}
void G4VDecayChannel::SetParent(const G4ParticleDefinition * parent_type)
{
if (parent_type != 0) SetParent(parent_type->GetParticleName());
}
G4int G4VDecayChannel::GetAngularMomentum()
{
// determine angular momentum
// fill pointers to daughter particles if not yet set
if (daughters == 0) FillDaughters();
const G4int PiSpin = parent->GetPDGiSpin();
const G4int PParity = parent->GetPDGiParity();
if (2==numberOfDaughters) { // up to now we can only handle two particle decays
const G4int D1iSpin = daughters[0]->GetPDGiSpin();
const G4int D1Parity = daughters[0]->GetPDGiParity();
const G4int D2iSpin = daughters[1]->GetPDGiSpin();
const G4int D2Parity = daughters[1]->GetPDGiParity();
const G4int MiniSpin = std::abs (D1iSpin - D2iSpin);
const G4int MaxiSpin = D1iSpin + D2iSpin;
const G4int lMax = (PiSpin+D1iSpin+D2iSpin)/2; // l is allways int
G4int lMin;
#ifdef G4VERBOSE
if (verboseLevel>1) {
G4cout << "iSpin: " << PiSpin << " -> " << D1iSpin << " + " << D2iSpin << G4endl;
G4cout << "2*jmin, 2*jmax, lmax " << MiniSpin << " " << MaxiSpin << " " << lMax << G4endl;
}
#endif
for (G4int j=MiniSpin; j<=MaxiSpin; j+=2){ // loop over all possible spin couplings
lMin = std::abs(PiSpin-j)/2;
#ifdef G4VERBOSE
if (verboseLevel>1)
G4cout << "-> checking 2*j=" << j << G4endl;
#endif
for (G4int l=lMin; l<=lMax; l++) {
#ifdef G4VERBOSE
if (verboseLevel>1)
G4cout << " checking l=" << l << G4endl;
#endif
if (l%2==0) {
if (PParity == D1Parity*D2Parity) { // check parity for this l
return l;
}
} else {
if (PParity == -1*D1Parity*D2Parity) { // check parity for this l
return l;
}
}
}
}
} else {
G4Exception ("G4VDecayChannel::GetAngularMomentum: Sorry, can't handle 3 particle decays (up to now)");
}
G4Exception ("G4VDecayChannel::GetAngularMomentum: Can't find angular momentum for this decay!");
return 0;
}
void G4VDecayChannel::DumpInfo()
{
G4cout << " BR: " << rbranch << " [" << kinematics_name << "]";
G4cout << " : " ;
for (G4int index=0; index < numberOfDaughters; index++)
{
if(daughters_name[index] != 0) {
G4cout << " " << *(daughters_name[index]);
} else {
G4cout << " not defined ";
}
}
G4cout << G4endl;
}
const G4String& G4VDecayChannel::GetNoName() const
{
return noName;
}

288
geant4/LEMuSR/G4Modified/G4VDecayChannel.hh Normal file
View File

@ -0,0 +1,288 @@
//
// ********************************************************************
// * DISCLAIMER *
// * *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in: *
// * http://cern.ch/geant4/license *
// * *
// * 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. *
// * *
// * This code implementation is the intellectual property of the *
// * GEANT4 collaboration. *
// * By copying, distributing or modifying the Program (or any work *
// * based on the Program) you indicate your acceptance of this *
// * statement, and all its terms. *
// ********************************************************************
//
//
// $Id: G4VDecayChannel.hh,v 1.11 2004/12/02 08:08:58 kurasige Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 27 July 1996 H.Kurashige
// 30 May 1997 H.Kurashige
// 23 Mar. 2000 H.Weber : add GetAngularMomentum()
// ------------------------------------------------------------
#ifndef G4VDecayChannel_h
#define G4VDecayChannel_h 1
#include "G4ios.hh"
#include "globals.hh"
#include <cmath>
#include "G4ThreeVector.hh"
class G4ParticleDefinition;
class G4DecayProducts;
class G4ParticleTable;
class G4VDecayChannel
{
// Class Description
// This class is a abstract class to describe decay kinematics
//
public:
//Constructors
G4VDecayChannel(const G4String &aName, G4int Verbose = 1);
G4VDecayChannel(const G4String &aName,
const G4String& theParentName,
G4double theBR,
G4int theNumberOfDaughters,
const G4String& theDaughterName1,
const G4String& theDaughterName2 = "",
const G4String& theDaughterName3 = "",
const G4String& theDaughterName4 = "" );
// Destructor
virtual ~G4VDecayChannel();
private:
// copy constructor and assignment operatotr
G4VDecayChannel(const G4VDecayChannel &);
G4VDecayChannel & operator=(const G4VDecayChannel &);
public:
// equality operators
G4int operator==(const G4VDecayChannel &right) const {return (this == &right);}
G4int operator!=(const G4VDecayChannel &right) const {return (this != &right);}
// less-than operator is defined for G4DecayTable
G4int operator<(const G4VDecayChannel &right) const;
public: // With Description
virtual G4DecayProducts* DecayIt(G4double parentMass = -1.0) = 0;
public: // With Description
//get kinematics name
G4String GetKinematicsName() const;
//get branching ratio
G4double GetBR() const;
//get number of daughter particles
G4int GetNumberOfDaughters() const;
//get the pointer to the parent particle
G4ParticleDefinition * GetParent();
//get the pointer to a daughter particle
G4ParticleDefinition * GetDaughter(G4int anIndex);
//get the angular momentum of the decay
G4int GetAngularMomentum();
//get the name of the parent particle
const G4String& GetParentName() const;
//get the name of a daughter particle
const G4String& GetDaughterName(G4int anIndex) const;
// get mass of parent
G4double GetParentMass() const;
G4double GetDaughterMass(G4int anIndex) const;
G4ThreeVector GetParentPolarization();
//set the parent particle (by name or by pointer)
void SetParent(const G4ParticleDefinition * particle_type);
void SetParent(const G4String &particle_name);
//set branching ratio
void SetBR(G4double value);
//set number of daughter particles
void SetNumberOfDaughters(G4int value);
//set a daughter particle (by name or by pointer)
void SetDaughter(G4int anIndex,
const G4ParticleDefinition * particle_type);
void SetDaughter(G4int anIndex,
const G4String &particle_name);
void SetParentPolarization(G4ThreeVector polar);
protected:
// kinematics name
G4String kinematics_name;
// branching ratio [0.0 - 1.0]
G4double rbranch;
// number of daughters
G4int numberOfDaughters;
// parent particle
G4String* parent_name;
//daughter particles
G4String** daughters_name;
protected: // With Description
// celar daughters array
void ClearDaughtersName();
protected:
// pointer to particle table
G4ParticleTable* particletable;
// temporary buffers of pointers to G4ParticleDefinition
G4ParticleDefinition* parent;
G4ParticleDefinition** daughters;
// parent mass
G4double parent_mass;
G4double* daughters_mass;
//parent theParentPolarization
G4ThreeVector theParentPolarization;
// fill daughters array
void FillDaughters();
// fill parent
void FillParent();
public: // With Description
void SetVerboseLevel(G4int value);
G4int GetVerboseLevel() const;
void DumpInfo();
private:
const G4String& GetNoName() const;
private:
// controle flag for output message
G4int verboseLevel;
// 0: Silent
// 1: Warning message
// 2: More
static const G4String noName;
};
inline
G4int G4VDecayChannel::operator<(const G4VDecayChannel &right) const
{
return (this->rbranch < right.rbranch);
}
inline
G4ParticleDefinition* G4VDecayChannel::GetDaughter(G4int anIndex)
{
//pointers to daughter particles are filled, if they are not set yet
if (daughters == 0) FillDaughters();
//get the pointer to a daughter particle
if ( (anIndex>=0) && (anIndex<numberOfDaughters) ) {
return daughters[anIndex];
} else {
if (verboseLevel>0)
G4cout << "G4VDecayChannel::GetDaughter index out of range "<<anIndex<<G4endl;
return 0;
}
}
inline
const G4String& G4VDecayChannel::GetDaughterName(G4int anIndex) const
{
if ( (anIndex>=0) && (anIndex<numberOfDaughters) ) {
return *daughters_name[anIndex];
} else {
if (verboseLevel>0){
G4cout << "G4VDecayChannel::GetDaughterName ";
G4cout << "index out of range " << anIndex << G4endl;
}
return GetNoName();
}
}
inline
G4double G4VDecayChannel::GetDaughterMass(G4int anIndex) const
{
if ( (anIndex>=0) && (anIndex<numberOfDaughters) ) {
return daughters_mass[anIndex];
} else {
if (verboseLevel>0){
G4cout << "G4VDecayChannel::GetDaughterMass ";
G4cout << "index out of range " << anIndex << G4endl;
}
return 0.0;
}
}
inline
G4ParticleDefinition* G4VDecayChannel::GetParent()
{
//the pointer to the parent particle is filled, if it is not set yet
if (parent == 0) FillParent();
//get the pointer to the parent particle
return parent;
}
inline
const G4String& G4VDecayChannel::GetParentName() const
{
return *parent_name;
}
inline
G4double G4VDecayChannel::GetParentMass() const
{
return parent_mass;
}
inline
void G4VDecayChannel::SetParent(const G4String &particle_name)
{
if (parent_name != 0) delete parent_name;
parent_name = new G4String(particle_name);
parent = 0;
}
inline
G4int G4VDecayChannel::GetNumberOfDaughters() const
{
return numberOfDaughters;
}
inline
G4String G4VDecayChannel::GetKinematicsName() const { return kinematics_name; }
inline
void G4VDecayChannel::SetBR(G4double value){ rbranch = value; }
inline
G4double G4VDecayChannel::GetBR() const { return rbranch; }
inline
void G4VDecayChannel::SetVerboseLevel(G4int value){ verboseLevel = value; }
inline
G4int G4VDecayChannel::GetVerboseLevel() const { return verboseLevel; }
inline G4ThreeVector G4VDecayChannel::GetParentPolarization(){return theParentPolarization;}
inline void G4VDecayChannel::SetParentPolarization(G4ThreeVector polar){theParentPolarization=polar;}
#endif

11
geant4/LEMuSR/G4Modified/Readme.txt Normal file
View File

@ -0,0 +1,11 @@
Geant4 ElectricField update files:
Copy files in the following directories:
G4El_EqRhs -----> $G4Install/source/geometry/magneticfield /src and /include
G4El_UsualEqRhs -----> $G4Install/source/geometry/magneticfield /src and /include
G4ChordFinder -----> $G4Install/source/geometry/magneticfield /src and /include
G4ParticleGun -----> $G4Install/source/event /src and /include
track/ is the geant4.6.2 version of $G4Install/source/track directory and should be updated

57
geant4/LEMuSR/G4Modified/update.sh Normal file
View File

@ -0,0 +1,57 @@
cp G4ChordFinder.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4ChordFinder.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4El_EqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4El_EqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4El_MagEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4El_MagEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4El_UsualEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4El_UsualEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4FieldManager.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4FieldManager.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4FieldManager.icc ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4MultipleScattering52.cc ~/geant4.8.0.p01/source/processes/electromagnetic/standard/src/
cp G4VDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4VDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4MuonDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4MuonDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4MuonDecayChannelWithSpin.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4MuonDecayChannelWithSpin.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4Muonium.cc ~/geant4.8.0.p01/source/particles/leptons/src/
cp G4Muonium.hh ~/geant4.8.0.p01/source/particles/leptons/include/
cp G4ParticleDefinition.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4ParticleDefinition.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4ParticleDefinition.icc ~/geant4.8.0.p01/source/particles/management/include/
cd ~/geant4.8.0.p01/source/geometry/magneticfield/
#gmake clean
gmake
cd ~/geant4.8.0.p01/source/processes/electromagnetic/standard/
#gmake clean
gmake
cd ~/geant4.8.0.p01/source/particles/management/
#gmake clean
gmake
cd ~/geant4.8.0.p01/source/particles/leptons/
#gmake clean
gmake
cd ~/LEMuSR/G4Modified

47
geant4/LEMuSR/G4Modified/update.sh~ Normal file
View File

@ -0,0 +1,47 @@
cp G4ChordFinder.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4ChordFinder.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4El_EqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4El_EqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4El_MagEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4El_MagEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4El_UsualEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4El_UsualEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4FieldManager.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/
cp G4FieldManager.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/
cp G4MultipleScattering52.cc ~/geant4.8.0.p01/source/processes/electromagnetic/standard/src/
cp G4VDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4VDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4MuonDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4MuonDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4MuonDecayChannelWithSpin.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4MuonDecayChannelWithSpin.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4Muonium.cc ~/geant4.8.0.p01/source/particles/leptons/src/
cp G4Muonium.hh ~/geant4.8.0.p01/source/particles/leptons/include/
cp G4ParticleDefinition.cc ~/geant4.8.0.p01/source/particles/management/src/
cp G4ParticleDefinition.hh ~/geant4.8.0.p01/source/particles/management/include/
cp G4ParticleDefinition.icc ~/geant4.8.0.p01/source/particles/management/include/
cd ~/geant4.8.0.p01/source/geometry/magneticfield/
#gmake clean
gmake
cd ~/geant4.8.0.p01/source/processes/electromagnetic/standard/
#gmake clean
gmake
cd ~/geant4.8.0.p01/source/particles/management/
#gmake clean
gmake
cd ~/geant4.8.0.p01/source/particles/leptons/
#gmake clean
gmake
cd ~/LEMuSR/G4Modified