diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index d9b5be4..e40e799 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index a0d73f8..3cd72f2 100644 --- a/doc/musrSim.tex +++ b/doc/musrSim.tex @@ -426,7 +426,11 @@ Three special volumes ``Target, M0, M1 and M2''. \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 - 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} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -950,6 +954,25 @@ These are typically not interesting for the analysis of the results of the simul \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)} 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. diff --git a/include/F04GlobalField.hh b/include/F04GlobalField.hh index 65501e3..c0d4704 100644 --- a/include/F04GlobalField.hh +++ b/include/F04GlobalField.hh @@ -141,9 +141,9 @@ public: // G4bool DoesAnyFieldValueNeedsToBeChanged(G4int eventNumber) {return globalChangeFieldInStepsMap[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 AddPointForFieldDerivativeTesting(G4ThreeVector point) {pointsAtWhichUserWantsToPrintFieldDerivative.push_back(point);} // Print field value at all points user wished to be print out: void PrintFieldAtRequestedPoints() const; @@ -164,6 +164,7 @@ private: const F04ElementField **fp; std::vector pointsAtWhichUserWantsToPrintFieldValue; + std::vector pointsAtWhichUserWantsToPrintFieldDerivative; private: diff --git a/src/F04GlobalField.cc b/src/F04GlobalField.cc index f0c024f..5f0720f 100644 --- a/src/F04GlobalField.cc +++ b/src/F04GlobalField.cc @@ -308,9 +308,14 @@ void F04GlobalField::CheckWhetherAnyNominalFieldValueNeedsToBeChanged(G4int even // Print field value at all points requested by the user: void F04GlobalField::PrintFieldAtRequestedPoints() const { - G4ThreeVector p; + G4ThreeVector p; 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 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++) { p = pointsAtWhichUserWantsToPrintFieldValue[i]; point[0] = p.x(); @@ -325,6 +330,34 @@ void F04GlobalField::PrintFieldAtRequestedPoints() const { Bfi[0]/tesla, Bfi[1]/tesla, Bfi[2]/tesla, Bfi[3]/(kilovolt/mm), Bfi[4]/(kilovolt/mm), Bfi[5]/(kilovolt/mm)); } + + if (pointsAtWhichUserWantsToPrintFieldDerivative.size()>0) { + G4cout< tesla = "<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 ); + } } diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 1b1688d..d4bc3f2 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -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; sscanf(&line[0],"%*s %*s %*s %g %g %g",&p0,&p1,&p2); if (F04GlobalField::Exists()) { @@ -660,7 +660,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { if (myGlobalField) { // 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 . - 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 { sprintf(eMessage,"musrDetectorConstruction.cc::Construct(): printFieldValueAtPoint requested, but field not found"); diff --git a/src/musrPhysicsList.cc b/src/musrPhysicsList.cc index f78845a..d634e4b 100644 --- a/src/musrPhysicsList.cc +++ b/src/musrPhysicsList.cc @@ -114,13 +114,15 @@ void musrPhysicsList::ConstructLeptons() //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)); + // 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)); + // 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); //csk // @@ -181,17 +183,21 @@ void musrPhysicsList::ConstructProcess() #include "G4ComptonScattering.hh" #include "G4GammaConversion.hh" #include "G4PhotoElectricEffect.hh" +#include "G4RayleighScattering.hh" #include "G4MultipleScattering.hh" +#include "G4eMultipleScattering.hh" #include "G4eIonisation.hh" #include "G4eBremsstrahlung.hh" #include "G4eplusAnnihilation.hh" +#include "G4MuMultipleScattering.hh" #include "G4MuIonisation.hh" #include "G4MuBremsstrahlung.hh" #include "G4MuPairProduction.hh" +#include "G4hMultipleScattering.hh" #include "G4hIonisation.hh" #include "G4UserSpecialCuts.hh" @@ -271,6 +277,7 @@ void musrPhysicsList::ConstructEM() 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=="G4RayleighScattering") pManager->AddDiscreteProcess(new G4RayleighScattering); else if (stringProcessName=="G4PenelopePhotoElectric") pManager->AddDiscreteProcess(new G4PenelopePhotoElectric); else if (stringProcessName=="G4PenelopeCompton") pManager->AddDiscreteProcess(new G4PenelopeCompton); else if (stringProcessName=="G4PenelopeGammaConversion") pManager->AddDiscreteProcess(new G4PenelopeGammaConversion); @@ -291,28 +298,30 @@ void musrPhysicsList::ConstructEM() 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); + 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=="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=="G4MuMultipleScattering") pManager->AddProcess(new G4MuMultipleScattering,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=="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=="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); + 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=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3); else if (stringProcessName=="MultipleAndCoulombScattering") { G4MultipleScattering* multScat = new G4MultipleScattering(); // G4CoulombScattering* coulScat = new G4CoulombScattering(); @@ -416,7 +425,7 @@ void musrPhysicsList::ConstructEM() (particle->GetPDGCharge() != 0.0) && (particle->GetParticleName() != "chargedgeantino")) { //all others charged particles except geantino - pmanager->AddProcess(new G4MultipleScattering,-1, 1,1); + pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1); // if (myTypeOfProcesses=="highenergy") { // pmanager->AddProcess(new G4hIonisation, -1, 2,2); // } diff --git a/src/musrPhysicsList.cc_Geant4.9.2p02_and_older b/src/musrPhysicsList.cc_Geant4.9.2p02_and_older new file mode 100644 index 0000000..f78845a --- /dev/null +++ b/src/musrPhysicsList.cc_Geant4.9.2p02_and_older @@ -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:"< FindParticle(stringParticleName); + // G4cout<<"particleDefinition of "<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 "<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 "<AddEmModel(0,coulScatModel,regionForCoulomb); + } + if (strcmp(charRegion3,"")!=0) { + G4Region* regionForCoulomb = FindG4Region(charRegion3,line); + G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<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."<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"<GetMyTypeOfProcesses(); + //del G4cout<<"musrPhysicsList::ConstructEM(): myTypeOfProcesses="<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 = "<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 :"<GetRegion(regionName,false); + if( myRegion != NULL ) { // G4Region found + return myRegion; + } + else { // G4Region not found + G4cout<<"musrPhysicsList: G4Region "<