diff --git a/include/musrDetectorConstruction.hh b/include/musrDetectorConstruction.hh index 214870a..b252e19 100644 --- a/include/musrDetectorConstruction.hh +++ b/include/musrDetectorConstruction.hh @@ -41,6 +41,7 @@ class G4VPhysicalVolume; class G4Material; class musrDetectorMessenger; class musrScintSD; +class musrMuEnergyLossLandau; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/include/musrMuEnergyLossLandau.hh b/include/musrMuEnergyLossLandau.hh new file mode 100644 index 0000000..95152db --- /dev/null +++ b/include/musrMuEnergyLossLandau.hh @@ -0,0 +1,95 @@ +/*************************************************************************** + * 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. * + ***************************************************************************/ + +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +// Muonium Formation according to yield.cc function (through GetYields method). +// Id : musrMuFormation.hh, v 1.4 +// Author: Taofiq PARAISO, T. Shiroka +// Date : 2007-12 +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +#ifndef musrMuEnergyLossLandau_h +#define musrMuEnergyLossLandau_h 1 + +#include "G4VDiscreteProcess.hh" +#include "G4ParticleTable.hh" + +#include +#include + + +/*! musrMuEnergyLossLandau class defines the energy loss given by the Landau distribution, + in, for example, the thin Carbon foil of the LEM beam line + */ + +class musrMuEnergyLossLandau : public G4VDiscreteProcess +{ + public: + + musrMuEnergyLossLandau(const G4String& name = "MuEnergyLossLandau", // process description + G4ProcessType aType = fElectromagnetic); + + ~musrMuEnergyLossLandau(); + + //! - Main method. muon energy loss process is executed at the END of a step. */ + G4VParticleChange* PostStepDoIt( + const G4Track&, + const G4Step&); + + G4double GetMeanFreePath(const G4Track& aTrack, + G4double previousStepSize, + G4ForceCondition* condition); + + + //! Condition for process application (step Object). + G4bool CheckCondition(const G4Step& aStep); + + //! Condition for process application (step Pointer). + G4bool CheckCondition(const G4Step* aStep); + + + G4String p_name; + G4bool condition; + + + void GetFinalEnergy( const G4Step* aStep); + TRandom *random; + static double landauMPV; + static double landauSigma; + + // model parameters + G4ParticleTable* particleTable; + G4ParticleDefinition* particle; + G4double rnd; + G4DynamicParticle *DP; + + //! The particle change object. + G4VParticleChange fParticleChange; + + void PrepareSecondary(const G4Track&); + G4Track* aSecondary; + + void InitializeSecondaries(const G4Track&); +}; + +#endif diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 828719b..309404a 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -29,6 +29,7 @@ #include "musrScintSD.hh" #include "musrEventAction.hh" +#include "musrMuEnergyLossLandau.hh" #include "G4GeometryManager.hh" #include "G4PhysicalVolumeStore.hh" @@ -74,6 +75,8 @@ #include "musrErrorMessage.hh" #include "musrSteppingAction.hh" +using namespace std; + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... musrDetectorConstruction::musrDetectorConstruction(G4String steeringFileName) @@ -1537,7 +1540,15 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { ||(strcmp(tmpString1,"materialPropertiesTable")==0)||(strcmp(tmpString1,"setMaterialPropertiesTable")==0)) { ; // processes are interpreded later in musrPhysicsList.cc } - + + else if (strcmp(tmpString1,"SetLandauMPV")==0){ + sscanf(&line[0],"%*s %*s %lf",&musrMuEnergyLossLandau::landauMPV); + cout << "landauMPV = " << musrMuEnergyLossLandau::landauMPV << endl; + } + else if (strcmp(tmpString1,"SetLandauSigma")==0){ + sscanf(&line[0],"%*s %*s %lf",&musrMuEnergyLossLandau::landauSigma); + cout << "landauSigma = " << musrMuEnergyLossLandau::landauSigma << endl; + } else ReportGeometryProblem(line); diff --git a/src/musrMuEnergyLossLandau.cc b/src/musrMuEnergyLossLandau.cc new file mode 100644 index 0000000..b79c4d5 --- /dev/null +++ b/src/musrMuEnergyLossLandau.cc @@ -0,0 +1,132 @@ +/*************************************************************************** + * 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. * + ***************************************************************************/ + +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +// Muonium Formation according to yield.cc function (through GetYields method). +// Id : musrMuFormation.cc, v 1.4 +// Author: Taofiq PARAISO, T. Shiroka +// Date : 2007-12 +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +#include "musrMuEnergyLossLandau.hh" + +using namespace std; + +musrMuEnergyLossLandau::musrMuEnergyLossLandau(const G4String& name, G4ProcessType aType) + : G4VDiscreteProcess(name, aType) +{ + random = new TRandom(); +} + +musrMuEnergyLossLandau::~musrMuEnergyLossLandau(){} + +double musrMuEnergyLossLandau::landauMPV = 0.1; +double musrMuEnergyLossLandau::landauSigma = 0.1; + + +G4VParticleChange* musrMuEnergyLossLandau::PostStepDoIt(const G4Track& trackData, + const G4Step& aStep) +{ // Initialize ParticleChange (by setting all its members equal to + // the corresponding members in G4Track) + fParticleChange.Initialize(trackData); + + G4Track theNewTrack; + if(CheckCondition(aStep)) + { + GetFinalEnergy(&aStep); + G4Step theStep; + PrepareSecondary( trackData); + fParticleChange.AddSecondary(aSecondary); + fParticleChange.ProposeTrackStatus(fStopAndKill) ; + } + else + { + fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ; + } + return &fParticleChange; +} + + +G4bool musrMuEnergyLossLandau::CheckCondition(const G4Step& aStep) +{ // Decide when to call the MuEnergyLossLandau - i.e. for muons going through the C foil. + condition=false; + p_name = aStep.GetTrack()->GetDefinition()->GetParticleName(); // particle name + std::string logVolName = aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName(); + if(p_name == "mu+" && ((logVolName=="log_coulombCFoil")||(logVolName=="log_CFoil"))) + { + condition=true; + } + return condition; +} + + +G4double musrMuEnergyLossLandau::GetMeanFreePath(const G4Track&, + G4double, + G4ForceCondition* condition1) +{ + *condition1 = Forced; + return DBL_MAX; +} + + +void musrMuEnergyLossLandau::GetFinalEnergy(const G4Step* aStep) +{ // Particle generation + particleTable=G4ParticleTable::GetParticleTable(); + G4double Eloss, Efinal; + G4double E = aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy()/CLHEP::keV; + + // Positive muon + if(p_name=="mu+"){ + do{ + Eloss = random->Landau(landauMPV,landauSigma); + } while (Eloss < 0.); + + Efinal = E - Eloss; +// cout << "E, Eloss, Efinal = " << E << ", " << Eloss << ", " << Efinal << endl; + + particle = particleTable->FindParticle(p_name) ; + + // Set the new dynamic particle DP + DP = new G4DynamicParticle(particle, + aStep->GetTrack()->GetDynamicParticle()->GetMomentumDirection(), + Efinal/1000.); + +/* IMPORTANT : COPY THOSE DATA TO GET THE SAME PARTICLE PROPERTIES!!! + SHOULD BE KEPT WHEN BUILDING A PARTICLE CHANGE */ + DP->SetProperTime( aStep->GetTrack()->GetDynamicParticle()->GetProperTime()); + DP->SetPolarization(aStep->GetTrack()->GetDynamicParticle()->GetPolarization().x(), + aStep->GetTrack()->GetDynamicParticle()->GetPolarization().y(), + aStep->GetTrack()->GetDynamicParticle()->GetPolarization().z()); + DP->SetPreAssignedDecayProperTime( + aStep->GetTrack()->GetDynamicParticle()->GetPreAssignedDecayProperTime()); + } +} + + +void musrMuEnergyLossLandau::PrepareSecondary(const G4Track& track) +{ + if(p_name=="mu+") + { + aSecondary = new G4Track(DP,track.GetGlobalTime(),track.GetPosition()); + } +} diff --git a/src/musrPhysicsList.cc b/src/musrPhysicsList.cc index 2486b2a..4e695cc 100644 --- a/src/musrPhysicsList.cc +++ b/src/musrPhysicsList.cc @@ -263,6 +263,9 @@ void musrPhysicsList::ConstructProcess() // For Muonium formation in the Carbon foil #include "musrMuFormation.hh" // includes the yield function Y = Y(E). +// For muon energy loss in thin Carbon foil +#include "musrMuEnergyLossLandau.hh" + // For a simple Muonium "scattering" when Mu hits solid materials #include "musrMuScatter.hh" @@ -381,6 +384,7 @@ void musrPhysicsList::ConstructEM() // 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=="musrMuEnergyLossLandau") pManager->AddProcess(new musrMuEnergyLossLandau,nr1,nr2,nr3); else if (stringProcessName=="G4Scintillation") { G4Scintillation* myScintillation = new G4Scintillation(); // if (musrParameters::finiteRiseTimeInScintillator) myScintillation->SetFiniteRiseTime(true);