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

4.3.2010 - Kamil Sedlak

Browse Source 
1) UPGRADE TO GEANT 4.9.3
      - Physics list had been modified
      - List of the physics processes in *.mac has to be modified too!
        (otherwise results of 4.9.3 do not reproduce well the results
         of 4.9.2).
        Example *.mac files still use the old implementation of the physics
        processes, i.e. they have to be corrected in future!

   2) POSSIBILITY TO WRITE OUT DERIVATIVES OF THE MAGNETIC FIELD 
      AT A GIVEN POINT (see the documentation,
      "/musr/command globalfield printFieldDerivativeAtPoint".)
This commit is contained in:
sedlak 2010-03-04 12:52:07 +00:00
parent 701f17857f
commit 610abf6591
8 changed files with 602 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/musrSim.pdf
View File

Binary file not shown.

25
doc/musrSim.tex
View File

@ -426,7 +426,11 @@ Three special volumes ``Target, M0, M1 and M2''.
\item{\bf /musr/command globalfield printFieldValueAtPoint \emph{x} \emph{y} \emph{z}} \\ \item{\bf /musr/command globalfield printFieldValueAtPoint \emph{x} \emph{y} \emph{z}} \\
Print out the field value at the point $(x, y, z)$ (given in the global Print out the field value at the point $(x, y, z)$ (given in the global
coordinate system. coordinate system).
\item{\bf /musr/command globalfield printFieldDerivativeAtPoint \emph{x} \emph{y} \emph{z}} \\
Print out the values of the derivative of the magnetic field at the point $(x, y, z)$
(given in the global coordinate system).
\end{description} \end{description}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -950,6 +954,25 @@ These are typically not interesting for the analysis of the results of the simul
\end{description} \end{description}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Upgrading Geant4 version}
\begin{description}
\item{\bf Geant 4.9.1 and older} - copy the file ``src/G4DecayWithSpin.cc\_for\_Geant4.9.1\_and\_older''
to the file ``src/G4DecayWithSpin.cc''.
Also copy ``src/G4EqEMFieldWithSpin.cc\_for\_Geant4.9.1\_and\_older'' to ``src/G4EqEMFieldWithSpin.cc''.
\item{\bf Geant 4.9.2.p02} -- (only this particular Geant version!) copy the file
``src/G4EqEMFieldWithSpin.cc\_for\_Geant4.9.2p02\_only'' to the file ``src/G4EqEMFieldWithSpin.cc''.
\item{\bf Geant 4.9.2 and older} - copy the file ``src/musrPhysicsList.cc\_Geant4.9.2p02\_and\_older''
to the file ``src/musrPhysicsList.cc''.
\item{\bf Geant 4.9.3}
From the point of view of musrSim, the largest change in version 4.9.3 with respect to 4.9.2
is the different treatment of the physics processes, therefore the musrPhysicsList.cc had to be
modified, and also the physics processes listed in the *.mac file have to be changed!
Note that using Geant~4.9.3 with old *.mac files will formally run OK, but the results of the
simulation (slightly) differ from the same simulation based on Geant~4.9.2. However after
updating the *.mac file, the simulation done with Geant~4.9.3 more or less reproduces results
obtained with Geant~4.9.2.
\end{description}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Example 1 -- Electrons passing through two scintillator tiles (101.mac)} \section{Example 1 -- Electrons passing through two scintillator tiles (101.mac)}
One of the easiest example to illustrate the basic features of the musrSim (and/or Geant4) is to One of the easiest example to illustrate the basic features of the musrSim (and/or Geant4) is to
shoot electrons into a scintillator block, and to observe the energy deposited inside it. shoot electrons into a scintillator block, and to observe the energy deposited inside it.

5
include/F04GlobalField.hh
View File

@ -141,9 +141,9 @@ public:
// G4bool DoesAnyFieldValueNeedsToBeChanged(G4int eventNumber) {return globalChangeFieldInStepsMap[eventNumber];} // G4bool DoesAnyFieldValueNeedsToBeChanged(G4int eventNumber) {return globalChangeFieldInStepsMap[eventNumber];}
void CheckWhetherAnyNominalFieldValueNeedsToBeChanged(G4int eventNumber); void CheckWhetherAnyNominalFieldValueNeedsToBeChanged(G4int eventNumber);
// Add point, at which user wishes to print out the field value // Add point, at which user wishes to print out the field value or field derivative
void AddPointForFieldTesting(G4ThreeVector point) {pointsAtWhichUserWantsToPrintFieldValue.push_back(point);} void AddPointForFieldTesting(G4ThreeVector point) {pointsAtWhichUserWantsToPrintFieldValue.push_back(point);}
void AddPointForFieldDerivativeTesting(G4ThreeVector point) {pointsAtWhichUserWantsToPrintFieldDerivative.push_back(point);}
// Print field value at all points user wished to be print out: // Print field value at all points user wished to be print out:
void PrintFieldAtRequestedPoints() const; void PrintFieldAtRequestedPoints() const;
@ -164,6 +164,7 @@ private:
const F04ElementField **fp; const F04ElementField **fp;
std::vector<G4ThreeVector> pointsAtWhichUserWantsToPrintFieldValue; std::vector<G4ThreeVector> pointsAtWhichUserWantsToPrintFieldValue;
std::vector<G4ThreeVector> pointsAtWhichUserWantsToPrintFieldDerivative;
private: private:

