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

Browse Source
This commit is contained in:
shiroka
2008-06-16 14:05:29 +00:00
parent b8cb746e20
commit e09fc2f057
38 changed files with 2691 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

181
geant4/LEMuSR/include/F04ElementField.hh Normal file
View File

@ -0,0 +1,181 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * 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. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
//
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
#ifndef F04ElementField_h
#define F04ElementField_h 1
#include "globals.hh"
#include "G4Navigator.hh"
#include "G4TransportationManager.hh"
#include "G4UserLimits.hh"
#include "G4VisAttributes.hh"
// class F04ElementField - interface for the EM field of one element
// This is the interface class used by GlobalField to compute the field
// value at a given point[].
// An element that represents an element with an EM field will
// derive a class from this one and implement the computation for the
// element. The construct() function will add the derived object into
// GlobalField.
class F04ElementField
{
private:
F04ElementField& operator=(const F04ElementField&);
public:
/// Constructor.
F04ElementField(const G4ThreeVector, G4LogicalVolume*);
/// the actual implementation constructs the F04ElementField
void construct();
/// Destructor.
virtual ~F04ElementField() { if (aNavigator) delete aNavigator; }
/// setMaxStep(G4double) sets the max. step size
void setMaxStep(G4double s)
{
maxStep = s;
userLimits->SetMaxAllowedStep(maxStep);
lvolume->SetUserLimits(userLimits);
}
/// getMaxStep() returns the max. step size
G4double getMaxStep() { return maxStep; }
/// setColor(G4String) sets the color
void setColor(G4String c)
{
color = c;
lvolume->SetVisAttributes(getVisAttribute(color));
}
/// getColor() returns the color
G4String getColor() { return color; }
/// getVisAttribute() returns the appropriate G4VisAttributes.
static G4VisAttributes* getVisAttribute(G4String color);
/// setGlobalPoint() ensures that the point is within the global
/// bounding box of this ElementField's global coordinates.
/// Normally called 8 times for the corners of the local bounding
/// box, after a local->global coordinate transform.
/// If never called, the global bounding box is infinite.
/// BEWARE: if called only once, the bounding box is just a point.
void setGlobalPoint(const G4double point[4])
{
if(minX == -DBL_MAX || minX > point[0]) minX = point[0];
if(minY == -DBL_MAX || minY > point[1]) minY = point[1];
if(minZ == -DBL_MAX || minZ > point[2]) minZ = point[2];
if(maxX == DBL_MAX || maxX < point[0]) maxX = point[0];
if(maxY == DBL_MAX || maxY < point[1]) maxY = point[1];
if(maxZ == DBL_MAX || maxZ < point[2]) maxZ = point[2];
}
/// isInBoundingBox() returns true if the point is within the
/// global bounding box - global coordinates.
bool isInBoundingBox(const G4double point[4]) const
{
if(point[2] < minZ || point[2] > maxZ) return false;
if(point[0] < minX || point[0] > maxX) return false;
if(point[1] < minY || point[1] > maxY) return false;
return true;
}
/// addFieldValue() will add the field value for this element to field[].
/// Implementations must be sure to verify that point[] is within
/// the field region, and do nothing if not.
/// point[] is in global coordinates and geant4 units; x,y,z,t.
/// field[] is in geant4 units; Bx,By,Bz,Ex,Ey,Ez.
/// For efficiency, the caller may (but need not) call
/// isInBoundingBox(point), and only call this function if that
/// returns true.
virtual void
addFieldValue(const G4double point[4], G4double field[6]) const = 0;
virtual G4double getLength() = 0;
virtual G4double getWidth() = 0;
virtual G4double getHeight() = 0;
///ADDED BY KAMIL. TS!!
void SetElementFieldName(G4String name) {elementFieldName=name;}
G4String GetElementFieldName() {return elementFieldName;}
void SetEventNrDependentField(G4double initialField, G4double finalField, G4int nrOfSteps);
std::map<G4int,G4double> GetEventNrDependentField() const {return changeFieldInStepsMap;}
// void SetElementFieldValueIfNeeded(G4int eventNr) const {
// std::map<G4int,G4double>::iterator itr;
// if ( (itr = changeFieldInStepsMap.find(eventNr)) != changeFieldInStepsMap.end() ) {
// G4double newNominalFieldValue = itr->second ;
//
// G4cout<<"Nominal Field changed for "<<elementFieldName<<" to value "<<newNominalFieldValue<<G4cout;
// }
// }
void SetElementFieldValueIfNeeded(G4int eventNr);
virtual G4double GetNominalFieldValue() = 0;
virtual void SetNominalFieldValue(G4double newFieldValue) =0;
protected:
G4LogicalVolume* lvolume;
G4AffineTransform global2local;
// F04ElementField(const F04ElementField&);
private:
static G4Navigator* aNavigator;
G4String color;
G4ThreeVector center;
G4double minX, minY, minZ, maxX, maxY,maxZ;
G4double maxStep;
G4UserLimits* userLimits;
G4String elementFieldName;
std::map<G4int,G4double> changeFieldInStepsMap;
};
#endif

73
geant4/LEMuSR/include/F04FieldMessenger.hh Normal file
View File

@ -0,0 +1,73 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * 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. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
//
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
#ifndef F04FieldMessenger_h
#define F04FieldMessenger_h 1
#include "globals.hh"
#include "G4UImessenger.hh"
class F04GlobalField;
class G4UIdirectory;
class G4UIcmdWithAString;
class G4UIcmdWithAnInteger;
class G4UIcmdWithADoubleAndUnit;
class G4UIcmdWithoutParameter;
class F04FieldMessenger: public G4UImessenger
{
public:
F04FieldMessenger(F04GlobalField* );
~F04FieldMessenger();
void SetNewValue(G4UIcommand*, G4String);
void SetNewValue(G4UIcommand*, G4int);
private:
F04GlobalField* fGlobalField;
G4UIdirectory* detDir;
G4UIcmdWithAnInteger* fStepperCMD;
G4UIcmdWithADoubleAndUnit* fMinStepCMD;
G4UIcmdWithADoubleAndUnit* fDeltaChordCMD;
G4UIcmdWithADoubleAndUnit* fDeltaOneStepCMD;
G4UIcmdWithADoubleAndUnit* fDeltaIntersectionCMD;
G4UIcmdWithADoubleAndUnit* fEpsMinCMD;
G4UIcmdWithADoubleAndUnit* fEpsMaxCMD;
G4UIcmdWithoutParameter* fUpdateCMD;
};
#endif

193
geant4/LEMuSR/include/F04GlobalField.hh Normal file
View File

@ -0,0 +1,193 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * 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. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
//
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
#ifndef F04GlobalField_h
#define F04GlobalField_h 1
#include <vector>
#include "G4FieldManager.hh"
#include "G4PropagatorInField.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4ChordFinder.hh"
#include "G4MagneticField.hh"
#include "G4ElectroMagneticField.hh"
#include "G4Mag_EqRhs.hh"
#include "G4Mag_SpinEqRhs.hh"
#include "G4EqMagElectricField.hh"
#include "G4EqEMFieldWithSpin.hh"
#include "F04FieldMessenger.hh"
#include "F04ElementField.hh"
// F04GlobalField - handles the global ElectroMagnetic field
//
// There is a single G04GlobalField object.
//
// The field from each individual beamline element is given by a
// ElementField object. Any number of overlapping ElementField
// objects can be added to the global field. Any element that
// represents an element with an EM field must add the appropriate
// ElementField to the global GlobalField object.
typedef std::vector<F04ElementField*> FieldList;
class F04GlobalField : public G4ElectroMagneticField {
//class F04GlobalField : public G4MagneticField {
private:
F04GlobalField();
F04GlobalField(const F04GlobalField&);
~F04GlobalField();
F04GlobalField& operator=(const F04GlobalField&);
void setupArray();
public:
/// getObject() returns the single F04GlobalField object.
/// It is constructed, if necessary.
static F04GlobalField* getObject();
//
static G4bool Exists() {if (object==NULL) {return false;} else {return true;}; }
/// GetFieldValue() returns the field value at a given point[].
/// field is really field[6]: Bx,By,Bz,Ex,Ey,Ez.
/// point[] is in global coordinates: x,y,z,t.
void GetFieldValue(const G4double* point, G4double* field) const;
/// DoesFieldChangeEnergy() returns true.
G4bool DoesFieldChangeEnergy() const { return true; }
/// addElementField() adds the ElementField object for a single
/// element to the global field.
void addElementField(F04ElementField* f) { if (fields) fields->push_back(f); }
/// clear() removes all ElementField-s from the global object,
/// and destroys them. Used before the geometry is completely
/// re-created.
void clear();
/// updates all field tracking objects and clear()
void updateField();
/// Set the Stepper types
void SetStepperType( G4int i ) { fStepperType = i; }
/// Set the Stepper
void SetStepper();
/// Set the minimum step length
void SetMinStep(G4double s) { minStep = s; }
/// Set the delta chord length
void SetDeltaChord(G4double s) { deltaChord = s; }
/// Set the delta one step length
void SetDeltaOneStep(G4double s) { deltaOneStep = s; }
/// Set the delta intersection length
void SetDeltaIntersection(G4double s) { deltaIntersection = s; }
/// Set the minimum eps length
void SetEpsMin(G4double s) { epsMin = s; }
/// Set the maximum eps length
void SetEpsMax(G4double s) { epsMax = s; }
/// Return the list of Element Fields
FieldList* getFields() { return fields; }
// Additions by KS
// G4bool DoesAnyFieldValueNeedsToBeChanged(G4int eventNumber) {return globalChangeFieldInStepsMap[eventNumber];}
void CheckWhetherAnyNominalFieldValueNeedsToBeChanged(G4int eventNumber);
// Add point, at which user wishes to print out the field value
void AddPointForFieldTesting(G4ThreeVector point) {pointsAtWhichUserWantsToPrintFieldValue.push_back(point);}
// Print field value at all points user wished to be print out:
void PrintFieldAtRequestedPoints() const;
protected:
/// Get the global field manager
G4FieldManager* GetGlobalFieldManager();
private:
static F04GlobalField* object;
G4int nfp;
G4bool first;
FieldList* fields;
const F04ElementField **fp;
std::vector<G4ThreeVector> pointsAtWhichUserWantsToPrintFieldValue;
private:
G4int fStepperType;
G4double minStep;
G4double deltaChord;
G4double deltaOneStep;
G4double deltaIntersection;
G4double epsMin;
G4double epsMax;
// G4Mag_EqRhs* fEquation;
// G4Mag_SpinEqRhs* fEquation;
// G4EqMagElectricField* fEquation;
G4EqEMFieldWithSpin* fEquation;
G4FieldManager* fFieldManager;
G4PropagatorInField* fFieldPropagator;
G4MagIntegratorStepper* fStepper;
G4ChordFinder* fChordFinder;
F04FieldMessenger* fFieldMessenger;
std::map<G4int,G4bool> globalChangeFieldInStepsMap;
};
#endif

