Moved Tao's code to TaoLEMuSR.
This commit is contained in:
shiroka
629
geant4/TaoLEMuSR/G4Modified/G4ChordFinder.cc
Normal file
629
geant4/TaoLEMuSR/G4Modified/G4ChordFinder.cc
Normal 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/TaoLEMuSR/G4Modified/G4ChordFinder.hh
Normal file
198
geant4/TaoLEMuSR/G4Modified/G4ChordFinder.hh
Normal 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/TaoLEMuSR/G4Modified/G4El_EqRhs.cc
Normal file
80
geant4/TaoLEMuSR/G4Modified/G4El_EqRhs.cc
Normal 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/TaoLEMuSR/G4Modified/G4El_EqRhs.hh
Normal file
87
geant4/TaoLEMuSR/G4Modified/G4El_EqRhs.hh
Normal 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/TaoLEMuSR/G4Modified/G4El_MagEqRhs.cc
Normal file
76
geant4/TaoLEMuSR/G4Modified/G4El_MagEqRhs.cc
Normal 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/TaoLEMuSR/G4Modified/G4El_MagEqRhs.hh
Normal file
42
geant4/TaoLEMuSR/G4Modified/G4El_MagEqRhs.hh
Normal 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/TaoLEMuSR/G4Modified/G4El_UsualEqRhs.cc
Normal file
110
geant4/TaoLEMuSR/G4Modified/G4El_UsualEqRhs.cc
Normal 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/TaoLEMuSR/G4Modified/G4El_UsualEqRhs.hh
Normal file
76
geant4/TaoLEMuSR/G4Modified/G4El_UsualEqRhs.hh
Normal 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/TaoLEMuSR/G4Modified/G4FieldManager.cc
Normal file
119
geant4/TaoLEMuSR/G4Modified/G4FieldManager.cc
Normal 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/TaoLEMuSR/G4Modified/G4FieldManager.hh
Normal file
191
geant4/TaoLEMuSR/G4Modified/G4FieldManager.hh
Normal 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/TaoLEMuSR/G4Modified/G4FieldManager.icc
Normal file
133
geant4/TaoLEMuSR/G4Modified/G4FieldManager.icc
Normal 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;
|
||||
}
|
||||
}
|
||||
209
geant4/TaoLEMuSR/G4Modified/G4MuonDecayChannel.cc
Normal file
209
geant4/TaoLEMuSR/G4Modified/G4MuonDecayChannel.cc
Normal file
@ -0,0 +1,209 @@
|
||||
//
|
||||
// ********************************************************************
|
||||
// * 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.13 2005/06/23 11:02:26 gcosmo Exp $
|
||||
// GEANT4 tag $Name: geant4-07-01 $
|
||||
//
|
||||
//
|
||||
// ------------------------------------------------------------
|
||||
// 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
|
||||
//2005
|
||||
// M. Melissas ( melissas AT cppm.in2p3.fr)
|
||||
// J. Brunner ( brunner AT cppm.in2p3.fr)
|
||||
// Adding V-A fluxes for neutrinos using a new algortithm :
|
||||
// ------------------------------------------------------------
|
||||
|
||||
#include "G4ParticleDefinition.hh"
|
||||
#include "G4DecayProducts.hh"
|
||||
#include "G4VDecayChannel.hh"
|
||||
#include "G4MuonDecayChannel.hh"
|
||||
#include "Randomize.hh"
|
||||
#include "G4LorentzVector.hh"
|
||||
#include "G4LorentzRotation.hh"
|
||||
#include "G4RotationMatrix.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,and electron mass
|
||||
// assumes the pure V-A coupling
|
||||
// the Neutrinos are correctly V-A.
|
||||
#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];
|
||||
// calcurate electron energy
|
||||
G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
|
||||
G4double x;
|
||||
|
||||
G4double Ee,Ene;
|
||||
|
||||
G4double gam;
|
||||
G4double EMax=parentmass/2-daughtermass[0];
|
||||
|
||||
|
||||
//Generating Random Energy
|
||||
do {
|
||||
Ee=G4UniformRand();
|
||||
do{
|
||||
x=xmax*G4UniformRand();
|
||||
gam=G4UniformRand();
|
||||
}while (gam >x*(1.-x));
|
||||
Ene=x;
|
||||
} while ( Ene < (1.-Ee));
|
||||
G4double Enm=(2.-Ee-Ene);
|
||||
|
||||
|
||||
//initialisation of rotation parameters
|
||||
|
||||
G4double costheta,sintheta,rphi,rtheta,rpsi;
|
||||
costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
|
||||
sintheta=sqrt(1.-costheta*costheta);
|
||||
|
||||
|
||||
rphi=twopi*G4UniformRand()*rad;
|
||||
rtheta=(acos(2.*G4UniformRand()-1.));
|
||||
rpsi=twopi*G4UniformRand()*rad;
|
||||
|
||||
G4RotationMatrix *rot= new G4RotationMatrix();
|
||||
rot->set(rphi,rtheta,rpsi);
|
||||
|
||||
//electron 0
|
||||
daughtermomentum[0]=sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
|
||||
G4ThreeVector *direction0 =new G4ThreeVector(0.0,0.0,1.0);
|
||||
|
||||
*direction0 *= *rot;
|
||||
|
||||
G4DynamicParticle * daughterparticle = new G4DynamicParticle ( daughters[0], *direction0 * daughtermomentum[0]);
|
||||
|
||||
products->PushProducts(daughterparticle);
|
||||
|
||||
//electronic neutrino 1
|
||||
|
||||
daughtermomentum[1]=sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]);
|
||||
G4ThreeVector *direction1 =new G4ThreeVector(sintheta,0.0,costheta);
|
||||
|
||||
*direction1 *= *rot;
|
||||
|
||||
G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( daughters[1], *direction1 * daughtermomentum[1]);
|
||||
products->PushProducts(daughterparticle1);
|
||||
|
||||
//muonnic neutrino 2
|
||||
|
||||
daughtermomentum[2]=sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]);
|
||||
G4ThreeVector *direction2 =new G4ThreeVector(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
|
||||
|
||||
*direction2 *= *rot;
|
||||
|
||||
G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( daughters[2],
|
||||
*direction2 * daughtermomentum[2]);
|
||||
products->PushProducts(daughterparticle2);
|
||||
|
||||
|
||||
|
||||
|
||||
// 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/TaoLEMuSR/G4Modified/G4MuonDecayChannel.hh
Normal file
63
geant4/TaoLEMuSR/G4Modified/G4MuonDecayChannel.hh
Normal 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-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
|
||||
460
geant4/TaoLEMuSR/G4Modified/G4VDecayChannel.cc
Normal file
460
geant4/TaoLEMuSR/G4Modified/G4VDecayChannel.cc
Normal 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;
|
||||
}
|
||||
290
geant4/TaoLEMuSR/G4Modified/G4VDecayChannel.hh
Normal file
290
geant4/TaoLEMuSR/G4Modified/G4VDecayChannel.hh
Normal 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: 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);
|
||||
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;
|
||||
|
||||
public:
|
||||
void SetParentPolarization(G4ThreeVector polar);
|
||||
};
|
||||
|
||||
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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
8
geant4/TaoLEMuSR/G4Modified/Readme.txt
Normal file
8
geant4/TaoLEMuSR/G4Modified/Readme.txt
Normal file
@ -0,0 +1,8 @@
|
||||
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
|
||||
|
||||
42
geant4/TaoLEMuSR/G4Modified/update.sh
Normal file
42
geant4/TaoLEMuSR/G4Modified/update.sh
Normal file
@ -0,0 +1,42 @@
|
||||
|
||||
|
||||
# In order to handle the LEMuSR new types of field
|
||||
cp G4ChordFinder.cc $G4INSTALL/source/geometry/magneticfield/src/
|
||||
cp G4ChordFinder.hh $G4INSTALL/source/geometry/magneticfield/include/
|
||||
|
||||
# For the FieldHasMagComponent boolean variable
|
||||
cp G4FieldManager.cc $G4INSTALL/source/geometry/magneticfield/src/
|
||||
cp G4FieldManager.hh $G4INSTALL/source/geometry/magneticfield/include/
|
||||
cp G4FieldManager.icc $G4INSTALL/source/geometry/magneticfield/include/
|
||||
|
||||
# Better equations of motions including time update
|
||||
cp G4El_EqRhs.cc $G4INSTALL/source/geometry/magneticfield/src/
|
||||
cp G4El_EqRhs.hh $G4INSTALL/source/geometry/magneticfield/include/
|
||||
|
||||
cp G4El_UsualEqRhs.cc $G4INSTALL/source/geometry/magneticfield/src/
|
||||
cp G4El_UsualEqRhs.hh $G4INSTALL/source/geometry/magneticfield/include/
|
||||
|
||||
# Enable the muon decay channel for "Mu" particle
|
||||
# Set parent polarization
|
||||
cp G4VDecayChannel.cc $G4INSTALL/source/particles/management/src/
|
||||
cp G4VDecayChannel.hh $G4INSTALL/source/particles/management/include/
|
||||
cp G4MuonDecayChannel.cc $G4INSTALL/source/particles/management/src/
|
||||
cp G4MuonDecayChannel.hh $G4INSTALL/source/particles/management/include/
|
||||
|
||||
|
||||
#cd $G4INSTALL/source/geometry/magneticfield/
|
||||
#gmake clean
|
||||
#gmake
|
||||
|
||||
|
||||
#cd $G4INSTALL/source/particles/management/
|
||||
#gmake clean
|
||||
#gmake
|
||||
|
||||
cd $G4INSTALL/source
|
||||
gmake
|
||||
|
||||
cd $G4WORKDIR
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user