diff --git a/geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc b/geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc deleted file mode 100644 index e7524a9..0000000 --- a/geant4/LEMuSR/G4Modified/G4MultipleScattering52.cc +++ /dev/null @@ -1,951 +0,0 @@ -// -// ******************************************************************** -// * 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/G4MuonDecayChannelWithSpin.cc b/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc deleted file mode 100644 index eb00503..0000000 --- a/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.cc +++ /dev/null @@ -1,264 +0,0 @@ -// -// ******************************************************************** -// * 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 deleted file mode 100644 index b3ee2c0..0000000 --- a/geant4/LEMuSR/G4Modified/G4MuonDecayChannelWithSpin.hh +++ /dev/null @@ -1,122 +0,0 @@ -// -// ******************************************************************** -// * 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 deleted file mode 100644 index 951836c..0000000 --- a/geant4/LEMuSR/G4Modified/G4Muonium.cc +++ /dev/null @@ -1,119 +0,0 @@ -// -// ******************************************************************** -// * 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 deleted file mode 100644 index 102f339..0000000 --- a/geant4/LEMuSR/G4Modified/G4Muonium.hh +++ /dev/null @@ -1,90 +0,0 @@ -// -// ******************************************************************** -// * 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/G4VDecayChannel.cc b/geant4/LEMuSR/G4Modified/G4VDecayChannel.cc deleted file mode 100644 index 8475c50..0000000 --- a/geant4/LEMuSR/G4Modified/G4VDecayChannel.cc +++ /dev/null @@ -1,460 +0,0 @@ -// -// ******************************************************************** -// * DISCLAIMER * -// * * -// * The following disclaimer summarizes all the specific disclaimers * -// * of contributors to this software. The specific disclaimers,which * -// * govern, are listed with their locations in: * -// * http://cern.ch/geant4/license * -// * * -// * Neither the authors of this software system, nor their employing * -// * institutes,nor the agencies providing financial support for this * -// * work make any representation or warranty, express or implied, * -// * regarding this software system or assume any liability for its * -// * use. * -// * * -// * This code implementation is the intellectual property of the * -// * GEANT4 collaboration. * -// * By copying, distributing or modifying the Program (or any work * -// * based on the Program) you indicate your acceptance of this * -// * statement, and all its terms. * -// ******************************************************************** -// -// -// $Id: G4VDecayChannel.cc,v 1.16 2004/12/02 08:09:00 kurasige Exp $ -// GEANT4 tag $Name: geant4-07-00-cand-03 $ -// -// -// ------------------------------------------------------------ -// GEANT 4 class header file -// -// History: first implementation, based on object model of -// 27 July 1996 H.Kurashige -// 30 May 1997 H.Kurashige -// 23 Mar. 2000 H.Weber : add GetAngularMomentum -// ------------------------------------------------------------ - -#include "G4ParticleDefinition.hh" -#include "G4ParticleTable.hh" -#include "G4DecayTable.hh" -#include "G4DecayProducts.hh" -#include "G4VDecayChannel.hh" - -const G4String G4VDecayChannel::noName = " "; - -G4VDecayChannel::G4VDecayChannel(const G4String &aName, G4int Verbose) - :kinematics_name(aName), - rbranch(0.0), - numberOfDaughters(0), - parent_name(0), daughters_name(0), - particletable(0), - parent(0), daughters(0), - parent_mass(0.0), daughters_mass(0), - verboseLevel(Verbose) -{ - // set pointer to G4ParticleTable (static and singleton object) - particletable = G4ParticleTable::GetParticleTable(); -} - -G4VDecayChannel::G4VDecayChannel(const G4String &aName, - const G4String& theParentName, - G4double theBR, - G4int theNumberOfDaughters, - const G4String& theDaughterName1, - const G4String& theDaughterName2, - const G4String& theDaughterName3, - const G4String& theDaughterName4 ) - :kinematics_name(aName), - rbranch(theBR), - numberOfDaughters(theNumberOfDaughters), - parent_name(0), daughters_name(0), - particletable(0), - parent(0), daughters(0), - parent_mass(0.0), daughters_mass(0), - verboseLevel(1) -{ - // set pointer to G4ParticleTable (static and singleton object) - particletable = G4ParticleTable::GetParticleTable(); - - // parent name - parent_name = new G4String(theParentName); - - // cleate array - daughters_name = new G4String*[numberOfDaughters]; - for (G4int index=0;index0) 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 deleted file mode 100644 index 01e18c2..0000000 --- a/geant4/LEMuSR/G4Modified/G4VDecayChannel.hh +++ /dev/null @@ -1,288 +0,0 @@ -// -// ******************************************************************** -// * 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/LEMuSR.cc b/geant4/LEMuSR/LEMuSR.cc index fa8915b..9a4ef1c 100644 --- a/geant4/LEMuSR/LEMuSR.cc +++ b/geant4/LEMuSR/LEMuSR.cc @@ -1,6 +1,6 @@ //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//* // LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION Geant4 SIMULATION -// ID : LEMuSR.cc , v 1.0 +// ID : LEMuSR.cc , v 1.2 // AUTHOR: Taofiq PARAISO // DATE : 2004-06-24 09:57 @@ -28,7 +28,7 @@ #include "G4UIterminal.hh" #include "G4UItcsh.hh" #include "G4ios.hh" -#include +#include // II LEMuSR CLASSES // a_ Mandatory Classes @@ -81,15 +81,17 @@ int main(int argc,char** argv)//argc:: defines the user interface LEMuSRPhysicsList* lemuPhysicsList = new LEMuSRPhysicsList(); - //! 2.2 LEMuSR Action class - LEMuSRPrimaryGeneratorAction* lemuPGA = new LEMuSRPrimaryGeneratorAction(); - - - //! 2.3 Setting the mandatory Initialization classes + + //! 2.2 Setting the mandatory Initialization classes runManager ->SetUserInitialization( lemuDetector ); runManager ->SetUserInitialization( lemuPhysicsList ); - + + + //! 2.3 LEMuSR Action class + LEMuSRPrimaryGeneratorAction* lemuPGA = new LEMuSRPrimaryGeneratorAction(); + + //! 2.4 Setting the mandatory Action class runManager ->SetUserAction( lemuPGA ); @@ -97,7 +99,7 @@ int main(int argc,char** argv)//argc:: defines the user interface //! 3 The optionnal classes runManager ->SetUserAction( new LEMuSRRunAction()); runManager ->SetUserAction( new LEMuSREventAction());// scintillators, sensitive detectors - runManager ->SetUserAction( new LEMuSRTrackingAction()); + runManager ->SetUserAction( new LEMuSRTrackingAction()); //! Optionnal stepping action: enable one at once only @@ -131,7 +133,7 @@ int main(int argc,char** argv)//argc:: defines the user interface G4UIsession* session = 0; - + #if defined G4UI_USE_ROOT // G4UIRoot is a ROOT based GUI @@ -143,17 +145,16 @@ int main(int argc,char** argv)//argc:: defines the user interface session = new G4UIterminal(); #endif - // UI->ApplyCommand("/control/execute visual.mac"); - UI->ApplyCommand("/run/verbose 2"); + UI->ApplyCommand("/vis/sceneHandler/create"); + UI->ApplyCommand("/control/verbose 1"); + UI->ApplyCommand("/run/verbose 1"); - // UI->ApplyCommand("/Detector/MagneticField on"); - // UI->ApplyCommand("/run/beamOn 1"); #if defined LEMU_TEST_ASYM - // UI->ApplyCommand("/Detector/ElectricField off"); UI->ApplyCommand("/Detector/MagneticField 50"); UI->ApplyCommand("/Detector/AsymCheck on"); UI->ApplyCommand("/lemuGun/gunPosition 0 0 0"); UI->ApplyCommand("/lemuGun/energy/defined 0"); + UI->ApplyCommand("/lemuGun/energy/offset 0"); G4cout<<"\n READY TO TEST ASYMETRY! "; #endif diff --git a/geant4/LEMuSR/include/FocalLengthTest.hh b/geant4/LEMuSR/include/FocalLengthTest.hh index d2c6a53..e6cd908 100644 --- a/geant4/LEMuSR/include/FocalLengthTest.hh +++ b/geant4/LEMuSR/include/FocalLengthTest.hh @@ -24,7 +24,7 @@ #include "G4UserSteppingAction.hh" #include "globals.hh" -#include +#include #include "LEMuSRVisManager.hh" #include "LEMuSRDetectorConstruction.hh" diff --git a/geant4/LEMuSR/include/TDCheck.hh b/geant4/LEMuSR/include/TDCheck.hh index 06b38d9..c6d83cd 100644 --- a/geant4/LEMuSR/include/TDCheck.hh +++ b/geant4/LEMuSR/include/TDCheck.hh @@ -24,7 +24,7 @@ #include "G4UserSteppingAction.hh" #include "globals.hh" -#include +#include #include "LEMuSRVisManager.hh" #include "LEMuSRDetectorConstruction.hh" diff --git a/geant4/LEMuSR/src/LEMuSRDecay.cc b/geant4/LEMuSR/src/LEMuSRDecay.cc index 4293427..6e8a6d7 100644 --- a/geant4/LEMuSR/src/LEMuSRDecay.cc +++ b/geant4/LEMuSR/src/LEMuSRDecay.cc @@ -131,7 +131,9 @@ G4VParticleChange* LEMuSRDecay::DecayIt(const G4Track& aTrack, const G4Step& aSt testb++; } #ifdef G4VERBOSE - G4cout<<"Decay Channel LEMuSR "<1) { + G4cout<<"Decay Channel LEMuSR "<