From 10aa47ae02f49e14083dce52b64e04b4799ac3a0 Mon Sep 17 00:00:00 2001 From: paraiso Date: Thu, 23 Feb 2006 14:12:45 +0000 Subject: [PATCH] Added to SVN Repository --- geant4/LEMuSR/G4Modified/G4ChordFinder.cc | 629 ++++++++++++ geant4/LEMuSR/G4Modified/G4ChordFinder.hh | 198 ++++ geant4/LEMuSR/G4Modified/G4El_EqRhs.cc | 80 ++ geant4/LEMuSR/G4Modified/G4El_EqRhs.hh | 87 ++ geant4/LEMuSR/G4Modified/G4El_MagEqRhs.cc | 76 ++ geant4/LEMuSR/G4Modified/G4El_MagEqRhs.hh | 42 + geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.cc | 110 ++ geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.hh | 76 ++ geant4/LEMuSR/G4Modified/G4FieldManager.cc | 119 +++ geant4/LEMuSR/G4Modified/G4FieldManager.hh | 191 ++++ geant4/LEMuSR/G4Modified/G4FieldManager.icc | 133 +++ .../G4Modified/G4MultipleScattering52.cc | 951 ++++++++++++++++++ .../LEMuSR/G4Modified/G4MuonDecayChannel.cc | 190 ++++ .../LEMuSR/G4Modified/G4MuonDecayChannel.hh | 63 ++ .../G4Modified/G4MuonDecayChannelWithSpin.cc | 264 +++++ .../G4Modified/G4MuonDecayChannelWithSpin.hh | 122 +++ geant4/LEMuSR/G4Modified/G4Muonium.cc | 119 +++ geant4/LEMuSR/G4Modified/G4Muonium.cc~ | 119 +++ geant4/LEMuSR/G4Modified/G4Muonium.hh | 90 ++ geant4/LEMuSR/G4Modified/G4Muonium.hh~ | 90 ++ .../LEMuSR/G4Modified/G4ParticleDefinition.cc | 290 ++++++ .../LEMuSR/G4Modified/G4ParticleDefinition.hh | 334 ++++++ .../G4Modified/G4ParticleDefinition.icc | 122 +++ geant4/LEMuSR/G4Modified/G4VDecayChannel.cc | 460 +++++++++ geant4/LEMuSR/G4Modified/G4VDecayChannel.hh | 288 ++++++ geant4/LEMuSR/G4Modified/Readme.txt | 11 + geant4/LEMuSR/G4Modified/update.sh | 57 ++ geant4/LEMuSR/G4Modified/update.sh~ | 47 + 28 files changed, 5358 insertions(+) create mode 100644 geant4/LEMuSR/G4Modified/G4ChordFinder.cc create mode 100644 geant4/LEMuSR/G4Modified/G4ChordFinder.hh create mode 100644 geant4/LEMuSR/G4Modified/G4El_EqRhs.cc create mode 100644 geant4/LEMuSR/G4Modified/G4El_EqRhs.hh create mode 100644 geant4/LEMuSR/G4Modified/G4El_MagEqRhs.cc create mode 100644 geant4/LEMuSR/G4Modified/G4El_MagEqRhs.hh create mode 100644 geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.cc create mode 100644 geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.hh create mode 100644 geant4/LEMuSR/G4Modified/G4FieldManager.cc create mode 100644 geant4/LEMuSR/G4Modified/G4FieldManager.hh create mode 100644 geant4/LEMuSR/G4Modified/G4FieldManager.icc create mode 100644 geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc create mode 100644 geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc create mode 100644 geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh create mode 100644 geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc create mode 100644 geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.hh create mode 100644 geant4/LEMuSR/G4Modified/G4Muonium.cc create mode 100644 geant4/LEMuSR/G4Modified/G4Muonium.cc~ create mode 100644 geant4/LEMuSR/G4Modified/G4Muonium.hh create mode 100644 geant4/LEMuSR/G4Modified/G4Muonium.hh~ create mode 100644 geant4/LEMuSR/G4Modified/G4ParticleDefinition.cc create mode 100644 geant4/LEMuSR/G4Modified/G4ParticleDefinition.hh create mode 100644 geant4/LEMuSR/G4Modified/G4ParticleDefinition.icc create mode 100644 geant4/LEMuSR/G4Modified/G4VDecayChannel.cc create mode 100644 geant4/LEMuSR/G4Modified/G4VDecayChannel.hh create mode 100644 geant4/LEMuSR/G4Modified/Readme.txt create mode 100644 geant4/LEMuSR/G4Modified/update.sh create mode 100644 geant4/LEMuSR/G4Modified/update.sh~ diff --git a/geant4/LEMuSR/G4Modified/G4ChordFinder.cc b/geant4/LEMuSR/G4Modified/G4ChordFinder.cc new file mode 100644 index 0000000..e854cdc --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4ChordFinder.cc @@ -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 +#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 < 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" + <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="< 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 < 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); +} diff --git a/geant4/LEMuSR/G4Modified/G4ChordFinder.hh b/geant4/LEMuSR/G4Modified/G4ChordFinder.hh new file mode 100644 index 0000000..4788d65 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4ChordFinder.hh @@ -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 diff --git a/geant4/LEMuSR/G4Modified/G4El_EqRhs.cc b/geant4/LEMuSR/G4Modified/G4El_EqRhs.cc new file mode 100644 index 0000000..4f31978 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4El_EqRhs.cc @@ -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!"<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 +{} diff --git a/geant4/LEMuSR/G4Modified/G4El_MagEqRhs.hh b/geant4/LEMuSR/G4Modified/G4El_MagEqRhs.hh new file mode 100644 index 0000000..dff334f --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4El_MagEqRhs.hh @@ -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 */ diff --git a/geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.cc b/geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.cc new file mode 100644 index 0000000..b66ee6f --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4El_UsualEqRhs.cc @@ -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" + <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; +} diff --git a/geant4/LEMuSR/G4Modified/G4FieldManager.hh b/geant4/LEMuSR/G4Modified/G4FieldManager.hh new file mode 100644 index 0000000..cf39d17 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4FieldManager.hh @@ -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 */ diff --git a/geant4/LEMuSR/G4Modified/G4FieldManager.icc b/geant4/LEMuSR/G4Modified/G4FieldManager.icc new file mode 100644 index 0000000..8866475 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4FieldManager.icc @@ -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; + } +} diff --git a/geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc b/geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc new file mode 100644 index 0000000..e7524a9 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc @@ -0,0 +1,951 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4MultipleScattering52.cc,v 1.2 2004/12/01 19:37:14 vnivanch Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-03 $ +// +// ----------------------------------------------------------------------------- +// 16/05/01 value of cparm changed , L.Urban +// 18/05/01 V.Ivanchenko Clean up against Linux ANSI compilation +// 07/08/01 new methods Store/Retrieve PhysicsTable (mma) +// 23-08-01 new angle and z distribution,energy dependence reduced, +// Store,Retrieve methods commented out temporarily, L.Urban +// 27-08-01 in BuildPhysicsTable:aParticleType.GetParticleName()=="mu+" (mma) +// 28-08-01 GetContinuousStepLimit and AlongStepDoIt moved from .icc file (mma) +// 03-09-01 value of data member factlim changed, L.Urban +// 10-09-01 small change in GetContinuousStepLimit, L.Urban +// 11-09-01 G4MultipleScatteringx put as default G4MultipleScattering +// store/retrieve physics table reactivated (mma) +// 13-09-01 corr. in ComputeTransportCrossSection, L.Urban +// 14-09-01 protection in GetContinuousStepLimit, L.Urban +// 17-09-01 migration of Materials to pure STL (mma) +// 27-09-01 value of data member factlim changed, L.Urban +// 31-10-01 big fixed in PostStepDoIt,L.Urban +// 24-04-02 some minor changes in boundary algorithm, L.Urban +// 06-05-02 bug fixed in GetContinuousStepLimit, L.Urban +// 24-05-02 changes in angle distribution and boundary algorithm, L.Urban +// 11-06-02 bug fixed in ComputeTransportCrossSection, L.Urban +// 12-08-02 bug fixed in PostStepDoIt (lateral displacement), L.Urban +// 15-08-02 new angle distribution, L.Urban +// 26-09-02 angle distribution + boundary algorithm modified, L.Urban +// 15-10-02 temporary fix for proton scattering +// 30-10-02 modified angle distribution,mods in boundary algorithm, +// changes in data members, L.Urban +// 30-10-02 rename variable cm - Ecm, V.Ivanchenko +// 11-12-02 precision problem in ComputeTransportCrossSection +// for small Tkin/for heavy particles cured, L.Urban +// 05-02-03 changes in data members, new sampling for geom. +// path length, step dependence reduced with new +// method +// 17-03-03 cut per region, V.Ivanchenko +// 13-04-03 add initialisation in GetContinuesStepLimit +// + change table size (V.Ivanchenko) +// 26-04-03 fix problems of retrieve tables (M.Asai) +// 23-05-03 important change in angle distribution for muons/hadrons +// the central part now is similar to the Highland parametrization + +// minor correction in angle sampling algorithm (for all particles) +// (L.Urban) +// 24-05-03 bug in nuclear size corr.computation fixed thanks to Vladimir(L.Urban) +// 30-05-03 misprint in PostStepDoIt corrected(L.Urban) +// 08-08-03 This class is frozen at the release 5.2 (V.Ivanchenko) +// 08-11-04 Remove Store/Retrieve tables (V.Ivantchenko) +// ----------------------------------------------------------------------------- +// +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4MultipleScattering52.hh" +#include "G4StepStatus.hh" +#include "G4Navigator.hh" +#include "G4TransportationManager.hh" +#include "Randomize.hh" +#include "G4ProductionCutsTable.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +using namespace std; + +G4MultipleScattering52::G4MultipleScattering52(const G4String& processName) + : G4VContinuousDiscreteProcess(processName), + theTransportMeanFreePathTable(0), + taubig(8.0),tausmall(1.e-14),taulim(1.e-5), + LowestKineticEnergy(0.1*keV), + HighestKineticEnergy(100.*TeV), + TotBin(100), + materialIndex(0), + tLast (0.0), + zLast (0.0), + boundary(true), + facrange(0.199),tlimit(1.e10*mm),tlimitmin(1.e-7*mm), + cf(1.001), + stepno(0),stepnolastmsc(-1000000),nsmallstep(5), + laststep(0.), + valueGPILSelectionMSC(NotCandidateForSelection), + zmean(0.),samplez(true), + range(1.),T0(1.),T1(1.),lambda0(1.),lambda1(-1.), + Tlow(0.),alam(1.),blam(1.),dtrl(0.15), + lambdam(-1.),clam(1.),zm(1.),cthm(1.), + fLatDisplFlag(true), + NuclCorrPar (0.0615), + FactPar(0.40), + facxsi(1.) + { } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4MultipleScattering52::~G4MultipleScattering52() +{ + if(theTransportMeanFreePathTable) + { + theTransportMeanFreePathTable->clearAndDestroy(); + delete theTransportMeanFreePathTable; + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void G4MultipleScattering52::BuildPhysicsTable( + const G4ParticleDefinition& aParticleType) +{ + // set values of some data members + if((aParticleType.GetParticleName() == "e-") || + (aParticleType.GetParticleName() == "e+")) + { + // parameters for e+/e- + alfa1 = 1.45 ; + alfa2 = 0.60 ; + alfa3 = 0.30 ; + b = 1. ; + xsi = facxsi*2.22 ; + c0 = 2.30 ; + } + else + { + // parameters for heavy particles + alfa1 = 1.10 ; + alfa2 = 0.14 ; + alfa3 = 0.07 ; + b = 1. ; + xsi = facxsi*2.70 ; + c0 = 1.40 ; + } + + // .............................. + Tlow = aParticleType.GetPDGMass(); + + // tables are built for MATERIALS + const G4double sigmafactor = twopi*classic_electr_radius* + classic_electr_radius; + G4double KineticEnergy,AtomicNumber,AtomicWeight,sigma,lambda; + G4double density; + + // destroy old tables if any + if (theTransportMeanFreePathTable) + { + theTransportMeanFreePathTable->clearAndDestroy(); + delete theTransportMeanFreePathTable; + } + + // create table + const G4ProductionCutsTable* theCoupleTable= + G4ProductionCutsTable::GetProductionCutsTable(); + size_t numOfCouples = theCoupleTable->GetTableSize(); + + theTransportMeanFreePathTable = new G4PhysicsTable(numOfCouples); + + // loop for materials + for (size_t i=0; i + GetMaterialCutsCouple(i); + const G4Material* material = couple->GetMaterial(); + const G4ElementVector* theElementVector = material->GetElementVector(); + const G4double* NbOfAtomsPerVolume = + material->GetVecNbOfAtomsPerVolume(); + const G4int NumberOfElements = material->GetNumberOfElements(); + density = material->GetDensity(); + + // loop for kinetic energy values + for (G4int i=0; iGetLowEdgeEnergy(i); + sigma = 0.; + + // loop for element in the material + for (G4int iel=0; ielGetZ(); + AtomicWeight = (*theElementVector)[iel]->GetA(); + sigma += NbOfAtomsPerVolume[iel]* + ComputeTransportCrossSection(aParticleType,KineticEnergy, + AtomicNumber,AtomicWeight); + } + sigma *= sigmafactor; + lambda = 1./sigma; + aVector->PutValue(i,lambda); + } + + theTransportMeanFreePathTable->insert(aVector); + + } + + if((aParticleType.GetParticleName() == "e-" ) || + (aParticleType.GetParticleName() == "mu+" ) || + (aParticleType.GetParticleName() == "Mu" ) || + (aParticleType.GetParticleName() == "proton") ) PrintInfoDefinition(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4double G4MultipleScattering52::ComputeTransportCrossSection( + const G4ParticleDefinition& aParticleType, + G4double KineticEnergy, + G4double AtomicNumber,G4double AtomicWeight) +{ + const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2* + Bohr_radius*Bohr_radius/(hbarc*hbarc); + const G4double epsmin = 1.e-4 , epsmax = 1.e10; + + const G4double Zdat[15] = { 4., 6.,13.,20.,26.,29.,32.,38.,47., + 50.,56.,64.,74.,79.,82. }; + + const G4double Tdat[23] = {0.0001*MeV,0.0002*MeV,0.0004*MeV,0.0007*MeV, + 0.001*MeV,0.002*MeV,0.004*MeV,0.007*MeV, + 0.01*MeV,0.02*MeV,0.04*MeV,0.07*MeV, + 0.1*MeV,0.2*MeV,0.4*MeV,0.7*MeV, + 1.*MeV,2.*MeV,4.*MeV,7.*MeV,10.*MeV,20.*MeV, + 10000.0*MeV}; + + // corr. factors for e-/e+ lambda + G4double celectron[15][23] = + {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054, + 1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111, + 1.112,1.108,1.100,1.093,1.089,1.087,0.7235 }, + {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051, + 1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108, + 1.109,1.105,1.097,1.090,1.086,1.082,0.7925 }, + {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156, + 1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132, + 1.131,1.124,1.113,1.104,1.099,1.098,0.9147 }, + {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236, + 1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113, + 1.112,1.105,1.096,1.089,1.085,1.098,0.9700 }, + {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265, + 1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073, + 1.073,1.070,1.064,1.059,1.056,1.056,1.0022 }, + {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330, + 1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074, + 1.074,1.070,1.063,1.059,1.056,1.052,1.0158 }, + {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386, + 1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069, + 1.068,1.064,1.059,1.054,1.051,1.050,1.0284 }, + {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439, + 1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039, + 1.039,1.037,1.034,1.031,1.030,1.036,1.0515 }, + {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631, + 1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033, + 1.031,1.028,1.024,1.022,1.021,1.024,1.0834 }, + {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669, + 1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022, + 1.020,1.017,1.015,1.013,1.013,1.020,1.0937 }, + {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720, + 1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997, + 0.995,0.993,0.993,0.993,0.993,1.011,1.1140 }, + {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855, + 1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976, + 0.974,0.972,0.973,0.974,0.975,0.987,1.1410 }, + {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059, + 1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954, + 0.950,0.947,0.949,0.952,0.954,0.963,1.1750 }, + {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182, + 1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947, + 0.941,0.938,0.940,0.944,0.946,0.954,1.1922 }, + // {45.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239, // paper..... + {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239, + 1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939, + 0.933,0.930,0.933,0.936,0.939,0.949,1.2026 }}; + G4double cpositron[15][23] = { + {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110, + 1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131, + 1.131,1.126,1.117,1.108,1.103,1.100,0.7235 }, + {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145, + 1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137, + 1.138,1.132,1.122,1.113,1.108,1.102,0.7925 }, + {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451, + 1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205, + 1.203,1.190,1.173,1.159,1.151,1.145,0.9147 }, + {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715, + 1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228, + 1.225,1.210,1.191,1.175,1.166,1.174,0.9700 }, + {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820, + 1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219, + 1.217,1.203,1.184,1.169,1.160,1.151,1.0022 }, + {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996, + 1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241, + 1.237,1.222,1.201,1.184,1.174,1.159,1.0158 }, + {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155, + 1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256, + 1.252,1.234,1.212,1.194,1.183,1.170,1.0284 }, + {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348, + 2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258, + 1.254,1.237,1.214,1.195,1.185,1.179,1.0515 }, + {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808, + 2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320, + 1.312,1.288,1.258,1.235,1.221,1.205,1.0834 }, + {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917, + 2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327, + 1.320,1.294,1.264,1.240,1.226,1.214,1.0937 }, + {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066, + 2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336, + 1.328,1.302,1.270,1.245,1.231,1.233,1.1140 }, + {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498, + 2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371, + 1.361,1.330,1.294,1.267,1.251,1.239,1.1410 }, + {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155, + 3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423, + 1.409,1.372,1.330,1.298,1.280,1.258,1.1750 }, + {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407, + 3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459, + 1.442,1.400,1.354,1.319,1.299,1.272,1.1922 }, + {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542, + 3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474, + 1.456,1.412,1.364,1.328,1.307,1.282,1.2026 }}; + + G4double sigma; + G4double Z23 = 2.*log(AtomicNumber)/3.; Z23 = exp(Z23); + + G4double ParticleMass = aParticleType.GetPDGMass(); + G4double ParticleKineticEnergy = KineticEnergy ; + + // correction if particle .ne. e-/e+ + // compute equivalent kinetic energy + // lambda depends on p*beta .... + G4double Mass = ParticleMass ; + if((aParticleType.GetParticleName() != "e-") && + (aParticleType.GetParticleName() != "e+") ) + { + G4double TAU = KineticEnergy/Mass ; + G4double c = Mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ; + G4double w = c-2. ; + G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ; + KineticEnergy = electron_mass_c2*tau ; + Mass = electron_mass_c2 ; + } + + G4double Charge = aParticleType.GetPDGCharge(); + G4double ChargeSquare = Charge*Charge/(eplus*eplus); + + G4double TotalEnergy = KineticEnergy + Mass ; + G4double beta2 = KineticEnergy*(TotalEnergy+Mass) + /(TotalEnergy*TotalEnergy); + G4double bg2 = KineticEnergy*(TotalEnergy+Mass) + /(Mass*Mass); + + G4double eps = epsfactor*bg2/Z23; + + if (eps epsmax) w1=log(2.*eps)+1./eps-3./(8.*eps*eps); + else w1=log((a+1.)/(a-1.))-2./(a+1.); + w = 1./((1.-x0)*eps); + if (w < epsmin) w2=-log(w)-1.+2.*w-1.5*w*w; + else w2 = log((a-x0)/(a-1.))-(1.-x0)/(a-x0); + corrnuclsize = w1/w2; + corrnuclsize = exp(-FactPar*ParticleMass/ParticleKineticEnergy)* + (corrnuclsize-1.)+1.; + } + + // interpolate in AtomicNumber and beta2 + // get bin number in Z + G4int iZ = 14; + while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1; + if (iZ==14) iZ = 13; + if (iZ==-1) iZ = 0 ; + + G4double Z1 = Zdat[iZ]; + G4double Z2 = Zdat[iZ+1]; + G4double ratZ = (AtomicNumber-Z1)/(Z2-Z1); + + // get bin number in T (beta2) + G4int iT = 22; + while ((iT>=0)&&(Tdat[iT]>=KineticEnergy)) iT -= 1; + if(iT==22) iT = 21; + if(iT==-1) iT = 0 ; + + // calculate betasquare values + G4double T = Tdat[iT], E = T + electron_mass_c2; + G4double b2small = T*(E+electron_mass_c2)/(E*E); + T = Tdat[iT+1]; E = T + electron_mass_c2; + G4double b2big = T*(E+electron_mass_c2)/(E*E); + G4double ratb2 = (beta2-b2small)/(b2big-b2small); + + G4double c1,c2,cc1,cc2,corr; + + if (Charge < 0.) + { + c1 = celectron[iZ][iT]; + c2 = celectron[iZ+1][iT]; + cc1 = c1+ratZ*(c2-c1); + + c1 = celectron[iZ][iT+1]; + c2 = celectron[iZ+1][iT+1]; + cc2 = c1+ratZ*(c2-c1); + + corr = cc1+ratb2*(cc2-cc1); + sigma /= corr; + } + + if (Charge > 0.) + { + c1 = cpositron[iZ][iT]; + c2 = cpositron[iZ+1][iT]; + cc1 = c1+ratZ*(c2-c1); + + c1 = cpositron[iZ][iT+1]; + c2 = cpositron[iZ+1][iT+1]; + cc2 = c1+ratZ*(c2-c1); + + corr = cc1+ratb2*(cc2-cc1); + sigma /= corr; + } + + // nucl. size correction for particles other than e+/e- only at present !!!! + if((aParticleType.GetParticleName() != "e-") && + (aParticleType.GetParticleName() != "e+") ) + sigma /= corrnuclsize; + + return sigma; + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4double G4MultipleScattering52::GetContinuousStepLimit( + const G4Track& track, + G4double, + G4double currentMinimumStep, + G4double&) +{ + G4double zPathLength,tPathLength; + const G4DynamicParticle* aParticle; + G4double tau,zt,cz,cz1,grej,grej0; + const G4double expmax = 100., ztmax = (2.*expmax+1.)/(2.*expmax+3.) ; + const G4double tmax = 1.e20*mm ; + + G4bool isOut; + + // this process is not a candidate for selection by default + valueGPILSelectionMSC = NotCandidateForSelection; + + tPathLength = currentMinimumStep; + + const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); + materialIndex = couple->GetIndex(); + + aParticle = track.GetDynamicParticle(); + T0 = aParticle->GetKineticEnergy(); + + lambda0 = (*theTransportMeanFreePathTable) + (materialIndex)->GetValue(T0,isOut); + + + range = G4EnergyLossTables::GetRange(aParticle->GetDefinition(), + T0,couple); + //VI Initialisation at the beginning of the step + cthm = 1.; + lambda1 = -1.; + lambdam = -1.; + alam = range; + blam = 1.+alam/lambda0 ; + zm = 1.; + + // special treatment near boundaries ? + if (boundary && range >= currentMinimumStep) + { + // step limitation at boundary ? + stepno = track.GetCurrentStepNumber() ; + if(stepno == 1) + { + stepnolastmsc = -1000000 ; + tlimit = 1.e10 ; + } + + if(stepno > 1) + { + if(track.GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) + { + stepnolastmsc = stepno ; + // if : diff.treatment for small/not small Z + if(range > lambda0) + tlimit = facrange*range ; + else + tlimit = facrange*lambda0 ; + if(tlimit < tlimitmin) tlimit = tlimitmin ; + laststep = tlimit ; + if(tPathLength > tlimit) + { + tPathLength = tlimit ; + valueGPILSelectionMSC = CandidateForSelection; + } + } + else if(stepno > stepnolastmsc) + { + if((stepno - stepnolastmsc) < nsmallstep) + { + if(tPathLength > tlimit) + { + laststep *= cf ; + tPathLength = laststep ; + valueGPILSelectionMSC = CandidateForSelection; + } + } + } + } + } + + // do the true -> geom transformation + zmean = tPathLength; + + tau = tPathLength/lambda0 ; + + if (tau < tausmall || range < currentMinimumStep) zPathLength = tPathLength; + else + { + if(tPathLength/range < dtrl) zmean = lambda0*(1.-exp(-tau)); + else + { + T1 = G4EnergyLossTables::GetPreciseEnergyFromRange( + aParticle->GetDefinition(),range-tPathLength,couple); + lambda1 = (*theTransportMeanFreePathTable) + (materialIndex)->GetValue(T1,isOut); + if(T0 < Tlow) + alam = range ; + else + alam = lambda0*tPathLength/(lambda0-lambda1) ; + blam = 1.+alam/lambda0 ; + if(tPathLength/range < 2.*dtrl) + { + zmean = alam*(1.-exp(blam*log(1.-tPathLength/alam)))/blam ; + lambdam = -1. ; + } + else + { + G4double w = 1.-0.5*tPathLength/alam ; + lambdam = lambda0*w ; + clam = 1.+alam/lambdam ; + cthm = exp(alam*log(w)/lambda0) ; + zm = alam*(1.-exp(blam*log(w)))/blam ; + zmean = zm + alam*(1.-exp(clam*log(w)))*cthm/clam ; + } + } + + // sample z + zt = zmean/tPathLength ; + if (samplez && (zt < ztmax) && (zt > 0.5)) + { + cz = 0.5*(3.*zt-1.)/(1.-zt) ; + if(tPathLength < exp(log(tmax)/(2.*cz))) + { + cz1 = 1.+cz ; + grej0 = exp(cz1*log(cz*tPathLength/cz1))/cz ; + do + { + zPathLength = tPathLength*exp(log(G4UniformRand())/cz1) ; + grej = exp(cz*log(zPathLength))*(tPathLength-zPathLength)/grej0 ; + } while (grej < G4UniformRand()) ; + } + else zPathLength = zmean; + } + else zPathLength = zmean; + } + // protection against z > lambda + if(zPathLength > lambda0) + zPathLength = lambda0 ; + + tLast = tPathLength; + zLast = zPathLength; + + return zPathLength; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4VParticleChange* G4MultipleScattering52::AlongStepDoIt( + const G4Track& track,const G4Step& step) +{ + // only a geom path->true path transformation is performed + + fParticleChange.Initialize(track); + + G4double geomPathLength = step.GetStepLength(); + + G4double truePathLength = 0. ; + + //VI change order of if operators + if(geomPathLength == zLast) truePathLength = tLast; + else if(geomPathLength/lambda0 < tausmall) truePathLength = geomPathLength; + else + { + if(lambda1 < 0.) truePathLength = -lambda0*log(1.-geomPathLength/lambda0) ; + else if(lambdam < 0.) + { + if(blam*geomPathLength/alam < 1.) + truePathLength = alam*(1.-exp(log(1.-blam*geomPathLength/alam)/ + blam)) ; + else + truePathLength = tLast; + } + else + { + if(geomPathLength <= zm) + { + if(blam*geomPathLength/alam < 1.) + truePathLength = alam*(1.-exp(log(1.-blam*geomPathLength/alam)/ + blam)) ; + else + truePathLength = 0.5*tLast; + + lambdam = -1. ; + } + else + { + if(clam*(geomPathLength-zm)/(alam*cthm) < 1.) + truePathLength = 0.5*tLast + alam*(1.- + exp(log(1.-clam*(geomPathLength-zm)/(alam*cthm)))/clam) ; + else + truePathLength = tLast ; + } + } + // protection .... + if(truePathLength > tLast) + truePathLength = tLast ; + } + + //VI truePath length cannot be smaller than geomPathLength + if (truePathLength < geomPathLength) truePathLength = geomPathLength; + fParticleChange.ProposeTrueStepLength(truePathLength); + + return &fParticleChange; + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4VParticleChange* G4MultipleScattering52::PostStepDoIt( + const G4Track& trackData, + const G4Step& stepData) +{ + // angle distribution parameters + const G4double kappa = 2.5, kappapl1 = kappa+1., kappami1 = kappa-1. ; + + fParticleChange.Initialize(trackData); + G4double truestep = stepData.GetStepLength(); + + const G4DynamicParticle* aParticle = trackData.GetDynamicParticle(); + G4double KineticEnergy = aParticle->GetKineticEnergy(); + G4double Mass = aParticle->GetDefinition()->GetPDGMass() ; + + // do nothing for stopped particles ! + if(KineticEnergy > 0.) + { + // change direction first ( scattering ) + G4double cth = 1.0 ; + G4double tau = truestep/lambda0 ; + + if (tau < tausmall) cth = 1.; + else if(tau > taubig) cth = -1.+2.*G4UniformRand(); + else + { + if(lambda1 > 0.) + { + if(lambdam < 0.) + tau = -alam*log(1.-truestep/alam)/lambda0 ; + else + tau = -log(cthm)-alam*log(1.-(truestep-0.5*tLast)/alam)/lambdam ; + } + + if(tau > taubig) cth = -1.+2.*G4UniformRand(); + else + { + const G4double amax=25. ; + const G4double tau0 = 0.02 ; + const G4double c_highland = 13.6*MeV, corr_highland=0.038 ; + + const G4double x1fac1 = exp(-xsi) ; + const G4double x1fac2 = (1.-(1.+xsi)*x1fac1)/(1.-x1fac1) ; + const G4double x1fac3 = 1.3 ; // x1fac3 >= 1. !!!!!!!!! + + G4double a,x0,c,xmean1,xmean2, + xmeanth,prob,qprob ; + G4double ea,eaa,b1,bx,eb1,ebx,cnorm1,cnorm2,f1x0,f2x0,w ; + + // for heavy particles take the width of the cetral part + // from the Highland formula + // (Particle Physics Booklet, July 2002, eq. 26.10) + if(Mass > electron_mass_c2) // + other conditions (beta, x/X0,...?) + { + G4double Q = fabs(aParticle->GetDefinition()->GetPDGCharge()) ; + G4double X0 = trackData.GetMaterialCutsCouple()-> + GetMaterial()->GetRadlen() ; + G4double xx0 = truestep/X0 ; + G4double betacp = KineticEnergy*(KineticEnergy+2.*Mass)/ + (KineticEnergy+Mass) ; + G4double theta0=c_highland*Q*sqrt(xx0)* + (1.+corr_highland*log(xx0))/betacp ; + if(theta0 > tausmall) + a = 0.5/(1.-cos(theta0)) ; + else + a = 1./(theta0*theta0) ; + } + else + { + w = log(tau/tau0) ; + if(tau < tau0) + a = (alfa1-alfa2*w)/tau ; + else + a = (alfa1+alfa3*w)/tau ; + } + + xmeanth = exp(-tau) ; + + x0 = 1.-xsi/a ; + if(x0 < -1.) x0 = -1. ; + + if(x0 == -1.) + { + // 1 model fuction only + // in order to have xmean1 > xmeanth -> qprob < 1 + if((1.-1./a) < xmeanth) + a = 1./(1.-xmeanth) ; + + if(a*(1.-x0) < amax) + ea = exp(-a*(1.-x0)) ; + else + ea = 0. ; + eaa = 1.-ea ; + xmean1 = 1.-1./a+(1.-x0)*ea/eaa ; + + c = 2. ; + b1 = b+1. ; + bx = b1 ; + eb1 = b1 ; + ebx = b1 ; + xmean2 = 0. ; + + prob = 1. ; + qprob = xmeanth/xmean1 ; + } + else + { + // 2 model fuctions + // in order to have xmean1 > xmeanth + if((1.-x1fac2/a) < xmeanth) + { + a = x1fac3*x1fac2/(1.-xmeanth) ; + if(a*(1.-x0) < amax) + ea = exp(-a*(1.-x0)) ; + else + ea = 0. ; + eaa = 1.-ea ; + xmean1 = 1.-1./a+(1.-x0)*ea/eaa ; + } + else + { + ea = x1fac1 ; + eaa = 1.-x1fac1 ; + xmean1 = 1.-x1fac2/a ; + } + + // from continuity of the 1st derivatives + c = a*(b-x0) ; + if(a*tau < c0) + c = c0*(b-x0)/tau ; + + if(c == 1.) c=1.000001 ; + if(c == 2.) c=2.000001 ; + if(c == 3.) c=3.000001 ; + + b1 = b+1. ; + bx=b-x0 ; + eb1=exp((c-1.)*log(b1)) ; + ebx=exp((c-1.)*log(bx)) ; + xmean2 = (x0*eb1+ebx+(eb1*bx-b1*ebx)/(2.-c))/(eb1-ebx) ; + + cnorm1 = a/eaa ; + f1x0 = cnorm1*exp(-a*(1.-x0)) ; + cnorm2 = (c-1.)*eb1*ebx/(eb1-ebx) ; + f2x0 = cnorm2/exp(c*log(b-x0)) ; + + // from continuity at x=x0 + prob = f2x0/(f1x0+f2x0) ; + // from xmean = xmeanth + qprob = (f1x0+f2x0)*xmeanth/(f2x0*xmean1+f1x0*xmean2) ; + } + + // protection against prob or qprob > 1 and + // prob or qprob < 0 + // *************************************************************** + if((qprob > 1.) || (qprob < 0.) || (prob > 1.) || (prob < 0.)) + { + // this print possibility has been left intentionally + // for debugging purposes .......................... + G4bool pr = false ; + // pr = true ; + if(pr) + { + const G4double prlim = 0.10 ; + if((fabs((xmeanth-xmean2)/(xmean1-xmean2)-prob)/prob > prlim) || + ((xmeanth-xmean2)/(xmean1-xmean2) > 1.) || + ((xmeanth-xmean2)/(xmean1-xmean2) < 0.) ) + { + G4cout.precision(5) ; + G4cout << "\nparticle=" << aParticle->GetDefinition()-> + GetParticleName() << " in material " + << trackData.GetMaterialCutsCouple()-> + GetMaterial()->GetName() << " with kinetic energy " + << KineticEnergy << " MeV," << G4endl ; + G4cout << " step length=" + << truestep << " mm" << G4endl ; + G4cout << "p=" << prob << " q=" << qprob << " -----> " + << "p=" << (xmeanth-xmean2)/(xmean1-xmean2) + << " q=" << 1. << G4endl ; + } + } + qprob = 1. ; + prob = (xmeanth-xmean2)/(xmean1-xmean2) ; + } + // ************************************************************** + + // sampling of costheta + if(G4UniformRand() < qprob) + { + if(G4UniformRand() < prob) + cth = 1.+log(ea+G4UniformRand()*eaa)/a ; + else + cth = b-b1*bx/exp(log(ebx-G4UniformRand()*(ebx-eb1))/(c-1.)) ; + } + else + cth = -1.+2.*G4UniformRand() ; + } + } + + G4double sth = sqrt(1.-cth*cth); + G4double phi = twopi*G4UniformRand(); + G4double dirx = sth*cos(phi), diry = sth*sin(phi), dirz = cth; + + G4ParticleMomentum ParticleDirection = aParticle->GetMomentumDirection(); + + G4ThreeVector newDirection(dirx,diry,dirz); + newDirection.rotateUz(ParticleDirection); + fParticleChange.ProposeMomentumDirection(newDirection.x(), + newDirection.y(), + newDirection.z()); + + if (fLatDisplFlag) + { + // compute mean lateral displacement, only for safety > tolerance ! + G4double safetyminustolerance = stepData.GetPostStepPoint()->GetSafety(); + G4double rmean, etau; + + if (safetyminustolerance > 0.) + { + if (tau < tausmall) rmean = 0.; + else if(tau < taulim) rmean = kappa*tau*tau*tau*(1.-kappapl1*tau/4.)/6.; + else + { + if(tau0.) rmean = 2.*lambda0*sqrt(rmean/3.); + else rmean = 0.; + + // for rmean > 0) only + if (rmean > 0.) + { + if (rmean>safetyminustolerance) rmean = safetyminustolerance; + + // sample direction of lateral displacement + phi = twopi*G4UniformRand(); + dirx = cos(phi); diry = sin(phi); dirz = 0.; + + G4ThreeVector latDirection(dirx,diry,dirz); + latDirection.rotateUz(ParticleDirection); + + // compute new endpoint of the Step + G4ThreeVector newPosition = stepData.GetPostStepPoint()->GetPosition() + + rmean*latDirection; + + G4Navigator* navigator = + G4TransportationManager::GetTransportationManager() + ->GetNavigatorForTracking(); + navigator->LocateGlobalPointWithinVolume(newPosition); + + fParticleChange.ProposePosition(newPosition); + } + } + } + } + + return &fParticleChange; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void G4MultipleScattering52::PrintInfoDefinition() +{ + G4String comments = " Tables of transport mean free paths."; + comments += "\n New model of MSC , computes the lateral \n"; + comments += " displacement of the particle , too."; + + G4cout << G4endl << GetProcessName() << ": " << comments + << "\n PhysicsTables from " + << G4BestUnit(LowestKineticEnergy ,"Energy") + << " to " << G4BestUnit(HighestKineticEnergy,"Energy") + << " in " << TotBin << " bins. \n"; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc new file mode 100644 index 0000000..bdb9345 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.cc @@ -0,0 +1,190 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4MuonDecayChannel.cc,v 1.11 2004/12/10 18:02:04 gcosmo Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-05 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 30 May 1997 H.Kurashige +// +// Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige +// ------------------------------------------------------------ + +#include "G4ParticleDefinition.hh" +#include "G4DecayProducts.hh" +#include "G4VDecayChannel.hh" +#include "G4MuonDecayChannel.hh" +#include "Randomize.hh" +#include "G4LorentzVector.hh" +#include "G4LorentzRotation.hh" + + +G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName, + G4double theBR) + :G4VDecayChannel("Muon Decay",1) +{ + // set names for daughter particles + if (theParentName == "mu+") { + SetBR(theBR); + SetParent("mu+"); + SetNumberOfDaughters(3); + SetDaughter(0, "e+"); + SetDaughter(1, "nu_e"); + SetDaughter(2, "anti_nu_mu"); + } else if (theParentName == "Mu") { + SetBR(theBR); + SetParent("Mu"); + SetNumberOfDaughters(3); + SetDaughter(0, "e+"); + SetDaughter(1, "nu_e"); + SetDaughter(2, "anti_nu_mu"); + } else if (theParentName == "mu-") { + SetBR(theBR); + SetParent("mu-"); + SetNumberOfDaughters(3); + SetDaughter(0, "e-"); + SetDaughter(1, "anti_nu_e"); + SetDaughter(2, "nu_mu"); + } else { +#ifdef G4VERBOSE + if (GetVerboseLevel()>0) { + G4cout << "G4MuonDecayChannel:: constructor :"; + G4cout << " parent particle is not muon but "; + G4cout << theParentName << G4endl; + } +#endif + } +} + +G4MuonDecayChannel::~G4MuonDecayChannel() +{ +} + +G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double) +{ + // this version neglects muon polarization + // assumes the pure V-A coupling + // gives incorrect energy spectrum for neutrinos +#ifdef G4VERBOSE + if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt "; +#endif + + if (parent == 0) FillParent(); + if (daughters == 0) FillDaughters(); + + // parent mass + G4double parentmass = parent->GetPDGMass(); + + //daughters'mass + G4double daughtermass[3]; + G4double sumofdaughtermass = 0.0; + for (G4int index=0; index<3; index++){ + daughtermass[index] = daughters[index]->GetPDGMass(); + sumofdaughtermass += daughtermass[index]; + } + + //create parent G4DynamicParticle at rest + G4ThreeVector dummy; + G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); + //create G4Decayproducts + G4DecayProducts *products = new G4DecayProducts(*parentparticle); + delete parentparticle; + + // calculate daughter momentum + G4double daughtermomentum[3]; + G4double energy; + // calcurate electron energy + G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass); + G4double x; + G4double r; + do { + do { + r = G4UniformRand(); + x = xmax*G4UniformRand(); + } while (r > (3.0 - 2.0*x)*x*x); + energy = x*parentmass/2.0 - daughtermass[0]; + } while (energy <0.0); + //create daughter G4DynamicParticle + // daughter 0 (electron) + daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]); + G4double costheta, sintheta, phi, sinphi, cosphi; + costheta = 2.*G4UniformRand()-1.0; + sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); + phi = twopi*G4UniformRand()*rad; + sinphi = std::sin(phi); + cosphi = std::cos(phi); + G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta); + G4DynamicParticle * daughterparticle + = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]); + products->PushProducts(daughterparticle); + + // daughter 1 ,2 (nutrinos) + // create neutrinos in the C.M frame of two neutrinos + G4double energy2 = parentmass*(1.0 - x/2.0); + G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0])); + G4double beta = -1.0*daughtermomentum[0]/energy2; + G4double costhetan = 2.*G4UniformRand()-1.0; + G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); + G4double phin = twopi*G4UniformRand()*rad; + G4double sinphin = std::sin(phin); + G4double cosphin = std::cos(phin); + + G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan); + G4DynamicParticle * daughterparticle1 + = new G4DynamicParticle( daughters[1], direction1*(vmass/2.)); + G4DynamicParticle * daughterparticle2 + = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.)); + + // boost to the muon rest frame + G4LorentzVector p4; + p4 = daughterparticle1->Get4Momentum(); + p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); + daughterparticle1->Set4Momentum(p4); + p4 = daughterparticle2->Get4Momentum(); + p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); + daughterparticle2->Set4Momentum(p4); + products->PushProducts(daughterparticle1); + products->PushProducts(daughterparticle2); + daughtermomentum[1] = daughterparticle1->GetTotalMomentum(); + daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); + + // output message +#ifdef G4VERBOSE + if (GetVerboseLevel()>1) { + G4cout << "G4MuonDecayChannel::DecayIt "; + G4cout << " create decay products in rest frame " <DumpInfo(); + } +#endif + return products; +} + + + + + + diff --git a/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh new file mode 100644 index 0000000..0a3f6d7 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4MuonDecayChannel.hh @@ -0,0 +1,63 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4MuonDecayChannel.hh,v 1.5 2001/07/11 10:01:56 gunter Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-01 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 30 May 1997 H.Kurashige +// ------------------------------------------------------------ +#ifndef G4MuonDecayChannel_h +#define G4MuonDecayChannel_h 1 + +#include "G4ios.hh" +#include "globals.hh" +#include "G4VDecayChannel.hh" + +class G4MuonDecayChannel :public G4VDecayChannel +{ + // Class Decription + // This class describes muon decay kinemtics. + // This version neglects muon polarization + // assumes the pure V-A coupling + // gives incorrect energy spectrum for neutrinos + // + + public: // With Description + //Constructors + G4MuonDecayChannel(const G4String& theParentName, + G4double theBR); + // Destructor + virtual ~G4MuonDecayChannel(); + + public: // With Description + virtual G4DecayProducts *DecayIt(G4double); + +}; + + +#endif diff --git a/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc b/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc new file mode 100644 index 0000000..eb00503 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc @@ -0,0 +1,264 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: +// 17 August 2004 P.Gumplinger and T.MacPhail +// samples Michel spectrum including 1st order +// radiative corrections +// Reference: Florian Scheck "Muon Physics", in Physics Reports +// (Review Section of Physics Letters) 44, No. 4 (1978) +// 187-248. North-Holland Publishing Company, Amsterdam +// at page 210 cc. +// +// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25. +// +// ------------------------------------------------------------ +// +#include "G4MuonDecayChannelWithSpin.hh" + +#include "Randomize.hh" +#include "G4DecayProducts.hh" +#include "G4LorentzVector.hh" + +G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4String& theParentName, + G4double theBR) + : G4MuonDecayChannel(theParentName,theBR) +{ +} + + +G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin() +{ +} + +G4DecayProducts *G4MuonDecayChannelWithSpin::DecayIt(G4double) +{ + // This version assumes V-A coupling with 1st order radiative correctons, + // the standard model Michel parameter values, but + // gives incorrect energy spectrum for neutrinos + +#ifdef G4VERBOSE + if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt "; +#endif + + if (parent == 0) FillParent(); + if (daughters == 0) FillDaughters(); + + // parent mass + G4double parentmass = parent->GetPDGMass(); + + EMMU = parentmass; + + //daughters'mass + G4double daughtermass[3]; + G4double sumofdaughtermass = 0.0; + for (G4int index=0; index<3; index++){ + daughtermass[index] = daughters[index]->GetPDGMass(); + sumofdaughtermass += daughtermass[index]; + } + + EMASS = daughtermass[0]; + + //create parent G4DynamicParticle at rest + G4ThreeVector dummy; + G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); + //create G4Decayproducts + G4DecayProducts *products = new G4DecayProducts(*parentparticle); + delete parentparticle; + + // calcurate electron energy + + G4double michel_rho = 0.75; //Standard Model Michel rho + G4double michel_delta = 0.75; //Standard Model Michel delta + G4double michel_xsi = 1.00; //Standard Model Michel xsi + G4double michel_eta = 0.00; //Standard Model eta + + G4double rndm, x, ctheta; + + G4double FG; + G4double FG_max = 2.00; + + G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU); + G4double x0 = EMASS/W_mue; + + G4double x0_squared = x0*x0; + + // *************************************************** + // x0 <= x <= 1. and -1 <= y <= 1 + // + // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y + // *************************************************** + + // ***** sampling F(x,y) directly (brute force) ***** + + do{ + + // Sample the positron energy by sampling from F + + rndm = G4UniformRand(); + + x = x0 + rndm*(1.-x0); + + G4double x_squared = x*x; + + G4double F_IS, F_AS, G_IS, G_AS; + + F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared); + F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared)); + + G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared); + G_IS = G_IS + michel_eta*(1.-x)*x0; + + G_AS = 3.*(michel_xsi-1.)*(1.-x); + G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared)); + G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS; + + F_IS = F_IS + G_IS; + F_AS = F_AS + G_AS; + + // *** Radiative Corrections *** + + G4double R_IS = F_c(x,x0); + + G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared); + + // *** Radiative Corrections *** + + G4double R_AS = F_theta(x,x0); + + rndm = G4UniformRand(); + + ctheta = 2.*rndm-1.; + + G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared); + + FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta); + + if(FG>FG_max){ + G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<PushProducts(daughterparticle0); + + + // daughter 1 ,2 (neutrinos) + // create neutrinos in the C.M frame of two neutrinos + G4double energy2 = parentmass*(1.0 - x/2.0); + G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0])); + G4double beta = -1.0*daughtermomentum[0]/energy2; + G4double costhetan = 2.*G4UniformRand()-1.0; + G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); + G4double phin = twopi*G4UniformRand()*rad; + G4double sinphin = std::sin(phin); + G4double cosphin = std::cos(phin); + + G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan); + G4DynamicParticle * daughterparticle1 + = new G4DynamicParticle( daughters[1], direction1*(vmass/2.)); + G4DynamicParticle * daughterparticle2 + = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.)); + + // boost to the muon rest frame + G4LorentzVector p4; + p4 = daughterparticle1->Get4Momentum(); + p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); + daughterparticle1->Set4Momentum(p4); + p4 = daughterparticle2->Get4Momentum(); + p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); + daughterparticle2->Set4Momentum(p4); + products->PushProducts(daughterparticle1); + products->PushProducts(daughterparticle2); + daughtermomentum[1] = daughterparticle1->GetTotalMomentum(); + daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); + + // output message +#ifdef G4VERBOSE + if (GetVerboseLevel()>1) { + G4cout << "G4MuonDecayChannelWithSpin::DecayIt "; + G4cout << " create decay products in rest frame " <DumpInfo(); + } +#endif + return products; +} + +G4double G4MuonDecayChannelWithSpin::R_c(G4double x){ + + G4int n_max = (int)(100.*x); + + if(n_max<10)n_max=10; + + G4double L2 = 0.0; + + for(G4int n=1; n<=n_max; n++){ + L2 += std::pow(x,n)/(n*n); + } + + G4double omega = std::log(EMMU/EMASS); + + G4double r_c; + + r_c = 2.*L2-(pi*pi/3.)-2.; + r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x)); + r_c = r_c - std::log(x)*(2.*std::log(x)-1.); + r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x); + + return r_c; +} diff --git a/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.hh b/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.hh new file mode 100644 index 0000000..b3ee2c0 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.hh @@ -0,0 +1,122 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: +// 17 August 2004 P.Gumplinger and T.MacPhail +// samples Michel spectrum including 1st order +// radiative corrections +// Reference: Florian Scheck "Muon Physics", in Physics Reports +// (Review Section of Physics Letters) 44, No. 4 (1978) +// 187-248. North-Holland Publishing Company, Amsterdam +// at page 210 cc. +// +// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25. +// +// ------------------------------------------------------------ +#ifndef G4MuonDecayChannelWithSpin_hh +#define G4MuonDecayChannelWithSpin_hh 1 + +#include "globals.hh" +#include "G4ThreeVector.hh" +#include "G4MuonDecayChannel.hh" + +class G4MuonDecayChannelWithSpin : public G4MuonDecayChannel +{ + // Class Decription + // This class describes muon decay kinemtics. + // This version assumes V-A coupling with 1st order radiative correctons, + // the standard model Michel parameter values, but + // gives incorrect energy spectrum for neutrinos + +public: // With Description + + //Constructors + G4MuonDecayChannelWithSpin(const G4String& theParentName, + G4double theBR); + // Destructor + virtual ~G4MuonDecayChannelWithSpin(); + +public: // With Description + + virtual G4DecayProducts *DecayIt(G4double); + + void SetPolarization(G4ThreeVector); + const G4ThreeVector& GetPolarization() const; + +private: + + G4ThreeVector parent_polarization; + +// Radiative Correction Factors + + G4double F_c(G4double x, G4double x0); + G4double F_theta(G4double x, G4double x0); + G4double R_c(G4double x); + + G4double EMMU; + G4double EMASS; + +}; + +inline void G4MuonDecayChannelWithSpin::SetPolarization(G4ThreeVector polar) +{ + parent_polarization = polar; +} + +inline const G4ThreeVector& G4MuonDecayChannelWithSpin::GetPolarization() const +{ + return parent_polarization; +} + +inline G4double G4MuonDecayChannelWithSpin::F_c(G4double x, G4double x0) +{ + G4double omega = std::log(EMMU/EMASS); + + G4double f_c; + + f_c = (5.+17.*x-34.*x*x)*(omega+std::log(x))-22.*x+34.*x*x; + f_c = (1.-x)/(3.*x*x)*f_c; + f_c = (6.-4.*x)*R_c(x)+(6.-6.*x)*std::log(x) + f_c; + f_c = (fine_structure_const/twopi) * (x*x-x0*x0) * f_c; + + return f_c; +} + +inline G4double G4MuonDecayChannelWithSpin::F_theta(G4double x, G4double x0) +{ + G4double omega = std::log(EMMU/EMASS); + + G4double f_theta; + + f_theta = (1.+x+34*x*x)*(omega+std::log(x))+3.-7.*x-32.*x*x; + f_theta = f_theta + ((4.*(1.-x)*(1.-x))/x)*std::log(1.-x); + f_theta = (1.-x)/(3.*x*x) * f_theta; + f_theta = (2.-4.*x)*R_c(x)+(2.-6.*x)*std::log(x)-f_theta; + f_theta = (fine_structure_const/twopi) * (x*x-x0*x0) * f_theta; + + return f_theta; +} + +#endif diff --git a/geant4/LEMuSR/G4Modified/G4Muonium.cc b/geant4/LEMuSR/G4Modified/G4Muonium.cc new file mode 100644 index 0000000..951836c --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4Muonium.cc @@ -0,0 +1,119 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4Muonium.cc,v 1.8 2003/06/16 16:57:55 gunter Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-01 $ +// +// +// ---------------------------------------------------------------------- +// GEANT 4 class implementation file +// +// History: first implementation, based on object model of +// 4th April 1996, G.Cosmo +// ********************************************************************** +// Added particle definitions, H.Kurashige, 19 April 1996 +// Added SetCuts implementation, L.Urban, 10 May 1996 +// Operators (+=, *=, ++, -> etc.) correctly used, P. Urban, 26/6/96 +// Add MuonPlusDefinition(), H.Kurashige 4 July 1996 +// Add MuoniumDefinition(),T.K.Paraiso 7 April 2005 +// ---------------------------------------------------------------------- + +#include +#include + +#include "G4Muonium.hh" + +#include "G4MuonDecayChannel.hh" +#include "G4DecayTable.hh" +#include "Randomize.hh" + +// ###################################################################### +// ### MUONPLUS ### +// ###################################################################### + +G4Muonium::G4Muonium( + const G4String& aName, G4double mass, + G4double width, G4double charge, + G4int iSpin, G4int iParity, + G4int iConjugation, G4int iIsospin, + G4int iIsospin3, G4int gParity, + const G4String& pType, G4int lepton, + G4int baryon, G4int encoding, + G4bool stable, G4double lifetime, + G4DecayTable *decaytable ) + : G4ParticleDefinition( aName,mass,width,charge,iSpin,iParity, + iConjugation,iIsospin,iIsospin3,gParity,pType, + lepton,baryon,encoding,stable,lifetime,decaytable ) +{ + SetParticleSubType("mu"); + SetPDGStable(false); + SetMuoniumGamma(); + + //create Decay Table + G4DecayTable* table = GetDecayTable(); + if (table!=NULL) delete table; + table = new G4DecayTable(); + + SetDecayTable(table); + + +} + +// ...................................................................... +// ... static member definitions ... +// ...................................................................... +// +// Arguments for constructor are as follows +// name mass width charge +// 2*spin parity C-conjugation +// 2*Isospin 2*Isospin3 G-parity +// type lepton number baryon number PDG encoding +// stable lifetime decay table +G4Muonium G4Muonium::theMuonium( + "Mu", 0.1056584*GeV, 2.99591e-16*MeV, 0.*eplus, + 1, 0, 0, + 0, 0, 0, + "lepton", -1, 0, -1313, + true, 2197.03*ns, NULL +); + +G4Muonium* G4Muonium::MuoniumDefinition() {return &theMuonium;} + +G4Muonium* G4Muonium::Muonium() +{ + return &theMuonium; +} + + +void G4Muonium::SetMuoniumGamma() +{ + G4double gamma; + gamma = 0.5*((1.*eplus)/(0.1056584*GeV/(c_light*c_light))-(1.*eplus)/(0.51099906*MeV/(c_light*c_light))); + + // SetGyromagneticRatio(gamma); +} + + + + + diff --git a/geant4/LEMuSR/G4Modified/G4Muonium.cc~ b/geant4/LEMuSR/G4Modified/G4Muonium.cc~ new file mode 100644 index 0000000..a72c26f --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4Muonium.cc~ @@ -0,0 +1,119 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4Muonium.cc,v 1.8 2003/06/16 16:57:55 gunter Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-01 $ +// +// +// ---------------------------------------------------------------------- +// GEANT 4 class implementation file +// +// History: first implementation, based on object model of +// 4th April 1996, G.Cosmo +// ********************************************************************** +// Added particle definitions, H.Kurashige, 19 April 1996 +// Added SetCuts implementation, L.Urban, 10 May 1996 +// Operators (+=, *=, ++, -> etc.) correctly used, P. Urban, 26/6/96 +// Add MuonPlusDefinition(), H.Kurashige 4 July 1996 +// Add MuoniumDefinition(),T.K.Paraiso 7 April 2005 +// ---------------------------------------------------------------------- + +#include +#include + +#include "G4Muonium.hh" + +#include "G4MuonDecayChannel.hh" +#include "G4DecayTable.hh" +#include "Randomize.hh" + +// ###################################################################### +// ### MUONPLUS ### +// ###################################################################### + +G4Muonium::G4Muonium( + const G4String& aName, G4double mass, + G4double width, G4double charge, + G4int iSpin, G4int iParity, + G4int iConjugation, G4int iIsospin, + G4int iIsospin3, G4int gParity, + const G4String& pType, G4int lepton, + G4int baryon, G4int encoding, + G4bool stable, G4double lifetime, + G4DecayTable *decaytable ) + : G4ParticleDefinition( aName,mass,width,charge,iSpin,iParity, + iConjugation,iIsospin,iIsospin3,gParity,pType, + lepton,baryon,encoding,stable,lifetime,decaytable ) +{ + SetParticleSubType("mu"); + SetPDGStable(false); + SetMuoniumGamma(); + + //create Decay Table + G4DecayTable* table = GetDecayTable(); + if (table!=NULL) delete table; + table = new G4DecayTable(); + + SetDecayTable(table); + + +} + +// ...................................................................... +// ... static member definitions ... +// ...................................................................... +// +// Arguments for constructor are as follows +// name mass width charge +// 2*spin parity C-conjugation +// 2*Isospin 2*Isospin3 G-parity +// type lepton number baryon number PDG encoding +// stable lifetime decay table +G4Muonium G4Muonium::theMuonium( + "Mu", 0.1056584*GeV, 2.99591e-16*MeV, 0.*eplus, + 1, 0, 0, + 0, 0, 0, + "lepton", -1, 0, -1313, + true, 2197.03*ns, NULL +); + +G4Muonium* G4Muonium::MuoniumDefinition() {return &theMuonium;} + +G4Muonium* G4Muonium::Muonium() +{ + return &theMuonium; +} + + +void G4Muonium::SetMuoniumGamma() +{ + G4double gamma; + gamma = 0.5*((1.*eplus)/(0.1056584*GeV/(c_light*c_light))-(1.*eplus)/(0.51099906*MeV/(c_light*c_light))); + + SetGyromagneticRatio(gamma); +} + + + + + diff --git a/geant4/LEMuSR/G4Modified/G4Muonium.hh b/geant4/LEMuSR/G4Modified/G4Muonium.hh new file mode 100644 index 0000000..102f339 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4Muonium.hh @@ -0,0 +1,90 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4Muonium.hh,v 1.7 2001/10/16 08:16:15 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-01 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 4-th April 1996, G.Cosmo +// **************************************************************** +// Added particle definitions, H.Kurashige, 19 April 1996 +// Added SetCuts, L.Urban, 12 June 1996 +// Added not static GetEnergyCuts() and GetLengthCuts(), G.Cosmo, 11 July 1996 +// Muonium defined as a Mu, T.K.Paraiso 07-04-05 +// ---------------------------------------------------------------- + +// Each class inheriting from G4VLepton +// corresponds to a particle type; one and only one +// instance for each class is guaranteed. + +#ifndef G4Muonium_h +#define G4Muonium_h 1 + +#include "globals.hh" +#include "G4ios.hh" +#include "G4ParticleDefinition.hh" + +// ###################################################################### +// ### MUONPLUS ### +// ###################################################################### + +class G4Muonium : public G4ParticleDefinition +{ +private: // constructors are hide as private + G4Muonium( + const G4String& aName, G4double mass, + G4double width, G4double charge, + G4int iSpin, G4int iParity, + G4int iConjugation, G4int iIsospin, + G4int iIsospin3, G4int gParity, + const G4String& pType, G4int lepton, + G4int baryon, G4int encoding, + G4bool stable, G4double lifetime, + G4DecayTable *decaytable + ); + +public: + virtual ~G4Muonium(){;}; + + void SetMuoniumGamma(); + + static G4Muonium* MuoniumDefinition(); + static G4Muonium* Muonium(); + +private: + static G4Muonium theMuonium; + +}; + + +#endif + + + + + + diff --git a/geant4/LEMuSR/G4Modified/G4Muonium.hh~ b/geant4/LEMuSR/G4Modified/G4Muonium.hh~ new file mode 100644 index 0000000..7b8f99a --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4Muonium.hh~ @@ -0,0 +1,90 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4Muonium.hh,v 1.7 2001/10/16 08:16:15 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-01 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 4-th April 1996, G.Cosmo +// **************************************************************** +// Added particle definitions, H.Kurashige, 19 April 1996 +// Added SetCuts, L.Urban, 12 June 1996 +// Added not static GetEnergyCuts() and GetLengthCuts(), G.Cosmo, 11 July 1996 +// Muonium defined as a Mu0, T.K.Paraiso 07-04-05 +// ---------------------------------------------------------------- + +// Each class inheriting from G4VLepton +// corresponds to a particle type; one and only one +// instance for each class is guaranteed. + +#ifndef G4Muonium_h +#define G4Muonium_h 1 + +#include "globals.hh" +#include "G4ios.hh" +#include "G4ParticleDefinition.hh" + +// ###################################################################### +// ### MUONPLUS ### +// ###################################################################### + +class G4Muonium : public G4ParticleDefinition +{ + private: + static G4Muonium theMuonium; + + private: // constructors are hide as private + G4Muonium( + const G4String& aName, G4double mass, + G4double width, G4double charge, + G4int iSpin, G4int iParity, + G4int iConjugation, G4int iIsospin, + G4int iIsospin3, G4int gParity, + const G4String& pType, G4int lepton, + G4int baryon, G4int encoding, + G4bool stable, G4double lifetime, + G4DecayTable *decaytable + ); + + public: + virtual ~G4Muonium(){} + + void SetMuoniumGamma(); + + + static G4Muonium* MuoniumDefinition(); + static G4Muonium* Muonium(); +}; + + +#endif + + + + + + diff --git a/geant4/LEMuSR/G4Modified/G4ParticleDefinition.cc b/geant4/LEMuSR/G4Modified/G4ParticleDefinition.cc new file mode 100644 index 0000000..fb4819b --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4ParticleDefinition.cc @@ -0,0 +1,290 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4ParticleDefinition.cc,v 1.23 2005/01/14 03:00:39 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-patch-01 $ +// +// +// -------------------------------------------------------------- +// GEANT 4 class implementation file +// +// History: first implementation, based on object model of +// 2nd December 1995, G.Cosmo +// ---------------- G4ParticleDefinition ----------------- +// first implementation by Makoto Asai, 29 January 1996 +// revised by G.Cosmo, 29 February 1996 +// revised by H.Kurashige, 19 April 1996 +// Code uses operators (+=, *=, ++, -> etc.) correctly, P. Urban, 26/6/96 +// revised by H.Kurashige, 4 July 1996 +// revised by H.Kurashige, 16 Feb 1997 +// revised by H.Kurashige, 10 Nov 1997 +// remove new/delete G4ProcessManager by H.Kurashige 06 June 1998 +// added Resonance flag and ApplyCuts flag H.Kurashige 27 June 1998 +// modify FillQuarkContents() for quarks/diquarks H.Kurashige 30 June 1998 +// modify encoding rule H.Kurashige 23 Oct. 98 +// modify FillQuarkContents() for deltas 25 Nov.,98 H.Kurashige +// +// modify FillQuarkContents() to use G4PDGCodeChecker 17 Aug. 99 H.Kurashige +// added Gyromagnetic Factors T.K.Paraïso Okt. 05 +// -------------------------------------------------------------- + + +#include "G4ParticleDefinition.hh" +#include "G4ParticleTable.hh" +#include "G4DecayTable.hh" +#include "G4PDGCodeChecker.hh" + +G4ParticleDefinition::G4ParticleDefinition( + const G4String& aName, + G4double mass, + G4double width, + G4double charge, + G4int iSpin, + G4int iParity, + G4int iConjugation, + G4int iIsospin, + G4int iIsospin3, + G4int gParity, + const G4String& pType, + G4int lepton, + G4int baryon, + G4int encoding, + G4bool stable, + G4double lifetime, + G4DecayTable *decaytable, + G4bool shortlived, + const G4String& subType, + G4int anti_encoding, + G4double gFactor) + + : theParticleName(aName), + thePDGMass(mass), + thePDGWidth(width), + thePDGCharge(charge), + thePDGiSpin(iSpin), + thePDGSpin(iSpin*0.5), + thePDGiParity(iParity), + thePDGiConjugation(iConjugation), + thePDGiGParity(gParity), + thePDGiIsospin(iIsospin), + thePDGiIsospin3(iIsospin3), + thePDGIsospin(iIsospin*0.5), + thePDGIsospin3(iIsospin3*0.5), + theLeptonNumber(lepton), + theBaryonNumber(baryon), + theParticleType(pType), + theParticleSubType(subType), + thePDGEncoding(encoding), + theAntiPDGEncoding(-1*encoding), + fShortLivedFlag(shortlived), + thePDGStable(stable), + thePDGLifeTime(lifetime), + theDecayTable(decaytable), + theProcessManager(0), + theAtomicNumber(0), + theAtomicMass(0), + theGyromagneticFactor(gFactor),//absolute value tao + theGyromagneticRatio(gFactor*charge/(2*mass/(c_light*c_light))), + verboseLevel(1), + fApplyCutsFlag(false) +{ + // tao + + // check name and register this particle into ParticleTable + theParticleTable = G4ParticleTable::GetParticleTable(); + theParticleTable->Insert(this); + + if (anti_encoding !=0) theAntiPDGEncoding = anti_encoding; + + // check quark contents +#ifdef G4VERBOSE + if (this->FillQuarkContents() != thePDGEncoding) { + if (verboseLevel>0) { + // cerr bnot G4cerr is used intentionally + // because G4ParticleDefinition constructor may be called + // before G4cerr object is instantiated !! + G4cerr << "Particle " << aName << " has a strange PDGEncoding " <theParticleName == right.theParticleName); +} + +G4int G4ParticleDefinition::operator!=(const G4ParticleDefinition &right) const +{ + return (this->theParticleName != right.theParticleName); +} + + + +G4int G4ParticleDefinition::FillQuarkContents() + // calculate quark and anti-quark contents + // return value is PDG encoding for this particle. + // It means error if the return value is differnt from + // this->thePDGEncoding. +{ + G4int flavor; + for (flavor= 0; flavor1) { + G4cout << " illegal charge (" << thePDGCharge/eplus; + G4cout << " PDG code=" << thePDGEncoding <1) { + G4cout << " illegal SPIN (" << thePDGiSpin << "/2"; + G4cout << " PDG code=" << thePDGEncoding <GetAntiPDGEncoding() << "]"<< G4endl; + G4cout << " Mass [GeV/c2] : " << thePDGMass/GeV ; + G4cout << " Width : " << thePDGWidth/GeV << G4endl; + G4cout << " Lifetime [nsec] : " << thePDGLifeTime/ns << G4endl; + G4cout << " Charge [e]: " << thePDGCharge/eplus << G4endl; + G4cout << " Spin : " << thePDGiSpin << "/2" << G4endl; + G4cout << " Parity : " << thePDGiParity << G4endl; + G4cout << " Charge conjugation : " << thePDGiConjugation << G4endl; + G4cout << " Isospin : (I,Iz): (" << thePDGiIsospin <<"/2"; + G4cout << " , " << thePDGiIsospin3 << "/2 ) " << G4endl; + G4cout << " GParity : " << thePDGiGParity << G4endl; + G4cout << " Quark contents (d,u,s,c,b,t) : " << theQuarkContent[0]; + G4cout << ", " << theQuarkContent[1]; + G4cout << ", " << theQuarkContent[2]; + G4cout << ", " << theQuarkContent[3]; + G4cout << ", " << theQuarkContent[4]; + G4cout << ", " << theQuarkContent[5] << G4endl; + G4cout << " AntiQuark contents : " << theAntiQuarkContent[0]; + G4cout << ", " << theAntiQuarkContent[1]; + G4cout << ", " << theAntiQuarkContent[2]; + G4cout << ", " << theAntiQuarkContent[3]; + G4cout << ", " << theAntiQuarkContent[4]; + G4cout << ", " << theAntiQuarkContent[5] << G4endl; + G4cout << " Lepton number : " << theLeptonNumber; + G4cout << " Baryon number : " << theBaryonNumber << G4endl; + G4cout << " Particle type : " << theParticleType ; + G4cout << " [" << theParticleSubType << "]" << G4endl; + + if ( fShortLivedFlag ){ + G4cout << " ShortLived : ON" << G4endl; + } + + if ( thePDGStable ){ + G4cout << " Stable : stable" << G4endl; + } else { + if( theDecayTable != 0 ){ + theDecayTable->DumpInfo(); + } else { + G4cout << "Decay Table is not defined !!" < + +class G4ProcessManager; +class G4DecayTable; +class G4ParticleTable; +class G4ParticlePropertyTable; + +class G4ParticleDefinition +{ + // Class Description + // This class containes all the static data of a particle. + // It also has uses a process manager in order to collect + // all the processes this kind of particle can undertake. + // + + friend class G4ParticlePropertyTable; + + public: // With Description + // Only one type of constructor can be used for G4ParticleDefinition. + // If you want to create new particle, you must set name of the particle + // at construction. Most of members seen as arguments of the constructor + // (except last 3 arguments concerning with decay ) are "constant" + // and can not be changed later. (No "SET" methods are available) + // Each type of particle must be constructed as a unique static object + // of special class derived from G4ParticleDefinition. + // see G4ParticleTypes for detail + + G4ParticleDefinition(const G4String& aName, + G4double mass, + G4double width, + G4double charge, + G4int iSpin, + G4int iParity, + G4int iConjugation, + G4int iIsospin, + G4int iIsospinZ, + G4int gParity, + const G4String& pType, + G4int lepton, + G4int baryon, + G4int encoding, + G4bool stable, + G4double lifetime, + G4DecayTable *decaytable, + G4bool shortlived = false, + const G4String& subType ="", + G4int anti_encoding =0, + G4double gFactor=2.0023193); + + virtual ~G4ParticleDefinition(); + + public: // With Description + // By these following Getxxxx methods, you can get values + // for members which can not be changed + // + const G4String& GetParticleName() const { return theParticleName; } + + G4double GetPDGMass() const { return thePDGMass; } + G4double GetPDGWidth() const { return thePDGWidth; } + G4double GetPDGCharge() const { return thePDGCharge; } + + G4double GetPDGSpin() const { return thePDGSpin; } + G4int GetPDGiSpin() const { return thePDGiSpin; } + G4int GetPDGiParity() const { return thePDGiParity; } + G4int GetPDGiConjugation() const { return thePDGiConjugation; } + G4double GetPDGIsospin() const { return thePDGIsospin; } + G4double GetPDGIsospin3() const { return thePDGIsospin3; } + G4int GetPDGiIsospin() const { return thePDGiIsospin; } + G4int GetPDGiIsospin3() const { return thePDGiIsospin3; } + G4int GetPDGiGParity() const { return thePDGiGParity; } + + const G4String& GetParticleType() const { return theParticleType; } + const G4String& GetParticleSubType() const { return theParticleSubType; } + G4int GetLeptonNumber() const { return theLeptonNumber; } + G4int GetBaryonNumber() const { return theBaryonNumber; } + + G4int GetPDGEncoding() const { return thePDGEncoding; } + G4int GetAntiPDGEncoding() const { return theAntiPDGEncoding; } + void SetAntiPDGEncoding(G4int aEncoding); + + + G4int GetQuarkContent(G4int flavor) const; + G4int GetAntiQuarkContent(G4int flavor) const; + // return the number of quark with flavor contained in this particle. + // The value of flavor is assigned as follows + // 1:d, 2:u, 3:s, 4:c, 5:b, 6:t + + public: // With Description + // ShortLived flag + G4bool IsShortLived() const { return fShortLivedFlag; } + + G4bool GetPDGStable() const { return thePDGStable; } + void SetPDGStable(const G4bool aFlag) { thePDGStable=aFlag; } + + G4double GetPDGLifeTime() const { return thePDGLifeTime; } + void SetPDGLifeTime(G4double aLifeTime) { thePDGLifeTime = aLifeTime; } + + public:// With Description + G4DecayTable* GetDecayTable(); + void SetDecayTable(G4DecayTable* aDecayTable); + // Set/Get Decay Table + // !! Decay Table can be modified !! + + public: // With Description + G4ProcessManager* GetProcessManager() const; + void SetProcessManager(G4ProcessManager* aProcessManager); + // Set/Get Process Manager + // !! Process Manager can be modified !! + + G4ParticleTable* GetParticleTable(); + // get pointer to the particle table + + void DumpTable() const; + // Prints information of data members. + + protected: + G4int FillQuarkContents(); + // calculate quark and anti-quark contents + // return value is PDG encoding for this particle. + // It means error if the return value is deffernt from + // this->thePDGEncoding. + + void SetParticleSubType(const G4String& subtype); + + public: + void SetAtomicNumber(G4int); + G4int GetAtomicNumber() const; + void SetAtomicMass(G4int); + G4int GetAtomicMass() const; + + G4double GetGyromagneticFactor() const { return theGyromagneticFactor; } + void SetGyromagneticFactor(G4double ge) { theGyromagneticFactor = ge; } + + G4double GetGyromagneticRatio() const { return theGyromagneticRatio; } + void SetGyromagneticRatio(G4double gamma) { theGyromagneticRatio = gamma; } + + public: + void SetVerboseLevel(G4int value); + G4int GetVerboseLevel() const; + // controle flag for output message + // 0: Silent + // 1: Warning message + // 2: More + + protected: + // !!! can not use "copy constructor" nor "default constructor" !!!! + G4ParticleDefinition(const G4ParticleDefinition &right); + G4ParticleDefinition(); + + private: + // !!! Assignment operation is forbidden !!! + const G4ParticleDefinition & operator=(const G4ParticleDefinition &right); + + public: + G4int operator==(const G4ParticleDefinition &right) const; + G4int operator!=(const G4ParticleDefinition &right) const; + + private: + // Values following can not be changed + // i.e. No Setxxxx Methods for them + + G4String theParticleName; + // The name of the particle. + // Each object must have its specific name!! + + // --- following member values must be defined with Units + G4double thePDGMass; + // The mass of the particle, in units of equivalent energy. + + G4double thePDGWidth; + // The decay width of the particle, usually the width of a + // Breit-Wigner function, assuming that you are near the + // mass center anyway. (in units of equivalent energy) + + G4double thePDGCharge; + // The charge of the particle.(in units of Coulomb) + + // ---- following members are quantum number + // i.e. discrete numbers can be allowded + // So, you can defined only by using integer in constructor + G4int thePDGiSpin; + // The total spin of the particle, also often denoted as + // capital J, in units of 1/2. + G4double thePDGSpin; + // The total spin of the particle, in units of 1. + + G4int thePDGiParity; + // The parity quantum number, in units of 1. If the parity + // is not defined for this particle, we will set this to 0. + + G4int thePDGiConjugation; + // This charge conjugation quantum number in units of 1. + + G4int thePDGiGParity; + // The value of the G-parity quantum number. + + G4int thePDGiIsospin; + G4int thePDGiIsospin3; + // The isospin and its 3rd-component in units of 1/2. + G4double thePDGIsospin; + G4double thePDGIsospin3; + // The isospin quantum number in units of 1. + + G4int theLeptonNumber; + // The lepton quantum number. + + G4int theBaryonNumber; + // The baryon quantum number. + + G4String theParticleType; + // More general textual type description of the particle. + + G4String theParticleSubType; + // Textual type description of the particle + // eg. pion, lamda etc. + + G4int thePDGEncoding; + // The Particle Data Group integer identifier of this particle + + G4int theAntiPDGEncoding; + // The Particle Data Group integer identifier of the anti-particle + + protected: + enum {NumberOfQuarkFlavor = 6}; + G4int theQuarkContent[NumberOfQuarkFlavor]; + G4int theAntiQuarkContent[NumberOfQuarkFlavor]; + // the number of quark (minus Sign means anti-quark) contents + // The value of flavor is assigned as follows + // 0:d, 1:u, 2:s, 3:c, 4:b, 5:t + + private: + // Following members can be changed after construction + + G4bool fShortLivedFlag; + // Particles which have true value of this flag + // will not be tracked by TrackingManager + + G4bool thePDGStable; + // Is an indicator that this particle is stable. It must + // not decay. If the user tries to assign a kind of decay + // object to it, it will refuse to take it. + + G4double thePDGLifeTime; + // Is related to the decay width of the particle. The mean + // life time is given in seconds. + + class G4DecayTable *theDecayTable; + // Points DecayTable + + private: + class G4ProcessManager *theProcessManager; + // Points to G4ProcessManager + + G4ParticleTable* theParticleTable; + + private: + G4int theAtomicNumber; + G4int theAtomicMass; + + G4double theGyromagneticFactor; + // The gyromanetic factor of the particle + G4double theGyromagneticRatio; + // The gyromanetic ratio of the particle + + private: + G4int verboseLevel; + + private: + G4bool fApplyCutsFlag; + public: + + void SetApplyCutsFlag(G4bool); + G4bool GetApplyCutsFlag() const; + +}; + +#include "G4ParticleDefinition.icc" + +#endif diff --git a/geant4/LEMuSR/G4Modified/G4ParticleDefinition.icc b/geant4/LEMuSR/G4Modified/G4ParticleDefinition.icc new file mode 100644 index 0000000..9046be6 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4ParticleDefinition.icc @@ -0,0 +1,122 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4ParticleDefinition.icc,v 1.8 2004/12/02 08:08:58 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-03 $ +// + +inline + G4ParticleTable* G4ParticleDefinition::GetParticleTable() +{ + return theParticleTable; +} + +inline + G4DecayTable* G4ParticleDefinition::GetDecayTable() +{ + return theDecayTable; +} + +inline + void G4ParticleDefinition::SetDecayTable(G4DecayTable* aDecayTable) +{ + theDecayTable = aDecayTable; +} + + +inline + void G4ParticleDefinition::SetVerboseLevel(G4int value) +{ + verboseLevel = value; +} + +inline + G4int G4ParticleDefinition::GetVerboseLevel() const +{ + return verboseLevel; +} + +inline + G4ProcessManager* G4ParticleDefinition::GetProcessManager() const +{ + return theProcessManager; +} + +inline + void G4ParticleDefinition::SetProcessManager(G4ProcessManager *aProcessManager) +{ + theProcessManager = aProcessManager; +} + +inline + G4int G4ParticleDefinition::GetQuarkContent(G4int flavor) const +{ + G4int content = 0; + if ((flavor>0) && (flavor<=NumberOfQuarkFlavor)){ + content = theQuarkContent[flavor-1]; + }else { +#ifdef G4VERBOSE + if (verboseLevel >0) { + G4cout << "Invalid Quark Flavor for G4ParticleDefinition::GetQuarkContent"; + G4cout << ": flavor=" << flavor <0) && (flavor<=NumberOfQuarkFlavor)){ + content = theAntiQuarkContent[flavor-1]; + }else { +#ifdef G4VERBOSE + if (verboseLevel >0) { + G4cout <<"Invalid Quark Flavor for G4ParticleDefinition::GetAntiQuarkContent"; + G4cout << ": flavor=" << flavor <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;index0) { + 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) ) { +#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]<GetParticleName()); +} + +void G4VDecayChannel::FillDaughters() +{ + G4int index; + +#ifdef G4VERBOSE + if (verboseLevel>1) G4cout << "G4VDecayChannel::FillDaughters()" <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 " <1) { + G4cout << " parent:" << *parent_name; + G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <GetPDGMass()/GeV; + G4cout << "[GeV/c/c]" <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; +} diff --git a/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh b/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh new file mode 100644 index 0000000..01e18c2 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh @@ -0,0 +1,288 @@ +// +// ******************************************************************** +// * DISCLAIMER * +// * * +// * The following disclaimer summarizes all the specific disclaimers * +// * of contributors to this software. The specific disclaimers,which * +// * govern, are listed with their locations in: * +// * http://cern.ch/geant4/license * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. * +// * * +// * This code implementation is the intellectual property of the * +// * GEANT4 collaboration. * +// * By copying, distributing or modifying the Program (or any work * +// * based on the Program) you indicate your acceptance of this * +// * statement, and all its terms. * +// ******************************************************************** +// +// +// $Id: G4VDecayChannel.hh,v 1.11 2004/12/02 08:08:58 kurasige Exp $ +// GEANT4 tag $Name: geant4-07-00-cand-03 $ +// +// +// ------------------------------------------------------------ +// GEANT 4 class header file +// +// History: first implementation, based on object model of +// 27 July 1996 H.Kurashige +// 30 May 1997 H.Kurashige +// 23 Mar. 2000 H.Weber : add GetAngularMomentum() +// ------------------------------------------------------------ +#ifndef G4VDecayChannel_h +#define G4VDecayChannel_h 1 + +#include "G4ios.hh" +#include "globals.hh" +#include +#include "G4ThreeVector.hh" + +class G4ParticleDefinition; +class G4DecayProducts; +class G4ParticleTable; + +class G4VDecayChannel +{ + // Class Description + // This class is a abstract class to describe decay kinematics + // + + public: + //Constructors + G4VDecayChannel(const G4String &aName, G4int Verbose = 1); + G4VDecayChannel(const G4String &aName, + const G4String& theParentName, + G4double theBR, + G4int theNumberOfDaughters, + const G4String& theDaughterName1, + const G4String& theDaughterName2 = "", + const G4String& theDaughterName3 = "", + const G4String& theDaughterName4 = "" ); + + // Destructor + virtual ~G4VDecayChannel(); + + private: + // copy constructor and assignment operatotr + G4VDecayChannel(const G4VDecayChannel &); + G4VDecayChannel & operator=(const G4VDecayChannel &); + + public: + // equality operators + G4int operator==(const G4VDecayChannel &right) const {return (this == &right);} + G4int operator!=(const G4VDecayChannel &right) const {return (this != &right);} + + // less-than operator is defined for G4DecayTable + G4int operator<(const G4VDecayChannel &right) const; + + public: // With Description + virtual G4DecayProducts* DecayIt(G4double parentMass = -1.0) = 0; + + public: // With Description + //get kinematics name + G4String GetKinematicsName() const; + //get branching ratio + G4double GetBR() const; + //get number of daughter particles + G4int GetNumberOfDaughters() const; + + //get the pointer to the parent particle + G4ParticleDefinition * GetParent(); + //get the pointer to a daughter particle + G4ParticleDefinition * GetDaughter(G4int anIndex); + + //get the angular momentum of the decay + G4int GetAngularMomentum(); + //get the name of the parent particle + const G4String& GetParentName() const; + //get the name of a daughter particle + const G4String& GetDaughterName(G4int anIndex) const; + + // get mass of parent + G4double GetParentMass() const; + G4double GetDaughterMass(G4int anIndex) const; + G4ThreeVector GetParentPolarization(); + //set the parent particle (by name or by pointer) + void SetParent(const G4ParticleDefinition * particle_type); + void SetParent(const G4String &particle_name); + //set branching ratio + void SetBR(G4double value); + //set number of daughter particles + void SetNumberOfDaughters(G4int value); + //set a daughter particle (by name or by pointer) + void SetDaughter(G4int anIndex, + const G4ParticleDefinition * particle_type); + void SetDaughter(G4int anIndex, + const G4String &particle_name); + void SetParentPolarization(G4ThreeVector polar); + protected: + // kinematics name + G4String kinematics_name; + // branching ratio [0.0 - 1.0] + G4double rbranch; + // number of daughters + G4int numberOfDaughters; + // parent particle + G4String* parent_name; + //daughter particles + G4String** daughters_name; + + protected: // With Description + // celar daughters array + void ClearDaughtersName(); + + protected: + // pointer to particle table + G4ParticleTable* particletable; + + // temporary buffers of pointers to G4ParticleDefinition + G4ParticleDefinition* parent; + G4ParticleDefinition** daughters; + + // parent mass + G4double parent_mass; + G4double* daughters_mass; + + //parent theParentPolarization + G4ThreeVector theParentPolarization; + + // fill daughters array + void FillDaughters(); + // fill parent + void FillParent(); + + public: // With Description + void SetVerboseLevel(G4int value); + G4int GetVerboseLevel() const; + void DumpInfo(); + + private: + const G4String& GetNoName() const; + + private: + // controle flag for output message + G4int verboseLevel; + // 0: Silent + // 1: Warning message + // 2: More + + static const G4String noName; +}; + +inline + G4int G4VDecayChannel::operator<(const G4VDecayChannel &right) const +{ + return (this->rbranch < right.rbranch); +} + +inline + G4ParticleDefinition* G4VDecayChannel::GetDaughter(G4int anIndex) + { + //pointers to daughter particles are filled, if they are not set yet + if (daughters == 0) FillDaughters(); + + //get the pointer to a daughter particle + if ( (anIndex>=0) && (anIndex0) + G4cout << "G4VDecayChannel::GetDaughter index out of range "<=0) && (anIndex0){ + G4cout << "G4VDecayChannel::GetDaughterName "; + G4cout << "index out of range " << anIndex << G4endl; + } + return GetNoName(); + } +} + +inline + G4double G4VDecayChannel::GetDaughterMass(G4int anIndex) const +{ + if ( (anIndex>=0) && (anIndex0){ + 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 + + + + + + diff --git a/geant4/LEMuSR/G4Modified/Readme.txt b/geant4/LEMuSR/G4Modified/Readme.txt new file mode 100644 index 0000000..d341a20 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/Readme.txt @@ -0,0 +1,11 @@ +Geant4 ElectricField update files: + +Copy files in the following directories: + +G4El_EqRhs -----> $G4Install/source/geometry/magneticfield /src and /include +G4El_UsualEqRhs -----> $G4Install/source/geometry/magneticfield /src and /include +G4ChordFinder -----> $G4Install/source/geometry/magneticfield /src and /include +G4ParticleGun -----> $G4Install/source/event /src and /include + + +track/ is the geant4.6.2 version of $G4Install/source/track directory and should be updated \ No newline at end of file diff --git a/geant4/LEMuSR/G4Modified/update.sh b/geant4/LEMuSR/G4Modified/update.sh new file mode 100644 index 0000000..74dee34 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/update.sh @@ -0,0 +1,57 @@ + + + +cp G4ChordFinder.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4ChordFinder.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ + +cp G4El_EqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4El_EqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ + +cp G4El_MagEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4El_MagEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ + +cp G4El_UsualEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4El_UsualEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4FieldManager.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4FieldManager.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4FieldManager.icc ~/geant4.8.0.p01/source/geometry/magneticfield/include/ + +cp G4MultipleScattering52.cc ~/geant4.8.0.p01/source/processes/electromagnetic/standard/src/ + +cp G4VDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4VDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4MuonDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4MuonDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4MuonDecayChannelWithSpin.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4MuonDecayChannelWithSpin.hh ~/geant4.8.0.p01/source/particles/management/include/ + +cp G4Muonium.cc ~/geant4.8.0.p01/source/particles/leptons/src/ +cp G4Muonium.hh ~/geant4.8.0.p01/source/particles/leptons/include/ +cp G4ParticleDefinition.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4ParticleDefinition.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4ParticleDefinition.icc ~/geant4.8.0.p01/source/particles/management/include/ + + +cd ~/geant4.8.0.p01/source/geometry/magneticfield/ +#gmake clean +gmake + +cd ~/geant4.8.0.p01/source/processes/electromagnetic/standard/ +#gmake clean +gmake + +cd ~/geant4.8.0.p01/source/particles/management/ +#gmake clean +gmake + +cd ~/geant4.8.0.p01/source/particles/leptons/ +#gmake clean +gmake + + + + +cd ~/LEMuSR/G4Modified + + + diff --git a/geant4/LEMuSR/G4Modified/update.sh~ b/geant4/LEMuSR/G4Modified/update.sh~ new file mode 100644 index 0000000..db882c7 --- /dev/null +++ b/geant4/LEMuSR/G4Modified/update.sh~ @@ -0,0 +1,47 @@ + + + +cp G4ChordFinder.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4ChordFinder.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4El_EqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4El_EqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4El_MagEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4El_MagEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4El_UsualEqRhs.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4El_UsualEqRhs.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4FieldManager.cc ~/geant4.8.0.p01/source/geometry/magneticfield/src/ +cp G4FieldManager.hh ~/geant4.8.0.p01/source/geometry/magneticfield/include/ +cp G4MultipleScattering52.cc ~/geant4.8.0.p01/source/processes/electromagnetic/standard/src/ +cp G4VDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4VDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4MuonDecayChannel.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4MuonDecayChannel.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4MuonDecayChannelWithSpin.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4MuonDecayChannelWithSpin.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4Muonium.cc ~/geant4.8.0.p01/source/particles/leptons/src/ +cp G4Muonium.hh ~/geant4.8.0.p01/source/particles/leptons/include/ +cp G4ParticleDefinition.cc ~/geant4.8.0.p01/source/particles/management/src/ +cp G4ParticleDefinition.hh ~/geant4.8.0.p01/source/particles/management/include/ +cp G4ParticleDefinition.icc ~/geant4.8.0.p01/source/particles/management/include/ + + +cd ~/geant4.8.0.p01/source/geometry/magneticfield/ +#gmake clean +gmake + +cd ~/geant4.8.0.p01/source/processes/electromagnetic/standard/ +#gmake clean +gmake + +cd ~/geant4.8.0.p01/source/particles/management/ +#gmake clean +gmake + +cd ~/geant4.8.0.p01/source/particles/leptons/ +#gmake clean +gmake + + + + +cd ~/LEMuSR/G4Modified \ No newline at end of file