33
src/F04GlobalField.cc
View File

@ -310,7 +310,12 @@ void F04GlobalField::CheckWhetherAnyNominalFieldValueNeedsToBeChanged(G4int even
void F04GlobalField::PrintFieldAtRequestedPoints() const { void F04GlobalField::PrintFieldAtRequestedPoints() const {
G4ThreeVector p; G4ThreeVector p;
G4double point[4]; G4double point[4];
G4double delta=0.1*mm;
G4double point2[4]; G4double point3[4]; G4double point4[4];
G4double Bfi[6]={0,0,0,0,0,0}; G4double Bfi[6]={0,0,0,0,0,0};
G4double BfiX[6]={0,0,0,0,0,0};
G4double BfiY[6]={0,0,0,0,0,0};
G4double BfiZ[6]={0,0,0,0,0,0};
for (unsigned int i=0; i < pointsAtWhichUserWantsToPrintFieldValue.size(); i++) { for (unsigned int i=0; i < pointsAtWhichUserWantsToPrintFieldValue.size(); i++) {
p = pointsAtWhichUserWantsToPrintFieldValue[i]; p = pointsAtWhichUserWantsToPrintFieldValue[i];
point[0] = p.x(); point[0] = p.x();
@ -325,6 +330,34 @@ void F04GlobalField::PrintFieldAtRequestedPoints() const {
Bfi[0]/tesla, Bfi[1]/tesla, Bfi[2]/tesla, Bfi[0]/tesla, Bfi[1]/tesla, Bfi[2]/tesla,
Bfi[3]/(kilovolt/mm), Bfi[4]/(kilovolt/mm), Bfi[5]/(kilovolt/mm)); Bfi[3]/(kilovolt/mm), Bfi[4]/(kilovolt/mm), Bfi[5]/(kilovolt/mm));
} }
if (pointsAtWhichUserWantsToPrintFieldDerivative.size()>0) {
G4cout<<G4endl<<" ===="<<G4endl<<"Derivatives of the magnetic field"<<G4endl;
G4cout<<" x y z (mm) dBx/dx dBx/dy dBx/dz dBy/dx dBy/dy dBy/dz dBz/dx dBz/dy dBz/dz (G/mm)"<<G4endl;
// G4cout<<"tesla = "<<tesla<<" gauss="<<gauss<<" ==> tesla = "<<tesla/gauss<<"gauss"<<G4endl;
}
for (unsigned int i=0; i < pointsAtWhichUserWantsToPrintFieldDerivative.size(); i++) {
p = pointsAtWhichUserWantsToPrintFieldDerivative[i];
point[0] = p.x();
point[1] = p.y();
point[2] = p.z();
point[3] = 0.;
point2[0] = point[0]+delta; point2[1] = point[1]; point2[2] = point[2]; point2[3] = point[3];
point3[0] = point[0]; point3[1] = point[1]+delta; point3[2] = point[2]; point3[3] = point[3];
point4[0] = point[0]; point4[1] = point[1]; point4[2] = point[2]+delta; point4[3] = point[3];
object->GetFieldValue(point,Bfi);
object->GetFieldValue(point2,BfiX);
object->GetFieldValue(point3,BfiY);
object->GetFieldValue(point4,BfiZ);
// printf (" Magnetic Field at %f, %f, %f mm is B= %10.10f, %10.10f, %10.10f tesla.\n",
// point[0]/mm,point[1]/mm,point[2]/mm,Bfi[0]/tesla,Bfi[1]/tesla,Bfi[2]/tesla);
printf (" %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",
point[0]/mm, point[1]/mm, point[2]/mm,
(BfiX[0]-Bfi[0])/gauss/delta*mm, (BfiY[0]-Bfi[0])/gauss/delta*mm, (BfiZ[0]-Bfi[0])/gauss/delta*mm,
(BfiX[1]-Bfi[1])/gauss/delta*mm, (BfiY[1]-Bfi[1])/gauss/delta*mm, (BfiZ[1]-Bfi[1])/gauss/delta*mm,
(BfiX[2]-Bfi[2])/gauss/delta*mm, (BfiY[2]-Bfi[2])/gauss/delta*mm, (BfiZ[2]-Bfi[2])/gauss/delta*mm );
}
} }

5
src/musrDetectorConstruction.cc
View File