66
geant4/LEMuSR/include/MuDecayChannel.hh Normal file
View File

@ -0,0 +1,66 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * 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. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: G4MuonDecayChannel.hh,v 1.6 2006/06/29 19:23:35 gunter Exp $
// GEANT4 tag $Name: geant4-09-00 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 30 May 1997 H.Kurashige
// ------------------------------------------------------------
#ifndef MuDecayChannel_h
#define MuDecayChannel_h 1
#include "G4ios.hh"
#include "globals.hh"
#include "G4VDecayChannel.hh"
class MuDecayChannel :public G4VDecayChannel
{
// Class Decription
// This class describes muon decay kinemtics.
// This version neglects muon polarization
// assumes the pure V-A coupling
// gives incorrect energy spectrum for neutrinos
//
public: // With Description
//Constructors
MuDecayChannel(const G4String& theParentName,
G4double theBR);
// Destructor
virtual ~MuDecayChannel();
public: // With Description
virtual G4DecayProducts *DecayIt(G4double);
};
#endif

125
geant4/LEMuSR/include/MuDecayChannelWithSpin.hh Normal file
View File

@ -0,0 +1,125 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * 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. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// ------------------------------------------------------------
// 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 MuDecayChannelWithSpin_hh
#define MuDecayChannelWithSpin_hh 1
#include "globals.hh"
#include "G4ThreeVector.hh"
#include "MuDecayChannel.hh"
class MuDecayChannelWithSpin : public MuDecayChannel
{
// 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
MuDecayChannelWithSpin(const G4String& theParentName,
G4double theBR);
// Destructor
virtual ~MuDecayChannelWithSpin();
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 MuDecayChannelWithSpin::SetPolarization(G4ThreeVector polar)
{
parent_polarization = polar;
}
inline const G4ThreeVector& MuDecayChannelWithSpin::GetPolarization() const
{
return parent_polarization;
}
inline G4double MuDecayChannelWithSpin::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 MuDecayChannelWithSpin::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

72
geant4/LEMuSR/include/lem4Axial2DElField.hh Normal file
View File

@ -0,0 +1,72 @@
#ifndef LEM3AXIAL2DELFIELD_H
#define LEM3AXIAL2DELFIELD_H 1
#include "F04ElementField.hh"
#include "F04GlobalField.hh" /// ATTENTION: USE GENERAL EM GLOBAL AND NOT MAGN. FIELD ONLY
///#include "G4MagneticField.hh"
///#include "G4ElectricField.hh"
#include "globals.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
class lem4Axial2DElField : public F04ElementField ///G4ElectricField
{
//! Contructor from a single field map
// lem4ElectricField(G4double fieldval,G4String file,G4String map_length_unit,
//G4double Offset, G4double nx, G4double ny, G4double nz);
///lem4Axial2DElField(const G4String filename, double fieldValue, double lenUnit, double fieldNormalisation, double offset);
public:
// Class constructor for 2D axial field map (x, z, Er, Ez)
lem4Axial2DElField(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
// "lenUnit" is the unit in which the grid coordinates are specified in the table
// "fieldNormalisation" is the normalisation that has to be applied on the field values in the table
// such that when applying V_L3 = 1 kV the E values coincide with those in the table
// "fieldValue" is the field value (in kV) set by the user (i.e. values normalised to 1 kV will be
// multiplied by this value).
// Virtual destructor
virtual ~lem4Axial2DElField() {};
// addFieldValue() adds the field for THIS particular map into field[].
// point[] is expressed in GLOBAL coordinates.
void addFieldValue(const G4double point[4], G4double* field) const;
// Usual Set and Get functions
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
// getWidth(), getHeight(), getLength(), return the dimensions of the field
// (used to define the boundary of the field)
virtual G4double getWidth() { return 2*dr; } // x coordinate
virtual G4double getHeight() { return 2*dr; } // y coordinate
virtual G4double getLength() { return 2*dz; } // z coordinate
///void GetFieldValue(const double Point[4], double *Efield) const;
///void SetFieldValue(double newFieldValue);
/// G4double GetFieldSetValue();
private:
// Storage space for the table
std::vector < std::vector < double > > rField;
std::vector < std::vector < double > > zField;
// The dimensions of the table
int nr, nz;
// The physical limits of the defined region and the global offset
double minr, maxr, minz, maxz;
// The physical extent of the defined region
double dr, dz;
double ffieldValue; ///, zOffset; included in log_vol offset!
///bool invertR, invertZ; // substituted by the Invert function
// Utility function for inverting field map
void Invert(const char* indexToInvert);
};
#endif

81
geant4/LEMuSR/include/lem4DetectorConstruction.hh Normal file
View File

@ -0,0 +1,81 @@
#ifndef lem4DetectorConstruction_h
#define lem4DetectorConstruction_h 1
#include "globals.hh"
#include "G4VUserDetectorConstruction.hh"
#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4FieldManager.hh"
#include <map>
class G4Tubs;
class G4Box;
class G4Cons;
class G4Trd;
class G4LogicalVolume;
class G4VPhysicalVolume;
class G4Material;
class lem4DetectorMessenger;
//class lem4TrackerSD;
//class lem4ShieldSD;
//class lem4IForwardSD;
//class lem4IBackwardSD;
class lem4ScintSD;
//class lem4UniformField; // Added by TS to implement uniform box field according to Gumplinger
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4DetectorConstruction : public G4VUserDetectorConstruction
{
public:
lem4DetectorConstruction();
~lem4DetectorConstruction();
public:
G4VPhysicalVolume* Construct();
void UpdateGeometry();
///ALL THESE ARE OLD FIELD SETTING METHODS. USE THAT OF GUMPLINGER INSTEAD! TS.
void SetMagField(const char* fieldTableType, const char* inputFileName, G4double fieldValue);
void SetUniformMagField(G4double);
void SetGaussianMagField(G4ThreeVector);
///void SetUniformElField(G4double); // Uniform Electric Field. TS
void SetInputParameterFileName(G4String fileName) {parameterFileName=fileName;};
void ReportGeometryProblem(char myString[501]);
void ReportProblemInStearingFile(char* myString); //!Changed -added
G4Material* CharToMaterial(char myString[100]);
G4LogicalVolume* FindLogicalVolume(G4String LogicalVolumeName);
void SetColourOfLogicalVolume(G4LogicalVolume* pLogVol,char* colour);
G4bool CheckIfStringContainsSubstring(char* string, char* substring); //!Changed -added
private:
G4String parameterFileName; // name of the file with the geometry parameters
G4bool checkOverlap; // parameter to check ovelaping volumes
G4double largestAcceptableStep; // parameter defining largest step in the magnetic field
//G4Material* TargetMater;// pointer to the target material
//G4Material* VacMater; // pointer to the target material
//G4Material* DetMater; // pointer to the target material
//G4Material* Al; // pointer to the target material
//G4Material* Pb; // pointer to the target material
//G4Material* AirMater; // pointer to the target material
//G4Material* Vacuum; // pointer to vacuum material
lem4ScintSD* aScintSD;
lem4DetectorMessenger* detectorMessenger; // pointer to the Messenger
std::map<std::string,G4RotationMatrix*> pointerToRotationMatrix;
std::map<std::string,G4FieldManager*> pointerToField;
private:
void DefineMaterials();
};
#endif

57
geant4/LEMuSR/include/lem4DetectorMessenger.hh Normal file
View File

@ -0,0 +1,57 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4DetectorMessenger_h
#define lem4DetectorMessenger_h 1
#include "globals.hh"
#include "G4UImessenger.hh"
class lem4DetectorConstruction;
class G4UIdirectory;
class G4UIcmdWithAString;
class G4UIcmdWithADoubleAndUnit;
class G4UIcmdWithAnInteger;
class G4UIcmdWithoutParameter;
class G4UIcmdWith3Vector;
class G4UIcmdWith3VectorAndUnit;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4DetectorMessenger: public G4UImessenger
{
public:
lem4DetectorMessenger(lem4DetectorConstruction*);
~lem4DetectorMessenger();
void SetNewValue(G4UIcommand*, G4String);
private:
lem4DetectorConstruction* myDetector;
G4UIdirectory* lem4Dir;
G4UIdirectory* detDir;
G4UIdirectory* runDir;
G4UIcmdWithAString* CommandCmd;
G4UIcmdWithAnInteger* RunIDSetCmd;
G4UIcmdWithAnInteger* RandomOptionCmd;
G4UIcmdWithAnInteger* HowOftenToPrintEventCmd;
// G4UIcmdWithAString* WhichProcessesCmd;
G4UIcmdWithoutParameter* UpdateCmd;
// G4UIcmdWithADoubleAndUnit* FieldCmd;
G4UIcmdWithADoubleAndUnit* UFieldCmd;
/// G4UIcmdWithADoubleAndUnit* UEFieldCmd; //Uniform Electric Field. TS
G4UIcmdWith3Vector* GFieldCmd;
public:
// cks: The two variables used for the random number initialisation each event (if required)
// long myEventNr;
// static G4bool setRandomNrSeedAccordingEventNr;
// void IncrementMyEventNr() {myEventNr++};
// long GetMyEventNr() {return myEventNr};
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

31
geant4/LEMuSR/include/lem4ErrorMessage.hh Normal file
View File

@ -0,0 +1,31 @@
#ifndef lem4ErrorMessage_h
#define lem4ErrorMessage_h 1
#include <map>
#include "globals.hh"
enum SEVERITY {INFO, WARNING, SERIOUS, FATAL};
typedef struct
{
SEVERITY mesSeverity;
int nTimes;
} ErrorStruct;
class lem4ErrorMessage {
public:
// enum SEVERITY {INFO, WARNING, SERIOUS, FATAL};
lem4ErrorMessage();
~lem4ErrorMessage();
static lem4ErrorMessage* GetInstance();
void lem4Error(SEVERITY severity, G4String message, G4bool silent);
void PrintErrorSummary();
private:
static lem4ErrorMessage* pointerToErrors;
G4int nErrors;
// std::map<G4String,int> ErrorMapping;
std::map<G4String,ErrorStruct> ErrorMapping;
std::map<SEVERITY,G4String> severityWord;
};
#endif

67
geant4/LEMuSR/include/lem4EventAction.hh Normal file
View File

@ -0,0 +1,67 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4EventAction_h
#define lem4EventAction_h 1
#include "globals.hh"
#include "G4UserEventAction.hh"
#include <vector>
using namespace std;
class G4Timer;
class G4Event;
class lem4MagneticField;
class G4ElectricField; // Unif. el. field. TS
class lem4TabulatedField3D;
class lem4TabulatedField2D;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4EventAction : public G4UserEventAction
{
public:
lem4EventAction();
~lem4EventAction();
public:
void BeginOfEventAction(const G4Event*);
void EndOfEventAction(const G4Event*);
static lem4EventAction* GetInstance();
private:
// pointer to this class
static lem4EventAction* pointer;
//
G4Timer* timer;
static vector<int> * pointerToSeedVector;
// Variables for the time-dependent magnetic field
G4bool timeDependentField;
G4double fieldValueStart, fieldValueEnd, fieldStep, lastFieldValue;
G4int fieldNrOfSteps;
G4int maxEventNr;
lem4MagneticField* pointerToMusrUniformField;
G4ElectricField* pointerToMusrUniformEField; // Pointer to uniform electric field. TS
lem4TabulatedField3D* pointerToTabulatedField3D;
lem4TabulatedField2D* pointerToTabulatedField2D;
G4int latestEventNr;
public:
static G4bool setRandomNrSeedAccordingEventNr;
static G4bool setRandomNrSeedFromFile;
static G4int nHowOftenToPrintEvent;
static vector<int> * GetPointerToSeedVector();
// void setMyEventNr(long n) {myEventNr=n;};
void SetTimeDependentField(G4bool setFieldToBeTimeDependend, G4double initialField,
G4double finalField, G4int nrOfSteps);
void SetNumberOfEventsInThisRun(G4int nrOfEventsInThisRun) {maxEventNr=nrOfEventsInThisRun;};
void SetPointerToMusrUniformField(lem4MagneticField* pointer) {pointerToMusrUniformField=pointer;};
void SetPointerToMusrUniformEField(G4ElectricField* pointer) {pointerToMusrUniformEField=pointer;}; // Unif. El. field. TS
void SetPointerToTabulatedField3D(lem4TabulatedField3D* pointer) {pointerToTabulatedField3D=pointer;};
void SetPointerToTabulatedField2D(lem4TabulatedField2D* pointer) {pointerToTabulatedField2D=pointer;};
G4int GetLatestEventNr(){return latestEventNr;};
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

16
geant4/LEMuSR/include/lem4GaussianField.hh Normal file
View File

@ -0,0 +1,16 @@
#include "globals.hh"
#include "G4MagneticField.hh"
class lem4GaussianField
#ifndef STANDALONE
: public G4MagneticField
#endif
{
private:
G4double ffieldValue;
G4double fieldMagnitude,RMS;
public:
lem4GaussianField(double fieldValue, double sigma);
void GetFieldValue( const double point[4], double *Bfield ) const;
G4double GetFieldSetValue();
};

38
geant4/LEMuSR/include/lem4MagneticField.hh Normal file
View File

@ -0,0 +1,38 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4MagneticField_H
#define lem4MagneticField_H
#include "G4UniformMagField.hh"
class G4FieldManager;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4MagneticField: public G4UniformMagField
{
public:
lem4MagneticField(G4ThreeVector); // The value of the field
lem4MagneticField(); // A zero field
~lem4MagneticField();
//Set the field (fieldValue,0,0)
void SetFieldValue(G4double fieldValue);
void SetFieldValue(G4ThreeVector fieldVector);
G4ThreeVector GetConstantFieldValue();
protected:
// Find the global Field Manager
G4FieldManager* GetGlobalFieldManager(); // static
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif
G4bool CheckOverlaps(G4int res=1000, G4double tol=0., G4bool verbose=true);

75
geant4/LEMuSR/include/lem4MuFormation.hh Normal file
View File

@ -0,0 +1,75 @@
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
// Muonium Formation according to yield.cc function (through GetYields method).
// Id : lem4MuFormation.hh, v 1.4
// Author: Taofiq PARAISO, T. Shiroka
// Date : 2007-12
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
#ifndef lem4MuFormation_h
#define lem4MuFormation_h 1
#include "G4VDiscreteProcess.hh"
#include "G4ParticleTable.hh"
#include "yields.h"
/*! lem4MuFormation class defines the muonium formation process in the Carbon foil
* according to yields from Gonin's paper Sci. Rev. Instrum. 65(3), 648-652 (1994).
* \image html yields3.gif The muonium formation yields.
* The main parameters are the foil thickness and muon energy. For a given energy,
* a corresponding proportion of the muons will be converted into Muonium.
* Concretely, the muon is eliminated and replaced by a Muonium with identical
* properties, including time, energy, momentum, position etc.
*
* The process is executed at the END of a step, i.e. the muon is converted into
* Muonium AFTER flying through the Carbon foil (see also yields.h). */
class lem4MuFormation : public G4VDiscreteProcess
{
public:
lem4MuFormation(const G4String& name = "MuFormation", // process description
G4ProcessType aType = fElectromagnetic);
~lem4MuFormation();
//! - Main method. Muonium formation 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 GetDatas( const G4Step* aStep);
// model parameters
G4ParticleTable* particleTable;
G4ParticleDefinition* particle;
Yields Gonin;
G4double yvector[3];
G4double rnd;
G4DynamicParticle *DP;
//! The particle change object.
G4VParticleChange fParticleChange;
void PrepareSecondary(const G4Track&);
G4Track* aSecondary;
void InitializeSecondaries(const G4Track&);
};
#endif

52
geant4/LEMuSR/include/lem4MuScatter.hh Normal file
View File

@ -0,0 +1,52 @@
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
// Muonium "Scattering"
// Id : lem4MuScatter.hh, v 1.4
// Author: Taofiq PARAISO, T. Shiroka
// Date : 2007-12
// Notes : Simplified model for Mu scattering. Spin effects have been excluded.
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
/*!
* lem4MuScatter class defines the Muonium scattering process. It implements a very
* basic model which assumes Muonium looses its electron as soon as it enters any
* material (except for vacuum and CFoil). The class has only a PostStepDoIt method.
* The in-flight Muonium spin precession has been supressed. */
#ifndef lem4MuScatter_h
#define lem4MuScatter_h 1
#include "G4VDiscreteProcess.hh"
class lem4MuScatter : public G4VDiscreteProcess
{
public:
lem4MuScatter(const G4String& name="MuScatt", // process description
G4ProcessType aType = fGeneral);
~lem4MuScatter();
//! \mm The actions are taken at the end of the step.
G4VParticleChange* PostStepDoIt(const G4Track&,
const G4Step&);
G4double GetMeanFreePath(const G4Track& aTrack,
G4double previousStepSize,
G4ForceCondition* condition);
//! The condition for applying the process.
G4bool CheckCondition(const G4Step& aStep);
G4bool condition;
G4double itime, gtime, ftime,deltatime;
G4String p_name;
G4DynamicParticle *DP;
G4ParticleChange fParticleChange;
void PrepareSecondary(const G4Track&);
G4Track* aSecondary;
void InitializeSecondaries(const G4Track&);
};
#endif

72
geant4/LEMuSR/include/lem4Muonium.hh Normal file
View File

@ -0,0 +1,72 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * 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. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: lem4Muonium.hh,v 1.10 2006/06/29 19:20:21 gunter Exp $
// GEANT4 tag $Name: geant4-09-00 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 4-th April 1996, G.Cosmo
// ****************************************************************
// New implementation as a utility class M.Asai, 26 July 2004
// ----------------------------------------------------------------
#ifndef lem4Muonium_h
#define lem4Muonium_h 1
#include "globals.hh"
#include "G4ios.hh"
#include "G4ParticleDefinition.hh"
// ######################################################################
// ### MUONIUM ###
// ######################################################################
class lem4Muonium : public G4ParticleDefinition
{
private:
static lem4Muonium* theInstance;
lem4Muonium(){}
~lem4Muonium(){}
public:
static lem4Muonium* Definition();
static lem4Muonium* MuoniumDefinition();
static lem4Muonium* Muonium();
};
#endif

47
geant4/LEMuSR/include/lem4Parameters.hh Normal file
View File

@ -0,0 +1,47 @@
#ifndef lem4Parameters_h
#define lem4Parameters_h 1
#include "globals.hh"
class lem4Parameters {
public:
///lem4Parameters();
lem4Parameters(G4String steeringFileName); // ADDED by TS - overloaded constructor
~lem4Parameters();
static lem4Parameters* GetInstance();
void SetMyTypeOfProcesses(G4String string) {typeOfProcesses=string;};
G4String GetMyTypeOfProcesses() {return typeOfProcesses;};
static G4bool storeOnlyEventsWithHits; // Variable specifying whether to store interesting
// or all events into the ROOT tree. (default = true)
static G4double signalSeparationTime; // Minimim time separation between two subsequent signal
static G4bool storeOnlyTheFirstTimeHit; // If true, only the hit that happened first will be
// stored, anything else will be ignored
// (useful in some special cases, not for a routine simulation)
static G4bool includeMuoniumProcesses; // If true, includes Muonium formation and all
// other Mu-related processes in the simulation
//ADDED by TS
static G4bool boolG4GeneralParticleSource;//If true, G4GeneralParticleSource will be initialised
//instead of G4ParticleGun - for simulating radioactive sources
static G4bool field_DecayWithSpin; // If true, then the routine for calculating the magnetic
// field will use a more precise argument.
// This variable is set to "true" by the SteppinAction
// and reset to "false" in the GetFieldValue.
// It is being changed on step by step basis.
static G4int nrOfEventsToBeGenerated; // Nr of events to be simulated in this run (set by /run/beamOn command)
private:
static lem4Parameters* pointerToParameters;
G4String typeOfProcesses; // Variable defining what kind of Physics processes to call:
// - Standard EM
// - LowEnergy (default)
// - Penelope
// - Coulomb
};
#endif

48
geant4/LEMuSR/include/lem4PhysicsList.hh Normal file
View File

@ -0,0 +1,48 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4PhysicsList_h
#define lem4PhysicsList_h 1
#include "G4VUserPhysicsList.hh"
#include "globals.hh"
//cks Added to have Geant default muon decay with spin
#include "G4DecayWithSpin.hh"
//#include "lem4DecayWithSpin.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4PhysicsList: public G4VUserPhysicsList
{
public:
lem4PhysicsList();
~lem4PhysicsList();
protected:
// Construct particle and physics
void ConstructParticle();
void ConstructProcess();
void SetCuts();
protected:
// these methods Construct particles
void ConstructBosons();
void ConstructLeptons();
void ConstructMesons();
void ConstructBaryons();
protected:
// these methods Construct physics processes and register them
void ConstructGeneral();
void ConstructEM();
// private:
// char myProcesses[100];
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

69
geant4/LEMuSR/include/lem4PrimaryGeneratorAction.hh Normal file
View File

@ -0,0 +1,69 @@
#ifndef lem4PrimaryGeneratorAction_h
#define lem4PrimaryGeneratorAction_h 1
#include "G4VUserPrimaryGeneratorAction.hh"
#include "globals.hh"
#include "Randomize.hh"
#include "G4ThreeVector.hh"
#include <stdio.h>
//cksdel #include "lem4ParticleGun.hh"
class G4ParticleGun;
class G4Event;
class lem4DetectorConstruction;
class lem4PrimaryGeneratorMessenger;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
class lem4PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
{
public:
lem4PrimaryGeneratorAction(lem4DetectorConstruction*);
~lem4PrimaryGeneratorAction();
public:
void GeneratePrimaries(G4Event*);
void SetRndmFlag(G4String val) {rndmFlag = val;}
void Setvertex(G4ThreeVector v) {x0=v[0]; y0=v[1]; z0=v[2];}
void SetvertexSigma(G4ThreeVector v) {xSigma=v[0]; ySigma=v[1]; zSigma=v[2];}
void SetvertexBoundary(G4ThreeVector v) {rMaxAllowed=v[0]; zMinAllowed=v[1]; zMaxAllowed=v[2];}
void SetEnergy(G4double val) {p0=std::sqrt(val*val + 2*mu_mass*val);} // {E0=val;}
void SetMomentum(G4double val) {p0=val;}
void SetMomentumSmearing(G4double val) {pSigma=val;}
void SetMomentumBoundary(G4ThreeVector v){pMinAllowed=v[0]; pMaxAllowed=v[1];}
void SetTilt(G4ThreeVector v) {xangle0=v[0]; yangle0=v[1];}
void SetSigmaTilt(G4ThreeVector v) {xangleSigma=v[0]; yangleSigma=v[1];}
void SetPitch(G4double val) {pitch=val;}
void SetInitialMuonPolariz(G4ThreeVector vIniPol);
void SetMuonDecayTimeLimits(G4ThreeVector decayTimeLimits);
void SetTurtleInput(G4String turtleFileName);
static G4String GetPrimaryName();
private:
G4ParticleGun* particleGun; // Pointer a to G4 service class
lem4DetectorConstruction* lem4Detector; // Pointer to the geometry
lem4PrimaryGeneratorMessenger* gunMessenger; // Messenger of this class
G4String rndmFlag; // Flag for a random impact point
G4double x0, y0, z0, xSigma, ySigma, zSigma, rMaxAllowed, zMinAllowed, zMaxAllowed;
G4double E0, p0, pSigma, pMinAllowed, pMaxAllowed;
G4double xangle0, yangle0, xangleSigma, yangleSigma, pitch;
G4bool UnpolarisedMuonBeam;
G4double xPolarisIni, yPolarisIni, zPolarisIni;
G4double muonDecayTimeMin;
G4double muonDecayTimeMax;
G4double muonMeanLife;
G4double mu_mass;
// For Turtle input
FILE* fTurtleFile;
G4bool takeMuonsFromTurtleFile;
static G4String thePrimaryParticleName;
public:
G4double decaytime;
};
#endif

41
geant4/LEMuSR/include/lem4PrimaryGeneratorMessenger.hh Normal file
View File

@ -0,0 +1,41 @@
#ifndef lem4PrimaryGeneratorMessenger_h
#define lem4PrimaryGeneratorMessenger_h 1
#include "G4UImessenger.hh"
#include "globals.hh"
class lem4PrimaryGeneratorAction;
class G4UIcmdWithAString;
class G4UIcmdWithADoubleAndUnit;
class G4UIcmdWithADouble;
class G4UIcmdWithAnInteger;
class G4UIcmdWith3VectorAndUnit;
class G4UIcmdWith3Vector;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
class lem4PrimaryGeneratorMessenger: public G4UImessenger
{
public:
lem4PrimaryGeneratorMessenger(lem4PrimaryGeneratorAction* lem4Gun);
~lem4PrimaryGeneratorMessenger();
void SetNewValue(G4UIcommand*, G4String newValue);
private:
lem4PrimaryGeneratorAction* lem4Action;
G4UIcmdWith3VectorAndUnit* setvertexCmd;
G4UIcmdWith3VectorAndUnit* setvertexSigmaCmd;
G4UIcmdWith3VectorAndUnit* setvertexBoundaryCmd;
G4UIcmdWithADoubleAndUnit* setEnergyCmd;
G4UIcmdWithADoubleAndUnit* setMomentumCmd;
G4UIcmdWithADoubleAndUnit* setMomentumSmearingCmd;
G4UIcmdWith3VectorAndUnit* setMomentumBoundaryCmd;
G4UIcmdWith3VectorAndUnit* setTiltAngleCmd;
G4UIcmdWith3VectorAndUnit* setSigmaTiltAngleCmd;
G4UIcmdWithADoubleAndUnit* setPitchCmd;
G4UIcmdWith3Vector* setMuonPolarizCmd;
G4UIcmdWith3VectorAndUnit* setMuonDecayTimeCmd;
G4UIcmdWithAString* setTurtleCmd;
};
#endif

283
geant4/LEMuSR/include/lem4RootOutput.hh Normal file
View File

@ -0,0 +1,283 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4RootOutput_h
#define lem4RootOutput_h 1
//#include "G4UserRunAction.hh"
#include "globals.hh"
#include "G4ThreeVector.hh"
// ROOT
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TVectorD.h"
//
#include <map>
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4RootOutput {
public:
lem4RootOutput();
~lem4RootOutput();
static lem4RootOutput* GetRootInstance();
public:
void BeginOfRunAction();
void EndOfRunAction();
void FillEvent();
void ClearAllRootVariables();
// void SetSensitiveDetectorMapping(std::string logivol, int volumeID);
void SetVolumeIDMapping(std::string logivol, int volumeID);
G4int ConvertVolumeToID(std::string logivol);
void SetSpecialSaveVolumeDefined() {boolIsAnySpecialSaveVolumeDefined=true;};
// Getting variables (just for debugging)
G4double GetDecayPositionZ() {return muDecayPosZ_t;};
G4double GetDecayTime() {return muDecayTime_t*microsecond;};
G4double GetTimeInTarget() {return muTargetTime_t*microsecond;};
// Setting variables common to the whole event:
void SetRunID (G4int id) {runID_t = id;};
void SetEventID (G4int id) {eventID_t = id;};
void SetDecayDetectorID (std::string detectorName) {muDecayDetID_t = SensDetectorMapping[detectorName];};
void SetBField (G4double F[6]) {B_t[0]=F[0]/tesla; B_t[1]=F[1]/tesla; B_t[2]=F[2]/tesla;
B_t[3]=F[3]/tesla; B_t[4]=F[4]/tesla; B_t[5]=F[5]/tesla;};
void SetDecayPolarisation (G4ThreeVector pol) {muDecayPolX_t=pol.x(); muDecayPolY_t=pol.y(); muDecayPolZ_t=pol.z();};
void SetDecayPosition (G4ThreeVector pos) {muDecayPosX_t=pos.x()/mm; muDecayPosY_t=pos.y()/mm;
muDecayPosZ_t=pos.z()/mm;};
void SetGlobTime (G4double gt) {globalTime_t = gt/microsecond;};
void SetDetectorInfo (G4int nDetectors, G4int ID, G4double edep, G4double edep_el, G4double edep_pos,
G4double edep_gam, G4double edep_mup,G4int nsteps, G4double length, G4double t1,
G4double t2, G4double x, G4double y, G4double z) {
if ((nDetectors<0)||(nDetectors>=(det_nMax-1))) {
G4cout<<"SERIOUS ERROR! lem4RootOutput.hh: nDetectors="<<nDetectors<<
" is larger than det_nMax="<<det_nMax<<G4endl;
return;
}
else {
det_n=nDetectors+1; det_ID[nDetectors]=ID; det_edep[nDetectors]=edep/MeV;
det_edep_el[nDetectors]=edep_el/MeV; det_edep_pos[nDetectors]=edep_pos/MeV;
det_edep_gam[nDetectors]=edep_gam/MeV; det_edep_mup[nDetectors]=edep_mup/MeV;
det_nsteps[nDetectors]=nsteps; det_length[nDetectors]=length/mm;
det_time_start[nDetectors]=t1/microsecond; det_time_end[nDetectors]=t2/microsecond;
det_x[nDetectors]=x/mm;det_y[nDetectors]=y/mm;det_z[nDetectors]=z/mm;
// G4cout << "Saved to root: eventID="<<eventID_t<<"; det_ID="<<det_ID[nDetectors]
// <<"; det_edep="<<det_edep[nDetectors]<<G4endl;
}
//G4cout << "eventID="<<eventID_t<<"; det_n="<<det_n<<"; det_ID="<<det_ID[nDetectors]<<"; "
// << "Deposited E (el, pos, gam, mu)="<<det_edep_el[nDetectors]<<", "
// << det_edep_pos[nDetectors]<<", "<< det_edep_gam[nDetectors]<<", "
// << det_edep_mup[nDetectors]<<G4endl;
//G4cout << " det_time_start="<<det_time_start[nDetectors]
// << " det_time_end="<<det_time_end[nDetectors] << G4endl;
};
void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z,
G4double px, G4double py, G4double pz) ;
void SetInitialMuonParameters(G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz,
G4double xpolaris, G4double ypolaris, G4double zpolaris) {
muIniPosX_t=x; muIniPosY_t=y; muIniPosZ_t=z;
muIniMomX_t=px; muIniMomY_t=py; muIniMomZ_t=pz;
muIniPolX_t=xpolaris; muIniPolY_t=ypolaris; muIniPolZ_t=zpolaris;
}
// void SetSubtrackInfo (G4int iSubtrack, G4int partID, G4int trID, G4int volID, G4double E,
// G4double xx, G4double yy, G4double zz, G4double xxPost, G4double yyPost, G4double zzPost,
// G4double Eki, G4double stepL, G4int nrOfSubst) {
// if ((iSubtrack<0)||(iSubtrack>=maxNSubTracks)) {
// G4cout<<"SERIOUS ERROR! lem4RootOutput.hh: iSubtrack="<<iSubtrack<<
// " is larger than maxNSubTracks="<<maxNSubTracks<<G4endl;
// return;
// }
// else {
// nSubTracks=iSubtrack+1;
// particleID_t[iSubtrack]=partID;
// trackID_t[iSubtrack]=trID;
// logicalVolumeID_t[iSubtrack]=volID;
// edep_t[iSubtrack]=E/MeV;
// x_t[iSubtrack]=xx/mm;
// y_t[iSubtrack]=yy/mm;
// z_t[iSubtrack]=zz/mm;
// post_x_t[iSubtrack]=xxPost/mm;
// post_y_t[iSubtrack]=yyPost/mm;
// post_z_t[iSubtrack]=zzPost/mm;
// kineticEnergy_t[iSubtrack]=Eki/MeV;
// stepLength_t[iSubtrack]=stepL/mm;
// nrOfSubsteps_t[iSubtrack]=nrOfSubst;
// }
// };
void SetPolInTarget(G4ThreeVector pol) {muTargetPolX_t=pol.x(); muTargetPolY_t=pol.y(); muTargetPolZ_t=pol.z();};
void SetTimeInTarget(G4double time) {muTargetTime_t = time/microsecond;};
void SetInitialPositronMomentum(G4ThreeVector mom) {posIniMomx_t=mom.x(); posIniMomy_t=mom.y(); posIniMomz_t=mom.z();};
void SetDecayTime(G4double time) {muDecayTime_t=time/microsecond;};
void SetFieldValue(G4double value) {fieldValue_t=value/tesla;};
void SetNrFieldNomVal(G4int n) {nFieldNomVal = n;}
void SetFieldNomVal(G4int i, G4double value);
G4int GetNrOfVolumes() {return det_nMax;}
void StoreGeantParameter(Int_t i, Double_t value) {
if (i<maxNGeantParameters) { GeantParametersD[i]=value; }
else {G4cout<<"lem4RootOutput.hh::StoreGeantParameter: index="<<i<<" out of range"
<<" (maxNGeantParameters=" <<maxNGeantParameters<<")"<<G4endl;}
};
TH2F *htest1, *htest2;
TH1F *htest3, *htest4, *htest5, *htest6, *htest7, *htest8;
public:
static G4bool store_runID;
static G4bool store_eventID;
static G4bool store_BFieldAtDecay;
static G4bool store_muIniPosX;
static G4bool store_muIniPosY;
static G4bool store_muIniPosZ;
static G4bool store_muIniMomX;
static G4bool store_muIniMomY;
static G4bool store_muIniMomZ;
static G4bool store_muIniPolX;
static G4bool store_muIniPolY;
static G4bool store_muIniPolZ;
static G4bool store_muDecayDetID;
static G4bool store_muDecayPosX;
static G4bool store_muDecayPosY;
static G4bool store_muDecayPosZ;
static G4bool store_muDecayTime;
static G4bool store_muDecayPolX;
static G4bool store_muDecayPolY;
static G4bool store_muDecayPolZ;
static G4bool store_muTargetTime;
static G4bool store_muTargetPolX;
static G4bool store_muTargetPolY;
static G4bool store_muTargetPolZ;
static G4bool store_posIniMomX;
static G4bool store_posIniMomY;
static G4bool store_posIniMomZ;
static G4bool store_globalTime;
static G4bool store_fieldValue;
static G4bool store_det_ID;
static G4bool store_det_edep;
static G4bool store_det_edep_el;
static G4bool store_det_edep_pos;
static G4bool store_det_edep_gam;
static G4bool store_det_edep_mup;
static G4bool store_det_nsteps;
static G4bool store_det_length;
static G4bool store_det_start;
static G4bool store_det_end;
static G4bool store_det_x;
static G4bool store_det_y;
static G4bool store_det_z;
static G4bool store_fieldNomVal;
private:
TFile* rootFile;
TTree* rootTree;
static lem4RootOutput* pointerToRoot;
static const Int_t maxNGeantParameters=30;
Double_t GeantParametersD[maxNGeantParameters]; // parameters transfered from GEANT to Root
// 0 ... fieldOption: 0 ... no field, 1 ... uniform, 2 ... gaussian, 3 ... from table
// 1 ... fieldValue: intensity of the magnetic field
// 2 ... minimum of the generated decay time of the muon (in microsecond)
// 3 ... maximum of the generated decay time of the muon (in microsecond)
// 4 ... muon mean life time (in microsecond)
// Variables common to the whole event:
Int_t runID_t;
Int_t eventID_t;
Double_t B_t[6];
Double_t muIniPosX_t, muIniPosY_t, muIniPosZ_t;
Double_t muIniMomX_t, muIniMomY_t, muIniMomZ_t;
Double_t muIniPolX_t, muIniPolY_t, muIniPolZ_t;
Int_t muDecayDetID_t;
Double_t muDecayPolX_t, muDecayPolY_t, muDecayPolZ_t;
Double_t muTargetPolX_t, muTargetPolY_t, muTargetPolZ_t;
Double_t muDecayPosX_t, muDecayPosY_t, muDecayPosZ_t;
Double_t muDecayTime_t;
Double_t posIniMomx_t, posIniMomy_t, posIniMomz_t;
Double_t globalTime_t, muTargetTime_t;
Double_t fieldValue_t;
public:
static const Int_t maxNFieldnNominalValues=30;
private:
Int_t nFieldNomVal;
Double_t fieldNomVal[maxNFieldnNominalValues];
// Variables for a particle in a given detector within the event
public:
static const Int_t maxNSubTracks=30;
private:
// Int_t nSubTracks;
// // Char_t particleName_t[6];
// Int_t particleID_t[maxNSubTracks];
// G4int trackID_t[maxNSubTracks];
// Double_t edep_t[maxNSubTracks];
// Double_t x_t[maxNSubTracks];
// Double_t y_t[maxNSubTracks];
// Double_t z_t[maxNSubTracks];
// Double_t post_x_t[maxNSubTracks];
// Double_t post_y_t[maxNSubTracks];
// Double_t post_z_t[maxNSubTracks];
// // Char_t logicalVolume_t[13];
// Int_t logicalVolumeID_t[maxNSubTracks];
// Double_t kineticEnergy_t[maxNSubTracks];
// Double_t stepLength_t[maxNSubTracks];
// Int_t nrOfSubsteps_t[maxNSubTracks];
// Variables for the activity inside a given detector
public:
static const Int_t det_nMax=100; // must be by 1 higher than the real number of detector "hits", because
// else the detector nr. 0 is counted (0 is used if no
// SensDetectorMapping correspond to a given logical volume).
private:
G4int det_n;
G4int det_ID[det_nMax];
G4double det_edep[det_nMax];
G4int det_nsteps[det_nMax];
G4double det_length[det_nMax];
G4double det_edep_el[det_nMax];
G4double det_edep_pos[det_nMax];
G4double det_edep_gam[det_nMax];
G4double det_edep_mup[det_nMax];
G4double det_time_start[det_nMax];
G4double det_time_end[det_nMax];
G4double det_x[det_nMax];
G4double det_y[det_nMax];
G4double det_z[det_nMax];
// G4bool firstStepInVolume_t;
// G4bool lastStepInVolume_t;
public:
static const Int_t save_nMax=1000;
private:
G4int save_n;
G4int save_detID[save_nMax];
G4int save_particleID[save_nMax];
G4double save_ke[save_nMax];
G4double save_x[save_nMax];
G4double save_y[save_nMax];
G4double save_z[save_nMax];
G4double save_px[save_nMax];
G4double save_py[save_nMax];
G4double save_pz[save_nMax];
G4bool boolIsAnySpecialSaveVolumeDefined;
std::map<std::string,int> SensDetectorMapping;
};
//G4int lem4RootOutput::maxNSubTracks=30;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

36
geant4/LEMuSR/include/lem4RunAction.hh Normal file
View File

@ -0,0 +1,36 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4RunAction_h
#define lem4RunAction_h 1
#include "G4UserRunAction.hh"
#include "globals.hh"
#include "lem4RootOutput.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class G4Timer;
class G4Run;
class lem4RunAction : public G4UserRunAction
{
public:
lem4RunAction();
~lem4RunAction();
public:
void BeginOfRunAction(const G4Run*);
void EndOfRunAction(const G4Run*);
private:
G4Timer* timer;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

128
geant4/LEMuSR/include/lem4ScintHit.hh Normal file
View File

@ -0,0 +1,128 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4ScintHit_h
#define lem4ScintHit_h 1
#include "G4VHit.hh"
#include "G4THitsCollection.hh"
#include "G4Allocator.hh"
#include "G4ThreeVector.hh"
#include "G4MagneticField.hh"
#include "globals.hh"
#include "G4ios.hh"
// ROOT
#include "TFile.h"
#include "TTree.h"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4ScintHit : public G4VHit
{
public:
lem4ScintHit();
~lem4ScintHit();
lem4ScintHit(const lem4ScintHit&);
const lem4ScintHit& operator=(const lem4ScintHit&);
G4int operator==(const lem4ScintHit&) const;
// bool operator() (lem4ScintHit hit1, lem4ScintHit hit2) { return (hit1.globalTime<hit2.globalTime);}
inline void* operator new(size_t);
inline void operator delete(void*);
void Draw();
void Print();
public:
void SetParticleName (G4String name) {particleName = name; };
void SetParticleID (G4int id) {particleID = id; };
void SetTrackID (G4int track) { trackID = track; };
void SetEdep (G4double de) { edep = de; };
void SetPrePos (G4ThreeVector xyz){ pre_pos = xyz; };
void SetPostPos (G4ThreeVector xyz){ post_pos = xyz; };
void SetPol (G4ThreeVector ijk){pol = ijk;};
void SetLogVolName (G4String logivol) {logicalVolume = logivol;};
void SetGlobTime (G4double gt) { globalTime = gt;};
void SetFirstStepInVolumeFlag (G4bool flag) { firstStepInVolume=flag;};
void SetLastStepInVolumeFlag (G4bool flag) { lastStepInVolume=flag;};
void SetKineticEnergy (G4double en) { kineticEnergy = en;};
void SetStepLength (G4double length) { stepLength = length;};
// void SetBField (G4double f[6]) {BF[0]=f[0]; BF[1]=f[1]; BF[2]=f[2]; BF[3]=f[3]; BF[4]=f[4]; BF[5]=f[5];};
void SetRunID(G4int i) {runID=i;};
void SetEventID(G4int i) {eventID=i;};
// void SetVerboseLevel (G4int n) { G4int lem4ScintHit::verboseLevel=n;};
G4String GetParticleName() {return particleName; };
G4int GetParticleID() {return particleID; };
G4int GetTrackID() { return trackID; };
G4double GetEdep() { return edep; };
G4ThreeVector GetPrePos(){ return pre_pos; };
G4ThreeVector GetPostPos(){ return post_pos; };
G4ThreeVector GetPol(){ return pol; };
G4String GetLogVolName() { return logicalVolume; };
G4double GetGlobalTime() { return globalTime; };
G4bool GetFirstStepInVolumeFlag() {return firstStepInVolume;};
G4bool GetLastStepInVolumeFlag() {return lastStepInVolume;};
G4double GetKineticEnergy() { return kineticEnergy; };
G4double GetStepLength() {return stepLength; };
G4double* GetBField () {return BF;};
G4int GetRunID() {return runID;};
G4int GetEventID() {return eventID;};
G4double point[4];
G4double B[6];
const G4Field *mfield;
private:
G4String particleName;
G4int particleID;
G4int trackID;
G4double edep;
G4ThreeVector pre_pos;
G4ThreeVector post_pos;
G4ThreeVector pol;
G4String logicalVolume;
G4double globalTime;
G4bool firstStepInVolume;
G4bool lastStepInVolume;
G4double kineticEnergy;
G4double stepLength;
G4int eventID;
G4int runID;
G4double BF[6];
static G4int ScintMultihit;
static G4int runIDoldScint;
static G4int eventIDoldScint;
static G4int NIS;
static G4int ScintChamberNbold;
static G4int verboseLevel;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
typedef G4THitsCollection<lem4ScintHit> lem4ScintHitsCollection;
extern G4Allocator<lem4ScintHit> lem4ScintHitAllocator;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
inline void* lem4ScintHit::operator new(size_t)
{
void *aHit;
aHit = (void *) lem4ScintHitAllocator.MallocSingle();
return aHit;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
inline void lem4ScintHit::operator delete(void *aHit)
{
lem4ScintHitAllocator.FreeSingle((lem4ScintHit*) aHit);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

42
geant4/LEMuSR/include/lem4ScintSD.hh Normal file
View File

@ -0,0 +1,42 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4ScintSD_h
#define lem4ScintSD_h 1
#include "G4VSensitiveDetector.hh"
#include "lem4ScintHit.hh"
#include "lem4RootOutput.hh"
//cksdel extern lem4RootOutput* myRootOutput;
class G4Step;
class G4HCofThisEvent;
//class lem4ScintHit;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4ScintSD : public G4VSensitiveDetector
{
public:
lem4ScintSD(G4String);
~lem4ScintSD();
void Initialize(G4HCofThisEvent*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
// bool timeOrdering (lem4ScintHit* hit1, lem4ScintHit* hit2);
// bool myREMOVEfunction (int i,int j) { return (i<j); }
private:
lem4ScintHitsCollection* scintCollection;
G4bool myStoreOnlyEventsWithHits;
G4double mySignalSeparationTime;
G4bool myStoreOnlyTheFirstTimeHit;
// G4bool Positron_momentum_already_stored;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

92
geant4/LEMuSR/include/lem4ShieldHit.hh Normal file
View File

@ -0,0 +1,92 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4ShieldHit_h
#define lem4ShieldHit_h 1
#include "G4VHit.hh"
#include "G4THitsCollection.hh"
#include "G4Allocator.hh"
#include "G4ThreeVector.hh"
#include "G4MagneticField.hh"
#include "globals.hh"
#include "G4ios.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4ShieldHit : public G4VHit
{
public:
lem4ShieldHit();
~lem4ShieldHit();
lem4ShieldHit(const lem4ShieldHit&);
const lem4ShieldHit& operator=(const lem4ShieldHit&);
G4int operator==(const lem4ShieldHit&) const;
inline void* operator new(size_t);
inline void operator delete(void*);
void Draw();
void Print();
public:
void SetTrackID (G4int track) { trackID = track; };
void SetParticleName (G4String name) {particle_name = name; };
void SetEdep (G4double de) { edep = de; };
void SetPos (G4ThreeVector xyz){ pos = xyz; };
void SetPol (G4ThreeVector ijk){pol = ijk;};
void SetLogVolName (G4String logivol) {logicalvolume = logivol;};
G4int GetTrackID() { return trackID; };
G4String GetParticleName() {return particle_name; };
G4double GetEdep() { return edep; };
G4ThreeVector GetPos(){ return pos; };
G4ThreeVector GetPol(){ return pol; };
G4String GetLogVolName() { return logicalvolume; };
G4double point[4];
G4double B[6];
const G4Field *mfield;
private:
G4int trackID;
G4String particle_name;
G4double edep;
G4ThreeVector pos;
G4ThreeVector pol;
G4String logicalvolume;
static G4int eventIDoldB;
static G4int Nb;
static G4int runIDoldB;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
typedef G4THitsCollection<lem4ShieldHit> lem4ShieldHitsCollection;
extern G4Allocator<lem4ShieldHit> lem4ShieldHitAllocator;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
inline void* lem4ShieldHit::operator new(size_t)
{
void *aHit;
aHit = (void *) lem4ShieldHitAllocator.MallocSingle();
return aHit;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
inline void lem4ShieldHit::operator delete(void *aHit)
{
lem4ShieldHitAllocator.FreeSingle((lem4ShieldHit*) aHit);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

32
geant4/LEMuSR/include/lem4ShieldSD.hh Normal file
View File

@ -0,0 +1,32 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4ShieldSD_h
#define lem4ShieldSD_h 1
#include "G4VSensitiveDetector.hh"
#include "lem4ShieldHit.hh"
class G4Step;
class G4HCofThisEvent;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4ShieldSD : public G4VSensitiveDetector
{
public:
lem4ShieldSD(G4String);
~lem4ShieldSD();
void Initialize(G4HCofThisEvent*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
private:
lem4ShieldHitsCollection* shieldCollection;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

47
geant4/LEMuSR/include/lem4SteppingAction.hh Normal file
View File

@ -0,0 +1,47 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#ifndef lem4SteppingAction_h
#define lem4SteppingAction_h 1
#include "G4UserSteppingAction.hh"
#include "G4ProcessManager.hh"
#include "globals.hh"
#include <fstream>
class G4Timer;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4SteppingAction : public G4UserSteppingAction
{
public:
static lem4SteppingAction* GetInstance();
lem4SteppingAction();
~lem4SteppingAction();
void UserSteppingAction(const G4Step *theStep);
void DoAtTheBeginningOfEvent();
void SetLogicalVolumeAsSpecialSaveVolume(G4String logicName, G4int volumeID);
private:
G4Timer* timer;
static lem4SteppingAction* pointer;
// G4int oldEventID;
G4ProcessManager* muPlusProcessManager;
G4bool multipleToCoulombScatteringIsPossible;
G4bool coulombScatteringIsActive;
size_t multipleScatteringIndex;
size_t coulombScatteringIndex;
G4bool muAlreadyWasInTargetInThisEvent;
G4bool radioactiveElectronAlreadySavedInThisEvent;
G4bool boolIsAnySpecialSaveVolumeDefined;
std::map<G4String,G4int> saveVolumeMapping;
G4String lastActualVolume;
// G4int debugOldEventID;
// G4int nrOfSteps;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

26
geant4/LEMuSR/include/lem4SteppingVerbose.hh Normal file
View File

@ -0,0 +1,26 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4SteppingVerbose;
#ifndef lem4SteppingVerbose_h
#define lem4SteppingVerbose_h 1
#include "G4SteppingVerbose.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4SteppingVerbose : public G4SteppingVerbose
{
public:
lem4SteppingVerbose();
~lem4SteppingVerbose();
void StepInfo();
void TrackingStarted();
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif

73
geant4/LEMuSR/include/lem4TabulatedElementField2D.hh Normal file
View File

@ -0,0 +1,73 @@
#ifndef lem4TabulatedElementField2D_h
#define lem4TabulatedElementField2D_h 1
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
#include "globals.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
class lem4TabulatedElementField2D : public F04ElementField
{
public: // with description
// Class constructor for 2D axial field map (r, z, EMr, EMz) - with EM = E or B; r - radial, z - longitudinal.
lem4TabulatedElementField2D(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
//
// "fieldType" is the type of EM field: electric - E, or magnetic - B
// "fieldValue" is the field to be applied (in T, or in kV/mm). The normalised field
// map values are multiplied by this value. The field-map itself has no units!
// "lenUnit" is the unit in which the grid coordinates of the field-map are specified
// "fieldNormalisation" is the normalisation factor that once applied to the tabulated field values
// satisfies the condition: (max. field value)*fieldNormalisation = 1
// To revert field direction, change its sign to negative.
// Virtual destructor
virtual ~lem4TabulatedElementField2D() {}
// addFieldValue() adds the field for THIS particular map into field[].
// point[] is expressed in GLOBAL coordinates.
void addFieldValue( const G4double Point[4], G4double* field) const;
// Usual Set and Get functions
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
// getWidth(), getHeight(), getLength(), return the dimensions of the field
// (used to define the boundary of the field)
virtual G4double getWidth() { return 2*dr; } // x coordinate
virtual G4double getHeight() { return 2*dr; } // y coordinate
virtual G4double getLength() { return dz; } // z coordinate
private:
// Storage space for the 2D table
std::vector< std::vector< double > > rField;
std::vector< std::vector< double > > zField;
// The field-map dimensions
int nr, nz;
// The field map Length unit (string and number)
G4String lUnit;
double lenUnit;
// The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla)
G4String fUnit;
double fieUnit;
// The field-map Field normalisation factor
double fieldNormalisation;
// The physical limits of the defined region
double minimumr, maximumr, minimumz, maximumz;
// The physical extent of the defined region
double dr, dz;
// See the description under the class constructor
char fldType;
double ffieldValue;
void Invert(const char* indexToInvert);
};
#endif

73
geant4/LEMuSR/include/lem4TabulatedElementField2Df.hh Normal file
View File

@ -0,0 +1,73 @@
#ifndef lem4TabulatedElementField2Df_h
#define lem4TabulatedElementField2Df_h 1
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
#include "globals.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
class lem4TabulatedElementField2Df : public F04ElementField
{
public: // with description
// Class constructor for r-axis FOLDED 2D axial field map (r, z, EMr, EMz) - with EM = E or B; r - radial, z - longit.
lem4TabulatedElementField2Df(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
//
// "fieldType" is the type of EM field: electric - E, or magnetic - B
// "fieldValue" is the field to be applied (in T, or in kV/mm). The normalised field
// map values are multiplied by this value. The field-map itself has no units!
// "lenUnit" is the unit in which the grid coordinates of the field-map are specified
// "fieldNormalisation" is the normalisation factor that once applied to the tabulated field values
// satisfies the condition: (max. field value)*fieldNormalisation = 1
// To revert field direction, change its sign to negative.
// Virtual destructor
virtual ~lem4TabulatedElementField2Df() {}
// addFieldValue() adds the field for THIS particular map into field[].
// point[] is expressed in GLOBAL coordinates.
void addFieldValue( const G4double Point[4], G4double* field) const;
// Usual Set and Get functions
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
// getWidth(), getHeight(), getLength(), return the dimensions of the field
// (used to define the boundary of the field)
virtual G4double getWidth() { return 2*dr; } // x coordinate
virtual G4double getHeight() { return 2*dr; } // y coordinate
virtual G4double getLength() { return 2*dz; } // z coordinate
private:
// Storage space for the 2D table
std::vector< std::vector< double > > rField;
std::vector< std::vector< double > > zField;
// The field-map dimensions
int nr, nz;
// The field map Length unit (string and number)
G4String lUnit;
double lenUnit;
// The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla)
G4String fUnit;
double fieUnit;
// The field-map Field normalisation factor
double fieldNormalisation;
// The physical limits of the defined region
double minimumr, maximumr, minimumz, maximumz;
// The physical extent of the defined region
double dr, dz;
// See the description under the class constructor
char fldType;
double ffieldValue;
void Invert(const char* indexToInvert);
};
#endif

76
geant4/LEMuSR/include/lem4TabulatedElementField3D.hh Normal file
View File

@ -0,0 +1,76 @@
#ifndef lem4TabulatedElementField3D_h
#define lem4TabulatedElementField3D_h 1
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
#include "globals.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
// Class for reading 3D electric and magnetic field map, either with or without coordinates.
class lem4TabulatedElementField3D : public F04ElementField
{
public: // with description
// Class constructor for 3D field map (x, y, z, EMx, EMy, EMz) - with EM = E or B
lem4TabulatedElementField3D(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
//
// "fieldType" is the type of EM field: electric - E, or magnetic - B
// "fieldValue" is the field to be applied (in T, or in kV/mm). The normalised field
// map values are multiplied by this value. The field-map itself has no units!
// "lenUnit" is the unit in which the grid coordinates of the field-map are specified
// "fieldNormalisation" is the normalisation factor that once applied to the tabulated field values
// satisfies the condition: (max. field value)*fieldNormalisation = 1
// To revert field direction, change its sign to negative.
// Virtual destructor
virtual ~lem4TabulatedElementField3D() {}
// addFieldValue() adds the field for THIS particular map into field[].
// point[] is expressed in GLOBAL coordinates.
void addFieldValue( const G4double Point[4], G4double* field) const;
// Usual Set and Get functions
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
// getWidth(), getHeight(), getLength(), return the dimensions of the field
// (used to define the boundary of the field)
virtual G4double getWidth() { return dx; } // x coordinate
virtual G4double getHeight() { return dy; } // y coordinate
virtual G4double getLength() { return dz; } // z coordinate
private:
// Storage space for the 3D table
std::vector< std::vector< std::vector< double > > > xField;
std::vector< std::vector< std::vector< double > > > yField;
std::vector< std::vector< std::vector< double > > > zField;
// The field-map dimensions
int nx, ny, nz;
// The field map Length unit (string and number)
///G4String lUnit;
char lUnit[50];
double lenUnit;
// The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla)
G4String fUnit;
double fieUnit;
// The field-map Field normalisation factor
double fieldNormalisation;
// The physical limits of the defined region
double minimumx, maximumx, minimumy, maximumy, minimumz, maximumz;
// The physical extent of the defined region
double dx, dy, dz;
// See the description under the class constructor
char fldType;
double ffieldValue;
void Invert(const char* indexToInvert);
};
#endif

52
geant4/LEMuSR/include/lem4TabulatedField2D.hh Normal file
View File

@ -0,0 +1,52 @@
#include "globals.hh"
#include "G4MagneticField.hh"
#include "G4ios.hh"
//cks For the special case of muon decay:
//#include "G4EventManager.hh"
//#include "G4RunManagerKernel.hh"
//csk
#include <fstream>
#include <vector>
#include <cmath>
class lem4TabulatedField2D
#ifndef STANDALONE
: public G4MagneticField
#endif
{
// Storage space for the table
std::vector< std::vector< double > > xField;
std::vector< std::vector< double > > zField;
// The dimensions of the table
int nx,nz;
// The physical limits of the defined region
double minx, maxx, minz, maxz;
// The physical extent of the defined region
double dx, dz;
double ffieldValue;
bool invertX, invertZ;
double positionOffset[4];
static lem4TabulatedField2D* pointerToTabulatedField2D;
public:
static lem4TabulatedField2D* GetInstance();
lem4TabulatedField2D(const char* filename, double fieldValue, double lenUnit, double fieldNormalisation );
// "lenUnit" is the unit in which the grid coordinates are specified in the table
// "fieldNormalisation" is the normalisation that has to be applied on the field values in the table
// such that the values correspond do 1T nominal value
// "fieldValue" is the field value (in T) that is required (i.e. values normalised to 1T will be
// multiplied by this value).
void GetFieldValue( const double Point[4],
double *Bfield ) const;
G4double GetFieldSetValue();
void SetFieldValue(double newFieldValue);
void SetPositionOffset(double offset[4]) { positionOffset[0]=offset[0]; positionOffset[1]=offset[1];
positionOffset[2]=offset[2]; positionOffset[3]=offset[3];}
};

46
geant4/LEMuSR/include/lem4TabulatedField3D.hh Normal file
View File

@ -0,0 +1,46 @@
#include "globals.hh"
#include "G4MagneticField.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
using namespace std;
class lem4TabulatedField3D
#ifndef STANDALONE
: public G4MagneticField
#endif
{
// Storage space for the table
vector< vector< vector< double > > > xField;
vector< vector< vector< double > > > yField;
vector< vector< vector< double > > > zField;
// The dimensions of the table
int nx,ny,nz;
// The physical limits of the defined region
double minx, maxx, miny, maxy, minz, maxz;
// The physical extent of the defined region
double dx, dy, dz;
double ffieldValue;
bool invertX, invertY, invertZ;
double positionOffset[4];
static lem4TabulatedField3D* pointerToTabulatedField3D;
public:
static lem4TabulatedField3D* GetInstance();
lem4TabulatedField3D(const char* filename, double fieldValue );
void GetFieldValue( const double Point[4],
double *Bfield ) const;
G4double GetFieldSetValue();
void SetFieldValue(double newFieldValue);
void SetPositionOffset(double offset[4]) { positionOffset[0]=offset[0]; positionOffset[1]=offset[1];
positionOffset[2]=offset[2]; positionOffset[3]=offset[3];}
};

43
geant4/LEMuSR/include/lem4UniformField.hh Normal file
View File

@ -0,0 +1,43 @@
#ifndef UNIFORM_BFIELD_HH
#define UNIFORM_BFIELD_HH
#include "G4LogicalVolume.hh"
#include "G4Box.hh"
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
// UniformField implements a constant electromagnetic field in any direction. TS
class lem4UniformField : public F04ElementField
{
public:
lem4UniformField(G4double EMF[6], G4LogicalVolume*, G4ThreeVector);
virtual ~lem4UniformField() {}
// TS: Define the two newly added VIRTUAL functions of F04ElementField
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
virtual G4double getLength() { return fieldLength; }
virtual G4double getWidth() { return fieldWidth; }
virtual G4double getHeight() { return fieldHeight; }
G4bool isOutside(G4ThreeVector& local) const;
G4bool isWithin (G4ThreeVector& local) const;
void addFieldValue(const G4double point[4], G4double field[6]) const;
private:
G4double EMfield[6];
G4double fieldLength;
G4double fieldWidth;
G4double fieldHeight;
};
#endif

26
geant4/LEMuSR/include/lem4VisManager.hh Normal file
View File

@ -0,0 +1,26 @@
#ifndef lem4VisManager_h
#define lem4VisManager_h 1
#ifdef G4VIS_USE
#include "G4VisManager.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class lem4VisManager: public G4VisManager {
public:
lem4VisManager ();
private:
void RegisterGraphicsSystems ();
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif
#endif

57
geant4/LEMuSR/include/musrTabulatedElementField2D.hh Normal file
View File

@ -0,0 +1,57 @@
#ifndef lem4TabulatedElementField2D_h
#define lem4TabulatedElementField2D_h 1
#include "globals.hh"
////#include "G4MagneticField.hh"
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
class lem4TabulatedElementField2D : public F04ElementField
{
public:
lem4TabulatedElementField2D(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
// "lenUnit" is the unit in which the grid coordinates are specified in the table
// "fieldNormalisation" is the normalisation that has to be applied on the field values in the table
// such that the values correspond do 1T nominal value
// "fieldValue" is the field value (in T) that is required (i.e. values normalised to 1T will be
// multiplied by this value).
/// Destructor.
virtual ~lem4TabulatedElementField2D() {}
/// addFieldValue() adds the field for this solenoid into field[].
/// point[] is in global coordinates.
void addFieldValue( const G4double Point[4], G4double* field) const;
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
// getWidth(), getHeight(), getLength(), return the dimensions of the field
// (used to define the boundary of the field)
virtual G4double getWidth() { return 2*dx; } // x coordinate
virtual G4double getHeight() { return 2*dx; } // y coordinate
virtual G4double getLength() { return 2*dz; } // z coordinate
private:
// Storage space for the table
std::vector< std::vector< double > > xField;
std::vector< std::vector< double > > zField;
// The dimensions of the table
int nx,nz;
// The physical limits of the defined region
double minimumx, maximumx, minimumz, maximumz;
// The physical extent of the defined region
double dx, dz;
double ffieldValue;
void Invert(const char* indexToInvert);
};
#endif

57
geant4/LEMuSR/include/musrTabulatedElementField3D.hh Normal file
View File

@ -0,0 +1,57 @@
#ifndef lem4TabulatedElementField3D_h
#define lem4TabulatedElementField3D_h 1
#include "globals.hh"
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
class lem4TabulatedElementField3D : public F04ElementField
{
public:
lem4TabulatedElementField3D(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter);
// "lenUnit" is the unit in which the grid coordinates are specified in the table
// "fieldNormalisation" is the normalisation that has to be applied on the field values in the table
// such that the values correspond do 1T nominal value
// "fieldValue" is the field value (in T) that is required (i.e. values normalised to 1T will be
// multiplied by this value).
/// Destructor.
virtual ~lem4TabulatedElementField3D() {}
/// addFieldValue() adds the field for this solenoid into field[].
/// point[] is in global coordinates.
void addFieldValue( const G4double Point[4], G4double* field) const;
G4double GetNominalFieldValue();
void SetNominalFieldValue(G4double newFieldValue);
// getWidth(), getHeight(), getLength(), return the dimensions of the field
// (used to define the boundary of the field)
virtual G4double getWidth() { return dx; } // x coordinate
virtual G4double getHeight() { return dy; } // y coordinate
virtual G4double getLength() { return dz; } // z coordinate
private:
// Storage space for the table
std::vector< std::vector< std::vector< double > > > xField;
std::vector< std::vector< std::vector< double > > > yField;
std::vector< std::vector< std::vector< double > > > zField;
// The dimensions of the table
int nx,ny,nz;
// The physical limits of the defined region
double minimumx, maximumx, minimumy, maximumy, minimumz, maximumz;
// The physical extent of the defined region
double dx, dy, dz;
double ffieldValue;
void Invert(const char* indexToInvert);
};
#endif

52
geant4/LEMuSR/include/musrTabulatedField2D.hh Normal file
View File

@ -0,0 +1,52 @@
#include "globals.hh"
#include "G4MagneticField.hh"
#include "G4ios.hh"
//cks For the special case of muon decay:
//#include "G4EventManager.hh"
//#include "G4RunManagerKernel.hh"
//csk
#include <fstream>
#include <vector>
#include <cmath>
class lem4TabulatedField2D
#ifndef STANDALONE
: public G4MagneticField
#endif
{
// Storage space for the table
std::vector< std::vector< double > > xField;
std::vector< std::vector< double > > zField;
// The dimensions of the table
int nx,nz;
// The physical limits of the defined region
double minx, maxx, minz, maxz;
// The physical extent of the defined region
double dx, dz;
double ffieldValue;
bool invertX, invertZ;
double positionOffset[4];
static lem4TabulatedField2D* pointerToTabulatedField2D;
public:
static lem4TabulatedField2D* GetInstance();
lem4TabulatedField2D(const char* filename, double fieldValue, double lenUnit, double fieldNormalisation );
// "lenUnit" is the unit in which the grid coordinates are specified in the table
// "fieldNormalisation" is the normalisation that has to be applied on the field values in the table
// such that the values correspond do 1T nominal value
// "fieldValue" is the field value (in T) that is required (i.e. values normalised to 1T will be
// multiplied by this value).
void GetFieldValue( const double Point[4],
double *Bfield ) const;
G4double GetFieldSetValue();
void SetFieldValue(double newFieldValue);
void SetPositionOffset(double offset[4]) { positionOffset[0]=offset[0]; positionOffset[1]=offset[1];
positionOffset[2]=offset[2]; positionOffset[3]=offset[3];}
};

46
geant4/LEMuSR/include/musrTabulatedField3D.hh Normal file
View File

@ -0,0 +1,46 @@
#include "globals.hh"
#include "G4MagneticField.hh"
#include "G4ios.hh"
#include <fstream>
#include <vector>
#include <cmath>
using namespace std;
class lem4TabulatedField3D
#ifndef STANDALONE
: public G4MagneticField
#endif
{
// Storage space for the table
vector< vector< vector< double > > > xField;
vector< vector< vector< double > > > yField;
vector< vector< vector< double > > > zField;
// The dimensions of the table
int nx,ny,nz;
// The physical limits of the defined region
double minx, maxx, miny, maxy, minz, maxz;
// The physical extent of the defined region
double dx, dy, dz;
double ffieldValue;
bool invertX, invertY, invertZ;
double positionOffset[4];
static lem4TabulatedField3D* pointerToTabulatedField3D;
public:
static lem4TabulatedField3D* GetInstance();
lem4TabulatedField3D(const char* filename, double fieldValue );
void GetFieldValue( const double Point[4],
double *Bfield ) const;
G4double GetFieldSetValue();
void SetFieldValue(double newFieldValue);
void SetPositionOffset(double offset[4]) { positionOffset[0]=offset[0]; positionOffset[1]=offset[1];
positionOffset[2]=offset[2]; positionOffset[3]=offset[3];}
};