@ -652,7 +652,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
} }
} }
else if (strcmp(tmpString2,"printFieldValueAtPoint")==0){ // Print the fieldvalue at the given point else if ((strcmp(tmpString2,"printFieldValueAtPoint")==0)||(strcmp(tmpString2,"printFieldDerivativeAtPoint")==0)) { // Print the fieldvalue at the given point
float p0, p1, p2; float p0, p1, p2;
sscanf(&line[0],"%*s %*s %*s %g %g %g",&p0,&p1,&p2); sscanf(&line[0],"%*s %*s %*s %g %g %g",&p0,&p1,&p2);
if (F04GlobalField::Exists()) { if (F04GlobalField::Exists()) {
@ -660,7 +660,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (myGlobalField) { if (myGlobalField) {
// The global field is not properly initialised at this point. Therefore the points have to // The global field is not properly initialised at this point. Therefore the points have to
// be stored in "F04GlobalField" object and printed later on in musrRunAction.cc . // be stored in "F04GlobalField" object and printed later on in musrRunAction.cc .
myGlobalField->AddPointForFieldTesting(G4ThreeVector(p0,p1,p2)); if (strcmp(tmpString2,"printFieldValueAtPoint")==0) myGlobalField->AddPointForFieldTesting(G4ThreeVector(p0,p1,p2));
if (strcmp(tmpString2,"printFieldDerivativeAtPoint")==0) myGlobalField->AddPointForFieldDerivativeTesting(G4ThreeVector(p0,p1,p2));
} }
else { else {
sprintf(eMessage,"musrDetectorConstruction.cc::Construct(): printFieldValueAtPoint requested, but field not found"); sprintf(eMessage,"musrDetectorConstruction.cc::Construct(): printFieldValueAtPoint requested, but field not found");

19
src/musrPhysicsList.cc
View File

@ -114,13 +114,15 @@ void musrPhysicsList::ConstructLeptons()
//cks: Trial to use Geant4 muon decay with spin //cks: Trial to use Geant4 muon decay with spin
G4DecayTable* MuonPlusDecayTable = new G4DecayTable(); G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",1.00)); // MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",1.00));
// MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",0.986)); MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",0.986));
// MuonPlusDecayTable -> Insert(new G4MuonRadiativeDecayChannelWithSpin("mu+",0.014)); MuonPlusDecayTable -> Insert(new G4MuonRadiativeDecayChannelWithSpin("mu+",0.014));
G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable); G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
// //
G4DecayTable* MuonMinusDecayTable = new G4DecayTable(); G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00)); // MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00));
MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",0.986));
MuonMinusDecayTable -> Insert(new G4MuonRadiativeDecayChannelWithSpin("mu-",0.014));
G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable); G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
//csk //csk
// //
@ -181,17 +183,21 @@ void musrPhysicsList::ConstructProcess()
#include "G4ComptonScattering.hh" #include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh" #include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh" #include "G4PhotoElectricEffect.hh"
#include "G4RayleighScattering.hh"
#include "G4MultipleScattering.hh" #include "G4MultipleScattering.hh"
#include "G4eMultipleScattering.hh"
#include "G4eIonisation.hh" #include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh" #include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh" #include "G4eplusAnnihilation.hh"
#include "G4MuMultipleScattering.hh"
#include "G4MuIonisation.hh" #include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh" #include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh" #include "G4MuPairProduction.hh"
#include "G4hMultipleScattering.hh"
#include "G4hIonisation.hh" #include "G4hIonisation.hh"
#include "G4UserSpecialCuts.hh" #include "G4UserSpecialCuts.hh"
@ -271,6 +277,7 @@ void musrPhysicsList::ConstructEM()
if (stringProcessName=="G4PhotoElectricEffect") pManager->AddDiscreteProcess(new G4PhotoElectricEffect); if (stringProcessName=="G4PhotoElectricEffect") pManager->AddDiscreteProcess(new G4PhotoElectricEffect);
else if (stringProcessName=="G4ComptonScattering") pManager->AddDiscreteProcess(new G4ComptonScattering); else if (stringProcessName=="G4ComptonScattering") pManager->AddDiscreteProcess(new G4ComptonScattering);
else if (stringProcessName=="G4GammaConversion") pManager->AddDiscreteProcess(new G4GammaConversion); else if (stringProcessName=="G4GammaConversion") pManager->AddDiscreteProcess(new G4GammaConversion);
else if (stringProcessName=="G4RayleighScattering") pManager->AddDiscreteProcess(new G4RayleighScattering);
else if (stringProcessName=="G4PenelopePhotoElectric") pManager->AddDiscreteProcess(new G4PenelopePhotoElectric); else if (stringProcessName=="G4PenelopePhotoElectric") pManager->AddDiscreteProcess(new G4PenelopePhotoElectric);
else if (stringProcessName=="G4PenelopeCompton") pManager->AddDiscreteProcess(new G4PenelopeCompton); else if (stringProcessName=="G4PenelopeCompton") pManager->AddDiscreteProcess(new G4PenelopeCompton);
else if (stringProcessName=="G4PenelopeGammaConversion") pManager->AddDiscreteProcess(new G4PenelopeGammaConversion); else if (stringProcessName=="G4PenelopeGammaConversion") pManager->AddDiscreteProcess(new G4PenelopeGammaConversion);
@ -292,6 +299,7 @@ void musrPhysicsList::ConstructEM()
char charRegion1[100]="", charRegion2[100]="", charRegion3[100]="", charControlString[10]=""; char charRegion1[100]="", charRegion2[100]="", charRegion3[100]="", charControlString[10]="";
sscanf(&line[0],"%*s %*s %*s %*s %*s %d %d %d %s %s %s %s",&nr1,&nr2,&nr3,charRegion1,charRegion2,charRegion3,charControlString); sscanf(&line[0],"%*s %*s %*s %*s %*s %d %d %d %s %s %s %s",&nr1,&nr2,&nr3,charRegion1,charRegion2,charRegion3,charControlString);
if (stringProcessName=="G4MultipleScattering") pManager->AddProcess(new G4MultipleScattering,nr1,nr2,nr3); if (stringProcessName=="G4MultipleScattering") pManager->AddProcess(new G4MultipleScattering,nr1,nr2,nr3);
else if (stringProcessName=="G4eMultipleScattering") pManager->AddProcess(new G4eMultipleScattering,nr1,nr2,nr3);
else if (stringProcessName=="G4eIonisation") pManager->AddProcess(new G4eIonisation,nr1,nr2,nr3); else if (stringProcessName=="G4eIonisation") pManager->AddProcess(new G4eIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4eBremsstrahlung") pManager->AddProcess(new G4eBremsstrahlung,nr1,nr2,nr3); else if (stringProcessName=="G4eBremsstrahlung") pManager->AddProcess(new G4eBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4eplusAnnihilation") pManager->AddProcess(new G4eplusAnnihilation,nr1,nr2,nr3); else if (stringProcessName=="G4eplusAnnihilation") pManager->AddProcess(new G4eplusAnnihilation,nr1,nr2,nr3);
@ -300,6 +308,7 @@ void musrPhysicsList::ConstructEM()
else if (stringProcessName=="G4PenelopeAnnihilation") pManager->AddProcess(new G4PenelopeAnnihilation,nr1,nr2,nr3); else if (stringProcessName=="G4PenelopeAnnihilation") pManager->AddProcess(new G4PenelopeAnnihilation,nr1,nr2,nr3);
else if (stringProcessName=="G4LowEnergyIonisation") pManager->AddProcess(new G4LowEnergyIonisation,nr1,nr2,nr3); else if (stringProcessName=="G4LowEnergyIonisation") pManager->AddProcess(new G4LowEnergyIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4LowEnergyBremsstrahlung") pManager->AddProcess(new G4LowEnergyBremsstrahlung,nr1,nr2,nr3); else if (stringProcessName=="G4LowEnergyBremsstrahlung") pManager->AddProcess(new G4LowEnergyBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuMultipleScattering") pManager->AddProcess(new G4MuMultipleScattering,nr1,nr2,nr3);
else if (stringProcessName=="G4MuIonisation") pManager->AddProcess(new G4MuIonisation,nr1,nr2,nr3); else if (stringProcessName=="G4MuIonisation") pManager->AddProcess(new G4MuIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4MuBremsstrahlung") pManager->AddProcess(new G4MuBremsstrahlung,nr1,nr2,nr3); else if (stringProcessName=="G4MuBremsstrahlung") pManager->AddProcess(new G4MuBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuPairProduction") pManager->AddProcess(new G4MuPairProduction,nr1,nr2,nr3); else if (stringProcessName=="G4MuPairProduction") pManager->AddProcess(new G4MuPairProduction,nr1,nr2,nr3);
@ -416,7 +425,7 @@ void musrPhysicsList::ConstructEM()
(particle->GetPDGCharge() != 0.0) && (particle->GetPDGCharge() != 0.0) &&
(particle->GetParticleName() != "chargedgeantino")) { (particle->GetParticleName() != "chargedgeantino")) {
//all others charged particles except geantino //all others charged particles except geantino
pmanager->AddProcess(new G4MultipleScattering,-1, 1,1); pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1);
// if (myTypeOfProcesses=="highenergy") { // if (myTypeOfProcesses=="highenergy") {
// pmanager->AddProcess(new G4hIonisation, -1, 2,2); // pmanager->AddProcess(new G4hIonisation, -1, 2,2);
// } // }

505
src/musrPhysicsList.cc_Geant4.9.2p02_and_older Normal file
View File

@ -0,0 +1,505 @@
/***************************************************************************
* musrSim - the program for the simulation of (mainly) muSR instruments. *
* More info on http://lmu.web.psi.ch/simulation/index.html . *
* musrSim is based od Geant4 (http://geant4.web.cern.ch/geant4/) *
* *
* Copyright (C) 2009 by Paul Scherrer Institut, 5232 Villigen PSI, *
* Switzerland *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the Free Software *
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *
***************************************************************************/
#include "globals.hh"
#include "G4ios.hh"
#include "musrPhysicsList.hh"
#include "G4VPhysicsConstructor.hh"
#include "G4ProcessManager.hh"
#include "G4ParticleTypes.hh"
#include "G4MuonDecayChannel.hh"
#include "G4DecayTable.hh"
//cks Added to have Geant default muon decay with spin
#include "G4MuonDecayChannelWithSpin.hh"
#include "G4MuonRadiativeDecayChannelWithSpin.hh"
#include "G4RadioactiveDecay.hh"
#include "G4IonConstructor.hh"
//TS Classes which account for Muonium as "particle" and its spin
#include "musrMuonium.hh"
#include "MuDecayChannel.hh"
#include "MuDecayChannelWithSpin.hh"
//
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
// cks 2009-06-08 G4StepLimiter and/or G4UserSpecialCuts are needed to activate the "G4UserLimits"
#include "G4StepLimiter.hh"
#include "G4UserSpecialCuts.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrPhysicsList::musrPhysicsList(): G4VUserPhysicsList()
{
defaultCutValue = 0.1*mm;
cutForGamma = 0.1*mm;
cutForElectron = 0.1*mm;
cutForMuon = 0.01*mm;
SetVerboseLevel(0);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrPhysicsList::~musrPhysicsList()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructParticle()
{
ConstructBosons();
ConstructLeptons();
ConstructMesons();
ConstructBaryons();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructBosons()
{
// pseudo-particles
G4Geantino::GeantinoDefinition();
G4ChargedGeantino::ChargedGeantinoDefinition();
// gamma
G4Gamma::GammaDefinition();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructLeptons()
{
// leptons
// e+/-
G4Electron::ElectronDefinition();
G4Positron::PositronDefinition();
// mu+/-
G4MuonPlus::MuonPlusDefinition();
G4MuonMinus::MuonMinusDefinition();
//cks
// G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
// MuonPlusDecayTable -> Insert(new musrMuonDecayChannel("mu+",1.00));
// G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
//csk
//
// Muonium - TS
musrMuonium::MuoniumDefinition();
//
// nu_e
G4NeutrinoE::NeutrinoEDefinition();
G4AntiNeutrinoE::AntiNeutrinoEDefinition();
// nu_mu
G4NeutrinoMu::NeutrinoMuDefinition();
G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
//cks: Trial to use Geant4 muon decay with spin
G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",1.00));
// MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",0.986));
// MuonPlusDecayTable -> Insert(new G4MuonRadiativeDecayChannelWithSpin("mu+",0.014));
G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
//
G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00));
G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
//csk
//
//TS: Using the muonium decay with and without spin
G4DecayTable* MuoniumDecayTable = new G4DecayTable();
MuoniumDecayTable -> Insert(new MuDecayChannel("Mu",0.50));
MuoniumDecayTable -> Insert(new MuDecayChannelWithSpin("Mu",0.5));
musrMuonium::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable);
//MuoniumDecayTable ->DumpInfo(); // Info on muonium decay channels
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructMesons()
{
// mesons
// light mesons
G4PionPlus::PionPlusDefinition();
G4PionMinus::PionMinusDefinition();
G4PionZero::PionZeroDefinition();
G4Eta::EtaDefinition();
G4EtaPrime::EtaPrimeDefinition();
G4KaonPlus::KaonPlusDefinition();
G4KaonMinus::KaonMinusDefinition();
G4KaonZero::KaonZeroDefinition();
G4AntiKaonZero::AntiKaonZeroDefinition();
G4KaonZeroLong::KaonZeroLongDefinition();
G4KaonZeroShort::KaonZeroShortDefinition();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructBaryons()
{
// baryons
G4Proton::ProtonDefinition();
G4AntiProton::AntiProtonDefinition();
G4Neutron::NeutronDefinition();
G4AntiNeutron::AntiNeutronDefinition();
// ions
G4IonConstructor iConstructor;
iConstructor.ConstructParticle();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructProcess()
{
AddTransportation();
ConstructEM();
ConstructGeneral();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"
#include "G4MultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4hIonisation.hh"
#include "G4UserSpecialCuts.hh"
//#include "musrAtRestSpinRotation.hh"
// For low energy physics processes:
#include "G4LowEnergyCompton.hh"
//#include "G4LowEnergyPolarizedCompton.hh"
#include "G4LowEnergyGammaConversion.hh"
#include "G4LowEnergyPhotoElectric.hh"
#include "G4LowEnergyRayleigh.hh"
#include "G4LowEnergyBremsstrahlung.hh"
#include "G4LowEnergyIonisation.hh"
#include "G4hLowEnergyIonisation.hh"
// For Penelope processes:
#include "G4PenelopeCompton.hh"
#include "G4PenelopeGammaConversion.hh"
#include "G4PenelopePhotoElectric.hh"
#include "G4PenelopeRayleigh.hh"
#include "G4PenelopeIonisation.hh"
#include "G4PenelopeBremsstrahlung.hh"
#include "G4PenelopeAnnihilation.hh"
// For Coulomb scattering instead of multiple scattering
#include "G4CoulombScattering.hh"
#include "G4CoulombScatteringModel.hh"
// For Muonium formation in the Carbon foil
#include "musrMuFormation.hh" // includes the yield function Y = Y(E).
// For a simple Muonium "scattering" when Mu hits solid materials
#include "musrMuScatter.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructEM()
{
// cks 2008.08.22. - Adding the possibility to define the processes from the steering file:
char charSteeringFileName[1000]; strcpy(charSteeringFileName,(musrParameters::mySteeringFileName).c_str());
FILE *fSteeringFile=fopen(charSteeringFileName,"r");
if (fSteeringFile==NULL) {
sprintf(eMessage,"musrPhysicsList::ConstructEM(): Failed to open macro file \"%s\" .",charSteeringFileName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
else {G4cout<<"musrPhysicsList: The Physics processes are being defined:"<<G4endl;}
char line[501];
while (!feof(fSteeringFile)) {
fgets(line,500,fSteeringFile);
if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) {
char tmpString0[100]="Unset", tmpString1[100]="Unset",tmpString2[100]="Unset";
sscanf(&line[0],"%s %s %s",tmpString0,tmpString1,tmpString2);
if ( (strcmp(tmpString0,"/musr/ignore")!=0)&&(strcmp(tmpString0,"/musr/command")!=0) ) continue;
if (strcmp(tmpString1,"process")!=0) continue;
float tmpCutValue;
if ((strcmp(tmpString2,"addProcess")==0)||(strcmp(tmpString2,"addDiscreteProcess")==0)) {
char charParticleName[100], charProcessName[100];
sscanf(&line[0],"%*s %*s %s %s %s",tmpString2,charParticleName,charProcessName);
G4cout<<"musrPhysicsList: Defining process "<<charProcessName<<" for "<<charParticleName<<G4endl;
G4String stringProcessName = charProcessName;
G4String stringParticleName = charParticleName;
G4ParticleDefinition* particleDefinition = G4ParticleTable::GetParticleTable() -> FindParticle(stringParticleName);
// G4cout<<"particleDefinition of "<<stringParticleName<<" = "<<particleDefinition<<G4endl;
if (particleDefinition==NULL) {
sprintf(eMessage,"musrPhysicsList: Partile \"%s\" not found in G4ParticleTable when trying to assign process \"%s\".",
charParticleName,charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
G4ProcessManager* pManager = particleDefinition->GetProcessManager();
if (strcmp(tmpString2,"addDiscreteProcess")==0) {
if (stringProcessName=="G4PhotoElectricEffect") pManager->AddDiscreteProcess(new G4PhotoElectricEffect);
else if (stringProcessName=="G4ComptonScattering") pManager->AddDiscreteProcess(new G4ComptonScattering);
else if (stringProcessName=="G4GammaConversion") pManager->AddDiscreteProcess(new G4GammaConversion);
else if (stringProcessName=="G4PenelopePhotoElectric") pManager->AddDiscreteProcess(new G4PenelopePhotoElectric);
else if (stringProcessName=="G4PenelopeCompton") pManager->AddDiscreteProcess(new G4PenelopeCompton);
else if (stringProcessName=="G4PenelopeGammaConversion") pManager->AddDiscreteProcess(new G4PenelopeGammaConversion);
else if (stringProcessName=="G4PenelopeRayleigh") pManager->AddDiscreteProcess(new G4PenelopeRayleigh);
else if (stringProcessName=="G4LowEnergyPhotoElectric") pManager->AddDiscreteProcess(new G4LowEnergyPhotoElectric);
else if (stringProcessName=="G4LowEnergyCompton") pManager->AddDiscreteProcess(new G4LowEnergyCompton);
else if (stringProcessName=="G4LowEnergyGammaConversion") pManager->AddDiscreteProcess(new G4LowEnergyGammaConversion);
else if (stringProcessName=="G4LowEnergyRayleigh") pManager->AddDiscreteProcess(new G4LowEnergyRayleigh);
else if (stringProcessName=="G4CoulombScattering") pManager->AddDiscreteProcess(new G4CoulombScattering);
else {
sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess. It can be easily added.",
charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
}
else if (strcmp(tmpString2,"addProcess")==0) {
G4int nr1, nr2, nr3;
char charRegion1[100]="", charRegion2[100]="", charRegion3[100]="", charControlString[10]="";
sscanf(&line[0],"%*s %*s %*s %*s %*s %d %d %d %s %s %s %s",&nr1,&nr2,&nr3,charRegion1,charRegion2,charRegion3,charControlString);
if (stringProcessName=="G4MultipleScattering") pManager->AddProcess(new G4MultipleScattering,nr1,nr2,nr3);
else if (stringProcessName=="G4eIonisation") pManager->AddProcess(new G4eIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4eBremsstrahlung") pManager->AddProcess(new G4eBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4eplusAnnihilation") pManager->AddProcess(new G4eplusAnnihilation,nr1,nr2,nr3);
else if (stringProcessName=="G4PenelopeIonisation") pManager->AddProcess(new G4PenelopeIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4PenelopeBremsstrahlung") pManager->AddProcess(new G4PenelopeBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4PenelopeAnnihilation") pManager->AddProcess(new G4PenelopeAnnihilation,nr1,nr2,nr3);
else if (stringProcessName=="G4LowEnergyIonisation") pManager->AddProcess(new G4LowEnergyIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4LowEnergyBremsstrahlung") pManager->AddProcess(new G4LowEnergyBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuIonisation") pManager->AddProcess(new G4MuIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4MuBremsstrahlung") pManager->AddProcess(new G4MuBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuPairProduction") pManager->AddProcess(new G4MuPairProduction,nr1,nr2,nr3);
// cks 2009-06-08 G4StepLimiter and/or G4UserSpecialCuts are needed to activate the "G4UserLimits"
else if (stringProcessName=="G4StepLimiter") pManager->AddProcess(new G4StepLimiter,nr1,nr2,nr3);
else if (stringProcessName=="G4UserSpecialCuts") pManager->AddProcess(new G4UserSpecialCuts,nr1,nr2,nr3);
// else if (stringProcessName=="G4DecayWithSpin") pManager->AddProcess(new G4DecayWithSpin,nr1,nr2,nr3);
// else if (stringProcessName=="G4hIonisation") pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
// else if (stringProcessName=="G4hLowEnergyIonisation") pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
else if (stringProcessName=="musrMuFormation") pManager->AddProcess(new musrMuFormation,nr1,nr2,nr3);
// cks: musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments
// in his original "musrPhysicsList.cc about implementing musrMuScatter.
// else if (stringProcessName=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3);
else if (stringProcessName=="MultipleAndCoulombScattering") {
G4MultipleScattering* multScat = new G4MultipleScattering();
// G4CoulombScattering* coulScat = new G4CoulombScattering();
G4CoulombScatteringModel* coulScatModel = new G4CoulombScatteringModel();
if (strcmp(charRegion1,"")!=0) {
G4Region* regionForCoulomb = FindG4Region(charRegion1,line);
G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<<charRegion1<<G4endl;
multScat->AddEmModel(0,coulScatModel,regionForCoulomb);
// multScat->AddEmModel(0,multScat,regionForCoulomb);
}
if (strcmp(charRegion2,"")!=0) {
G4Region* regionForCoulomb = FindG4Region(charRegion2,line);
G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<<charRegion2<<G4endl;
multScat->AddEmModel(0,coulScatModel,regionForCoulomb);
}
if (strcmp(charRegion3,"")!=0) {
G4Region* regionForCoulomb = FindG4Region(charRegion3,line);
G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<<charRegion3<<G4endl;
multScat->AddEmModel(0,coulScatModel,regionForCoulomb);
}
if (strcmp(charControlString,"")!=0) {
G4cout<<"More than 3 regions requested for Coulomb Scattering, but presently only up to 3 such regions are supported."<<G4endl;
G4cout<<"Please extend the number of supported regions in musrPhysicsList.cc to higher number."<<G4endl;
G4cout<<"The extention of the code to larger number of regions is not very difficult."<<G4endl;
G4cout<<" S T O P F O R C E D"<<G4endl;
exit(1);
}
pManager->AddProcess(multScat,nr1,nr2,nr3);
}
else {
sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addProcess. It can be easily added.",
charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
}
}
else if (strcmp(tmpString2,"cutForGamma")==0) {
sscanf(&line[0],"%*s %*s %*s %g",&tmpCutValue);
cutForGamma = tmpCutValue;
}
else if (strcmp(tmpString2,"cutForElectron")==0) {
sscanf(&line[0],"%*s %*s %*s %g",&tmpCutValue);
cutForElectron = tmpCutValue;
}
//cks UNFORTUNATELY cutForMuon does not work!!!
// else if (strcmp(tmpString2,"cutForMuon")==0) {
// sscanf(&line[0],"%*s %*s %*s %g",&tmpCutValue);
// cutForMuon = tmpCutValue;
// }
else ReportProblemWithProcessDefinition(line);
}
}
fclose(fSteeringFile);
G4cout<<"\n\n\n\n"<<G4endl;
// csk 2008.08.22.
//del G4String myTypeOfProcesses = musrParameters::GetInstance()->GetMyTypeOfProcesses();
//del G4cout<<"musrPhysicsList::ConstructEM(): myTypeOfProcesses="<<myTypeOfProcesses<<G4endl;
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if ((particleName == "gamma")||(particleName == "e-")||(particleName == "e+")) {
// do nothing
}
else if ((particleName=="mu+")||(particleName=="mu-")) { //muon
G4DecayWithSpin* theDecayProcess = new G4DecayWithSpin();
// theDecayProcess->SetVerboseLevel(2);
pmanager->AddProcess(theDecayProcess);
pmanager ->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
}
else if (particleName=="Mu") {
// TS:
// Muonium "scattering" Kamil: the following 3 lines could be replaced by reading the musrMuScatter
// process through the steering file
G4VProcess* aMuScatt = new musrMuScatter();
pmanager->AddProcess(aMuScatt);
pmanager->SetProcessOrdering(aMuScatt, idxPostStep, 1);
//
G4Decay* theDecayProcess = new G4Decay();
//musrDecayWithSpin* theDecayProcess = new musrDecayWithSpin();
pmanager->AddProcess(theDecayProcess);
pmanager->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
}
else if ((!particle->IsShortLived()) &&
(particle->GetPDGCharge() != 0.0) &&
(particle->GetParticleName() != "chargedgeantino")) {
//all others charged particles except geantino
pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
// if (myTypeOfProcesses=="highenergy") {
// pmanager->AddProcess(new G4hIonisation, -1, 2,2);
// }
// else {
pmanager->AddProcess(new G4hLowEnergyIonisation, -1, 2,2);
// }
///pmanager->AddProcess(new G4UserSpecialCuts, -1,-1,3);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4Decay.hh"
void musrPhysicsList::ConstructGeneral() {
if (musrParameters::boolG4GeneralParticleSource) {
G4RadioactiveDecay* theRadioactiveDecay = new G4RadioactiveDecay();
G4GenericIon* ion = G4GenericIon::GenericIon();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (particle == ion) {
pmanager->AddProcess(theRadioactiveDecay, 0, -1, 3);
}
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4Region.hh"
#include "G4RegionStore.hh"
#include "G4ProductionCuts.hh"
void musrPhysicsList::SetCuts()
{
//G4VUserPhysicsList::SetCutsWithDefault method sets
//the default cut value for all particle types
//
//cks 6.10.2009 SetCutsWithDefault();
// set cut values for gamma at first and for e- second and next for e+,
// because some processes for e+/e- need cut values for gamma
SetCutValue(cutForGamma, "gamma");
SetCutValue(cutForElectron, "e-");
SetCutValue(cutForElectron, "e+");
//cks - This command does not work for muons: SetCutValue(cutForMuon, "mu-");
//cks - This command does not work for muons: SetCutValue(cutForMuon, "mu+");
// G4cout<<"Kamil: cutForGamma = "<<cutForGamma/mm<<" mm"<<G4endl;
// G4cout<<"Kamil: cutForElectron = "<<cutForElectron/mm<<" mm"<<G4endl;
// G4cout<<"Kamil: cutForMuons = "<<cutForMuon/mm<<" mm"<<G4endl;
if (verboseLevel>0) DumpCutValuesTable();
DumpCutValuesTable();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ReportProblemWithProcessDefinition(char myString[501]) {
G4cout<<"\nE R R O R in musrPhysicsList.cc: "
<<"Unknown keyword requested in the steering (*.mac) file :"<<G4endl;
G4cout<<" "<<myString<<G4endl;
G4cout<<"S T O P F O R C E D!"<<G4endl;
exit(1);
}
G4Region* musrPhysicsList::FindG4Region(G4String regionName, char* lineOfSteeringFile) {
G4Region* myRegion = G4RegionStore::GetInstance()->GetRegion(regionName,false);
if( myRegion != NULL ) { // G4Region found
return myRegion;
}
else { // G4Region not found
G4cout<<"musrPhysicsList: G4Region "<<regionName<<" not found."<<G4endl;
G4cout<<" The critical command line of the steering file is:"<<G4endl;
G4cout<<" "<<lineOfSteeringFile<<G4endl;
G4cout<<" S T O P F O R C E D"<<G4endl;
exit(1);
}
}

2
src/musrRootOutput.cc
View File

@ -47,6 +47,8 @@ musrRootOutput::musrRootOutput() {
ProcessIDMapping["muIoni"]=11; ProcessIDMapping["muIoni"]=11;
ProcessIDMapping["MuFormation"]=12; ProcessIDMapping["MuFormation"]=12;
ProcessIDMapping["Decay"]=13; ProcessIDMapping["Decay"]=13;
ProcessIDMapping["conv"]=14;
ProcessIDMapping["compt"]=15;
ProcessIDMapping["initialParticle"]=100; ProcessIDMapping["initialParticle"]=100;
} }