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

11.2.2011 Kamil Sedlak

Browse Source 
1) correction of GNUmakefile - by mistake a musrSimTest was set there
   instead of musrSim
2) replacing odet_timeE variable by odet_timeMean
3) replacing one obsolete method from stepping action
   (Geant 4.9.4 was complaining) by a proper method
4) update of the first part of the documentation
5) implementing a possibility to write out output Root data
   into a different directory
6) perhaps some other minor changes
This commit is contained in:
sedlak 2011-02-11 16:45:24 +00:00
parent a5e99ed164
commit c2025fea97
9 changed files with 277 additions and 156 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
GNUmakefile
View File

@ -3,7 +3,7 @@
# GNUmakefile for examples module. Gabriele Cosmo, 06/04/98.
# --------------------------------------------------------------
name := musrSimTest
name := musrSim
G4TARGET := $(name)
G4EXLIB := true
##LDFLAGS := $(shell root-config --glibs)

BIN
doc/musrSim.pdf
View File

Binary file not shown.

291
doc/musrSim.tex
View File

@ -47,40 +47,62 @@ Geant4}. The root output variables are also described.
\section{Scope of the musrSim program}
The program ``musrSim'' is a relatively general program that can be used to simulate
the response of a $\mu$SR~\cite{Blundel:1999} instruments (detectors) to muons and their decay particles
(electrons, positrons and gammas). Even though musrSim is tailored to the needs of
the $\mu$SR technique~\cite{shirokaGeant}, it has been used also in the detector development
(electrons, positrons and gammas), optionaly including ``optical photons''.
Even though musrSim is tailored to the needs of
the $\mu$SR technique~\cite{shirokaGeant}, it has been used also
in the studies of beam-line elements like spin rotators, as well as
in the detector development
studies without any muons involved, e.g.\ to test the response of an APD-based
scintillator counters to the irradiation of Sr radioactive source~\cite{AlexeyTestAPD}.
It should be streitforward to apply musrSim also to the low energy particle-physics
experiments with muons.
The program is based on the Geant4~\cite{geant} and Root~\cite{root} libraries.
Geant4 is Monte Carlo toolkit used (not only) in particle physics to simulate
the passage of particles through the detectors.
Root is an analysis tool that allows to manipulate and analyse the simulated data,
the passage of particles through detectors.
Root is an analysis tool that allows one to manipulate and analyse the simulated data,
namely to plot histograms and other graphical output.
The simulation of an instrument consists of two steps -- simulation of the instrument
response, which is done by ``musrSim'' program, and the subsequent analysis of the
simulated output by either ``musrSimAna'' program (in case of \musr instruments) or
by Root (if a special analysis is required, e.g.\ for non-\musr apparatus).
The reason for this spit is purely
practical -- while it takes a lot of computer time to simulate large statistics data
with musrSim, the subsequent analysis of the simulated data is relatively quick
and allows the user to optimise and test different options of the analysis
(e.g.\ different thresholds in positron counters, different logic in connecting
veto and coincidence detectors, different rate of incoming muons, ...).
The ``musrSimAna'' program is described in a separate manual.
The aim of the musrSim is to provide an easy-to-use simulation program, which does
not require a deep knowledge of Geant4 in order to simulate a ($\mu$SR) detector.
In our view, the main advantages of musrSim are:
\begin{itemize}
\item Simple way how to define or modify the instrument geometry.
\item Limited (ideally no) need to modify and/or recompile the source code.
\item Simple way how to define or modify the instrument geometry, including
the sample environment, collimators and other parts.
\item Limited (ideally no) need to modify and/or recompile the source code,
because the parameters defining the instrument geometry, initial muon beam,
electromagnetic fields, and other parameters are defined in
a text file (the so called ``macro file'').
\item Implementation of the $\mu$SR-specific classes (muon spin rotation
in magnetic fields, muonium formation and decay, ...).
\item Possibility to read in the output files of the TURTLE~\cite{turtle}
program for the beam-line simulation.
\item Simple way how to define (overlapping) electromagnetic fields.
\item Output in the Root tree.
\item Possibility to analyse the output with a general ``musrSimAna'' program.
\item Possibility to use musrSim easily for calculating muon stopping profile
(also in e.g.\ sample cells) or in developments of detector components
(e.g.\ light propagation in the scintillator of a positron/muon counter
and the subsequent light collectin in a photomultiplier tube or APD).
\end{itemize}
%
On the other hand, there are also some drawbacks and limitations:
\begin{itemize}
\item The user has to have an installation of Geant4 and Root before installing musrSim.
\item It is supposed the user will analyse the data with Root, therefore
Root has to be installed and some knowledge of it is needed.
Even though is relatively easy to simulate a $\mu$SR instrument,
and to create the output Root file without even knowing the c++
programming language, some c++ programming is needed to
analyse the simulation output and to plot graphs.
Root has to be installed.
\item At present time the program does not simulate any muon-spin related
physics processes happening in the sample, except for the muon
spin rotation.
@ -95,17 +117,25 @@ On the other hand, there are also some drawbacks and limitations:
%=============================================================================================
\section{How to install and run musrSim}
To install and run musrSim, one has to install Geant4 and Root first.
The present version of musrSim has been tested with Geant version 4.9.1, patch no.~3
and Root version 5.20.00.
The present version of musrSim has been tested with Geant version 4.9.4,
and with Root version 5.24.00. While the version of Root should not be critical,
the users of musrSim are encouraged to always use the latest version of Geant4
due to continuous improvements of this package. A novice
user of Geant4 may consider reading 30 pages of chapter~2 of the
Geant4 User's Guide for Application Developers, called
``Getting Started with Geant4 - Running a Simple Example''.
Once Geant4 has been successfully installed and some of the default Geant4 exmples
Once Geant4 has been successfully installed and some of the default Geant4 examples
has been run, the musrSim installation package can be downloaded from the web page
http://lmu.web.psi.ch/simulation/index.html.
Usually the ``env.sh'' script has to be run to set-up the environment variables
appropriatelly before the musrSim or any other Geant4 application can be compiled
or run.
The simulation can be executed by \emph{``musrSim RUNNUMBER.mac''}, where RUNNUMBER.mac is a ``macro file''
containing the information about the instrument setup. The string ``RUNNUMBER'' represents the integer run number.
appropriatelly before the musrSim (and/or any other Geant4 application) can be compiled
and run.
The simulation is started by executing:
{\bf $>$ musrSim \emph{RUNNUMBER}.mac} \\
where \emph{RUNNUMBER}.mac is a ``macro file'' containing the information about
the instrument setup. The string \emph{RUNNUMBER} is an integer representing the run number.
In order to simulate a new instrument, the user has to define the following blocks of information
in the macro file:
@ -113,10 +143,12 @@ in the macro file:
\begin{itemize}
\item Define the geometry (the so-called ``volumes'') of the new instrument.
Note that in Geant4 volumes can be included inside other volumes
(a ``daughter'' volume is positioned inside the ``mother'' volume),
(a ``daughter'' volume is positioned inside its ``mother'' volume),
and it is therefore necessary to distinguish between the global (world)
coordinates and the coordinates of the daughter volumes (local coordinates).
It is not allowed to overlap any two different volumes partially.
It is not allowed to overlap the volumes (with the exception of the
daughter and mother volume, in which case the daughter volume has to be fully
contained within its mother volume).
\item Define the electric and magnetic fields.
\item Define physics processes relevant for your case.
\item Define the initial muon parameters (more generally -- initial particle parameters).
@ -126,10 +158,28 @@ in the macro file:
\end{itemize}
By default, the output of the simulation is written out in the subdirectory ``data'' with
the name ``musr\_RUNNUMBER.root''. The default ``data'' directory can be changed by the
command ``/musr/command rootOutputDirectoryName \emph{dirName}''.
(Note that the execution of the simulation can be terminated gently by creating a file ``RUNNUMBER.stop'' in the working directory.)
the name ``musr\_RUNNUMBER.root''. The default ``data'' subdirectory can be changed
(see Chapter~\ref{sec:otherParameters}). The simulated data are storred in a Root tree,
which is some kind of a table with many variables for every simulated event, inside
the musr\_RUNNUMBER.root file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{How to stop musrSim}
The execution of musrSim stops as soon as one of the following condition is fulfilled:
\begin{itemize}
\item Number of simulated events reaches the number of required events defined by
command ``{\bf /run/beamOn \emph{nrOfEvents}}'' (see Chapter~\ref{sec:otherParameters}).
\item Time for which the simulation is running exceeds the maximum allowed time defined
by command ``{\bf /musr/command maximumRunTimeAllowed \emph{timeMax}}''
(see Chapter~\ref{sec:otherParameters}).
\item Within a few seconds after a file ``RUNNUMBER.stop'' is created in the
directory, in which the musrSim program was started. This allows the user
to stop the simulation gently at any time.
\end{itemize}
If the musrSim program is killed in any other way, the output Root file does not close properly, and
subsequently can not be analysed. The simulated data are lost.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conventions}
The default units of the musrSim in both the macro file (RUNNUMBER.mac) and in the Root tree
are summarised in table~\ref{tab:units}.
@ -140,36 +190,61 @@ are summarised in table~\ref{tab:units}.
\lower 1mm \hbox{\textbf{Quantity}} && \lower 1mm \hbox{\textbf{Default unit}} \\[5pt]
\hline
Length && mm \\
Time && $\mu$s \\
Time && $\mu$s (sometimes ns) \\
Energy && MeV \\
Momentum && MeV/c\\
Magnetic field && tesla \\
Electric field && keV/mm \\
Anlgel && degree \\
\hline
\end{tabular}
\caption{The default units in musrSim.}
\caption{The default units of musrSim.}
\label{tab:units}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Tips and tricks}
\begin{itemize}
\item Visualise your instrument geometry during its construction.
\item Check the output file for error messages, especially for volume overlaps.
\item Note the the dimensions in volume definitions are often half-lenghts, not
the lenghts (e.g.\ half-lengths of the box edges).
\item The order of some commands in macro file matters -- e.g.\ one has to define
a mother volume before the daughter volume, etc.
\item Threre are some special volume names, namely \emph{World},
\emph{Target} (same as \emph{target}), \emph{M0}, \emph{M1}, \emph{M2} and
volumes starging with keywords \emph{Shield} (same as \emph{shield}), and
volumes containing string \emph{save}. All these volume influence the
behaviour of the simulation, so understand how they act (see different chapters
of this manual) before you use them.
\item When implementing a new magnetic or electric field always check that fields at
a few random space points correspond to the field you expect to observe
(use command /musr/command globalfield printFieldValueAtPoint~).
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Detector construction}
The user must first define the instrument geometry in the macro file. It should be realively
easy to understand how this is done from the example macro files distributed with musrSim
program.
An important parameter to note in the following description of ``/musr/command construct''
command is the \emph{idNumber} -- an integer number uniquely identifying every volume.
The ID volume numbers are later on used also in the musrSimAna program, when some identification
of a volume becomes necessary (i.e.\ volumes are described by \emph{idNumber} rather than
by their \emph{name} command later on).
Another crutial parameter in ``/musr/command construct'' command
is \emph{sensitiveClass}, which defines whether the volume
is sensitive (i.e.\ a signal can be detected in the volume) or not (i.e.\ the volume
is just a dead material influencing the penetrating particles but not detecting them).
\begin{description}
\item{\bf /musr/command rotation \emph{matrixName} $\alpha$ $\beta$ $\gamma$} \\
{\bf /musr/command rotation \emph{matrixName} \emph{vx} \emph{vy} \emph{vz} \emph{angle}}\\
These commands define a rotation matrix of the name \emph{matrixName} that can be used later on
during the definition of the detector geometry (see command ``/musr/command construct'').
It can be defined either by the Euler angles (if there are three float parameters behind the
\emph{matrixName}) or by the vector \emph{(vx,vy,vz)} and an \emph{angle} of rotation around this
vector (if the fourth float parameter behind the \emph{matrixName} is non-zero).
All angles are specified in degrees.
\item{\bf /musr/command construct \emph{solid} \emph{name} \emph{dimensions} ... \emph{material}
\emph{x} \emph{y} \emph{z} \emph{motherVolume} \emph{matrixName}
\emph{sensitiveClass} \emph{idNumber} }\\
This command defines a volume in {\sc Geant4} (It comprises three steps of {\sc Geant4}: defines a solid,
logical volume and physical volume. Details can to be found in {\sc Geant4} manual). \\
logical volume and physical volume. Details can be found in {\sc Geant4} manual). \\
\begin{itemize}
\item \emph{solid} (string) can be one of the G4VSolid.cc particular types, presently ``tubs'', ``cons'',
``box'', ``trd'', ``sphere'', ``para'',
@ -178,12 +253,15 @@ are summarised in table~\ref{tab:units}.
(a tube with a rectangular hole along its axis), "tubsboxsegm"
(a volume that looks like an intersection of tube and box) and
``trd90y'' (a trd volume rotated by 90 degrees around $y$ axis in addition
to the rotation requested by \emph{matrixName}). Not all G4VSolids are
to the rotation requested by \emph{matrixName}),
``cylpart'' (cylinder from which a box has been subtracted)
``GPDcollimator'' (collimator used at GPD instrument). Not all G4VSolids defined
in Geant4 are
presently supported, but it is relatively easy to implement a new kind of solids
in the musrDetectorConstruction.cc class.
\item \emph{name} (string) stands for the name of the volume. As the command
``/musr/command construct'' constructs
three kinds of classes/volumes (the solid, logical volume and physical
three kinds of Geant4 classes/volumes (the solid, logical volume and physical
volume), there are three names of the concrete volume used internally inside
musrSim: sol\_\emph{name}, log\_\emph{name} and phys\_\emph{name}.
The main volume, inside which all other volumes are positioned, has to be called ``World''.
@ -198,17 +276,17 @@ are summarised in table~\ref{tab:units}.
``G4\_PLASTIC\_SC\_VINYLTOLUENE'' for a scintillator, ...).
One can also define a new material inside the function
musrDetectorConstruction::DefineMaterials(). Presently ``Mylar'', ``Brass''
``Steel'', ``Macor'', ``MCPglass'', ``MgO'', ``SiO2'', ``K2O'' and ``B2O3'' are defined there.
``Steel'', ``Macor'', ``MCPglass'', ``MgO'', ``SiO2'', ``Al2O3'', ``K2O'' and ``B2O3'' are defined there.
\item \emph{x, y, z} (floats) -- coordinates of the volume, used to position the volume within
its mother volume (as used by the G4PVPlacement).
Thus these coordinates are interpreted in the local coordinate system of the \emph{motherVolume}.
\item \emph{motherVolume} (string) -- name of the mother volume, in which the given volume should be
positioned. Note that the mother volume has to be defined first (before its
daughter), and that the name of mother starts with a string {\bf log\_}\emph{name},
daughter), and that the name of mother starts with a string {\bf log\_}\emph{mothername},
following the naming convention defined above.
When the ``World'' volume is defined, its \emph{motherVolume} should be set to ``no\_logical\_volume''.
\item \emph{matrixName} (string) -- name of the rotation matrix that will be used to position
the volume inside its mother volume (as used in member function G4PVPlacement).
the volume inside its mother volume (as used by the G4PVPlacement).
Use string ``norot'' if no rotation is required for the given volume.
Otherwise the rotation matrix has to be defined by the command line
``/musr/command rotation'' {\bf before} the given volume is defined.
@ -227,6 +305,15 @@ are summarised in table~\ref{tab:units}.
(the variable ``det\_ID[det\_n]'').
\end{itemize}
\item{\bf /musr/command rotation \emph{matrixName} $\alpha$ $\beta$ $\gamma$} \\
{\bf /musr/command rotation \emph{matrixName} \emph{vx} \emph{vy} \emph{vz} \emph{angle}}\\
These commands define a rotation matrix of the name \emph{matrixName} that can be used
during the definition of the detector geometry (see the command ``/musr/command construct'').
It can be defined either by the Euler angles (if there are three float parameters behind the
\emph{matrixName}) or by the vector \emph{(vx,vy,vz)} and an \emph{angle} of rotation around this
vector (if the fourth float parameter behind the \emph{matrixName} is non-zero).
All angles are specified in degrees.
\item{\bf /musr/command region define \emph{regionName} \emph{logicalVolume}}\\
The ``G4Region'' can be created using this command, and a logical volume of the
name \emph{logicalVolume} will be assigned to it. If the G4Region of the name
@ -242,7 +329,8 @@ are summarised in table~\ref{tab:units}.
Set the so-called ``production cuts'' in the G4Region called \emph{regionName}.
The variables \emph{gammaCut, electronCut} and \emph{positronCut} are given in mm.
\end{description}
Three special volumes ``Target, M0, M1 and M2''.
{\huge !!! EDIT: !!! Three special volumes ``Target, M0, M1 and M2''.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Electric and magnetic fields}
@ -255,7 +343,7 @@ Three special volumes ``Target, M0, M1 and M2''.
\emph{[fieldValueFinal]} \emph{[fieldNrOfSteps]}} \\
This command specifies the electric and/or magnetic fields, which are (in some sense)
independent of any logical volume and can overlap with each other.
In the case of tabulated field read in from and external field map file the
In the case of tabulated field read in from an external field map file, the
field values used internally by the Geant4 are linearly interpolated using
eight (3D) or four (2D) grid points surrounding the point of interest.
%
@ -263,15 +351,15 @@ Three special volumes ``Target, M0, M1 and M2''.
\item \emph{fieldName} (string) -- name of the field (important mainly for the user and
print-out messages of the musrSim.
\item \emph{half\_x}, \emph{half\_y}, \emph{half\_z} (floats) -- the (half) dimensions
of the box, within which the uniform field is defined.
of the box, within which the uniform field is defined. Used only for the uniform fields.
\item {\bf uniform / fromfile} -- specifies whether the field is uniform within
some volume or whether it is read in from an external file as a field-map.
\item \emph{X}, \emph{Y}, \emph{Z} (floats) -- position of the centre of the field
in the {\bf global} coordinate system. IMPORTANT: For some technical internal
Geant4 reasons, this POSITION HAS TO LAY WITHIN THE \emph{logicalVolume}!
(Note that the logical volume may be positioned somewhere deep in a volume
structure, not directly within the ``World'' volume, and
therefore the (local) coordinates in the definition of the the logical volume
Geant4 reasons, this \bf{POSITION HAS TO LAY WITHIN THE \emph{logicalVolume}!}
Note that the logical volume may be positioned somewhere deep in a volume
structure, not directly within the ``World'' volume, and
therefore the (local) coordinates in the definition of the logical volume
do not have to match the (global) coordinates \emph{X}, \emph{Y} and \emph{Z}.
\item \emph{logicalVolume} (string) -- specifies the logical volume, to which
the field is ``assigned''. One may ask, why a logical volume is needed for a field
@ -289,15 +377,17 @@ Three special volumes ``Target, M0, M1 and M2''.
larger volume. The volume can also be made of vacuum (i.e.\ G4\_Galactic).
\item \emph{Bx}, \emph{By}, \emph{Bz}, \emph{Ex}, \emph{Ey}, \emph{Ez} (float) -- the vector
of the uniform electromagnetic field. The units are tesla, for the first three
components, and kilovolt/mm for the last three components.
components, and kilovolt/mm for the last three components.
Used only for the uniform fields.
\item \emph{fieldTableType} (string) -- specifies the format in which the field map is
written in the file. In general, the field is specified in a grid of
three space coordinates $x$, $y$ and $z$ (3D).
Sometimes it is convenient to use the symmetry of the field
and to reduce the field description to $R$ and $z$ (2D) only.
In the following, we use this terms: \\
\emph{nx, ny, nz} or \emph{nR, nz} -- the number of divisions of the grid in
\emph{x, y, z} or \emph{R, z}.
In the following, we use the following terms: \\
\emph{nx, ny, nz} -- the number of divisions of the grid in
\emph{x, y, z}.\\
\emph{nR, nz} -- the number of divisions of the grid in \emph{R, z}.\\
\emph{length unit} -- the unit in which the grid coordinates are specified, usually
cm or m.\\
\emph{field normalisation factor} -- a multiplicative factor applied to the values
@ -329,9 +419,9 @@ Three special volumes ``Target, M0, M1 and M2''.
However, a different \emph{field normalisation factor} can be specified
in the field map file using the keyword ``fieldNormalisation \emph{number}''
before the line started with 0.\\
It is expected that the we first loop over the $z$ coordinate of the field map,
then (when $z$ changed from minimum to maximum) it is looped over $y$ coordinate,
and the highest-lever loop goes over $x$ coordinate. Hoever, if the order
It is expected that the first loop goes over the $z$ coordinate of the field map,
then (after $z$ changed from minimum to maximum) it is looped over the $y$ coordinate,
and the highest-level loop goes over $x$ coordinate. However, if the order
of looping is reversed in the field map, it can be specified using the
keyword ``variableIncreasingOrder xyz'' placed in the field map before
the line started with 0.\\
@ -356,7 +446,7 @@ Three special volumes ``Target, M0, M1 and M2''.
Similar case is the ``symmetryType 2'', where the planes of symmetry are
(x,y) and (y,z).
These two symmetry types are realised in a the spin rotator oriented along the z axis.\\
Example of the beginning of the field map file:\\
Example of the beginning of the 3DBOpera-type field map file:\\
2 2 55\\
1 X\\
2 Y\\
@ -396,7 +486,7 @@ Three special volumes ``Target, M0, M1 and M2''.
\emph{R, z, dummy, Field\_R, Field\_z, dummy}\\
\item \emph{fieldInputFileName} (string) -- Name of the field map file.
\item \emph{fieldValue} (float) -- the value of the field at some reference point
(usually in the centre of the field). It serves as some multiplicative
(usually in the centre of the field). It acts as an multiplicative
factor. The units are tesla for the magnetic field and kV/mm
for the electric field.
\item \emph{[fieldValueFinal]} and \emph{[fieldNrOfSteps]} (floats)
@ -415,7 +505,7 @@ Three special volumes ``Target, M0, M1 and M2''.
See ``musrDetectorConstruction.cc'' and ``BLEngeFunction.hh'' for the details.
\item{\bf /musr/command globalfield setparameter \emph{parameterName} \emph{parameterValue} }\\
Set up some parameters used internally by Geant4 when calculating the motion
Set up some accuracy parameters used internally by Geant4 when calculating the motion
of charged particles in the magnetic field.\\
\emph{parameterName} (string) -- one of the following parameters: ``SetDeltaIntersection''
``SetDeltaOneStep'', ``SetMinimumEpsilonStep'', ``SetMaximumEpsilonStep'',
@ -444,24 +534,34 @@ Three special volumes ``Target, M0, M1 and M2''.
\emph{particle} (string) -- name of the particle to which a process is applied.\\
\emph{process} (string) -- name of the process to be assigned.\\
\emph{ordAtRestDoIt, ordAlongSteptDoIt, ordPostStepDoIt} (int) -- priority switches.\\
See the file musrPhysicsList.cc for the list of defined processes (e.g. G4MultipleScattering,
See the file musrPhysicsList.cc for the list of defined processes (e.g. G4eMultipleScattering,
G4eIonisation, ...) and Geant4 manual for the detail description of the processes.
There is one special process, combined from G4MultipleScattering and G4CoulombScattering,
defined by the following command:\\
{\bf /musr/command process addProcess \emph{particle} MultipleAndCoulombScattering \emph{ordAtRestDoIt} \emph{ordAlongSteptDoIt} \emph{ordPostStepDoIt} \emph{G4Region1} [\emph{G4Region2}] [\emph{G4Region3}]}\\
The G4MultipleScattering (rough but very fast approximation of scattering) will be applied
elsewhere in the detector, except for the \emph{G4Region1} (and eventually \emph{G4Region2}
and \emph{G4Region3}), where more precise but very slow process G4CoulombScattering
will be applied instead of G4MultipleScattering. Note that up to three
G4Regions are supported at the moment, but this limitation is not intrinsic to {\sc Geant4}
and it can be therefore changed in musrPhysicsList.cc, if needed. The G4Regions have to
be defined in the detector construction phase by the command ``/musr/command region define ...''.
% There is one special process, combined from G4MultipleScattering and G4CoulombScattering,
% defined by the following command:\\
%{\bf /musr/command process addProcess \emph{particle} MultipleAndCoulombScattering \emph{ordAtRestDoIt} \emph{ordAlongSteptDoIt} \emph{ordPostStepDoIt} \emph{G4Region1} [\emph{G4Region2}] [\emph{G4Region3}]}\\
% The G4MultipleScattering (rough but very fast approximation of scattering) will be applied
% elsewhere in the detector, except for the \emph{G4Region1} (and eventually \emph{G4Region2}
% and \emph{G4Region3}), where more precise but very slow process G4CoulombScattering
% will be applied instead of G4MultipleScattering. Note that up to three
% G4Regions are supported at the moment, but this limitation is not intrinsic to {\sc Geant4}
% and it can be therefore changed in musrPhysicsList.cc, if needed. The G4Regions have to
% be defined in the detector construction phase by the command ``/musr/command region define ...''.
\end{description}
%=============================================================================================
\section{Initial (muon) beam parameters}
%
In historical versions of musrSim
it had been implicitly assumed that the (muon) beam was oriented along the $z$-axis. Many
variables defined below are therefore related to this $z$-axis (e.g.\ beam tilt, vertex sigma, ...).
However, later on we implemented a possibility to change the direction of the initial beam
to a random vector by the command ``/gun/direction''. It now works like this:
the primary particles are first generated using all parameters of /gun/* commands
as if the beam went along the $z$-axis, and just in the last moment before Geant4 starts
to track them, they are (optionaly) rotated to the direction defined by the /gun/direction command.
This way the smearing of the vertex as well as beam tilt/pitch are propagated through the rotation.
\begin{description}
\item{\bf /gun/primaryparticle \emph{primaryParticleName}}\\
@ -472,20 +572,33 @@ Three special volumes ``Target, M0, M1 and M2''.
\item{\bf /gun/meanarrivaltime \emph{meanArrivalTime}}\\
(default: /gun/meanarrivaltime 33.33333 microsecond)\\
Set mean arrival time difference between two subsequent muons (at continuos beam).
The output variable ``timeToNextEvent'' is subsequently generated using
The output variable ``timeToNextEvent'' is subsequently randomly generated using
the value of \emph{meanArrivalTime} and filled into the Root tree.
\item{\bf /gun/starttime \emph{t0} \emph{unit}}\\
By default, muons are generated at time = 0. The time of generation of muons can be
set randomly according to the Gaussian or uniform distribution using variables
``/gun/starttime'' and ``/gun/starttimesigma''. See the description of the
``/gun/vertexsigma'' to understand how the choice is done. This command appears
useful for the simulation of the laser-induced muon beam at ISIS,
however it might be confusing or misleading when used together with
the /gun/meanarrivaltime command. It is therefore recommended not to use it
unless necessary.
\item{\bf /gun/starttimesigma \emph{t0} \emph{unit}}\\
See the description of ``/gun/starttime'' command.
\item{\bf /gun/vertex \emph{x0} \emph{y0} \emph{z0} \emph{unit}}\\
(default: /gun/vertex 0 0 -100 mm) \\
Set mean values of the $x$, $y$ and $z$ coordinates of the generated particles (muons).
The smearing around these mean values instead is set by /gun/vertexsigma and
restricted by /gun/vertexboundary (see below).\\
The smearing around these mean values can be set by /gun/vertexsigma and
restricted by /gun/vertexboundary and/or /gun/boxboundary (see below).\\
(Applicable also to TURTLE input).
\item{\bf /gun/vertexsigma \emph{xSigma} \emph{ySigma} \emph{zSigma} \emph{unit}}\\
(default: /gun/vertexsigma 0 0 0 mm) \\
If {\it xSigma} $>0$ ... set $\sigma$, i.e. the standard deviation (RMS), of the $x$ coordinate
of the generated particles (muons) to $\sigma=xSigma$. The $x$ coordinate of the initial
If {\it xSigma} $>0$ ... set $\sigma_x$, i.e. the standard deviation (RMS), of the $x$ coordinate
of the generated particles (muons) to $\sigma_x=xSigma$. The $x$ coordinate of the initial
muon is then generated according to the Gaussian distribution with the mean value of $x0$
and the standard deviation of {\it xSigma}. \\
If {\it xSigma} $<0$ ... the $x$ coordinate of the initial
@ -494,15 +607,6 @@ Three special volumes ``Target, M0, M1 and M2''.
Similar is true for {\it ySigma} and {\it zSigma}. \\
(Ignored by the TURTLE input).
\item{\bf /gun/starttime \emph{t0} \emph{unit}}\\
By default, muons are generated at time = 0. The time of generation of muons can be
set randomly according to the Gaussian or uniform distribution using variables
``/gun/starttime'' and ``/gun/starttimesigma''. See the description of the
``/gun/vertexsigma'' to understand how the choice is done.
\item{\bf /gun/starttimesigma \emph{t0} \emph{unit}}\\
See the description of ``/gun/starttime'' command.
\item{\bf /gun/vertexboundary \emph{R\_max} \emph{z\_min} \emph{z\_max} \emph{unit}}\\
Set maximum allowed radius, and minimum and maximum z coordinate of the generated particles (muons).
This command might be useful especially if the user wants to restrict the
@ -528,7 +632,7 @@ Three special volumes ``Target, M0, M1 and M2''.
The variables \emph{xMaxSource0}, \emph{yMaxSource0} and \emph{xMaxSource0} define the
centre of the box, the variables \emph{xMaxSource}, \emph{yMaxSource} and \emph{xMaxSource} define the
halfwidth, halfheight and halflength of the box. \\
This command can be useful for muonium excitations by laser (low energy muon beam induced by laser).\\
This command can be useful for the laser-induced muon beam at ISIS.\\
(Ignored by the TURTLE input).
\item{\bf /gun/kenergy \emph{kineticEnergy} \emph{unit}}\\
@ -599,6 +703,13 @@ Three special volumes ``Target, M0, M1 and M2''.
``/gun/muonPolarizVector'' is smaller than 1e-8!} \\
(Applicable also to TURTLE input).
\item{\bf /gun/direction \emph{xDirection} \emph{yDirection} \emph{zDirection}}\\
Set initial beam direction as a vector (without units).
The vector does not have to be normalised to 1. Both beam direction
and beam spot are rotated.
See comments at the beginning of this chapter.\\
(Applicable also to TURTLE input).
\item{\bf /gun/decaytimelimits \emph{muDecayTimeMin} \emph{muDecayTimeMax} \emph{muMeanLife} \emph{unit}}\\
(default: /gun/decaytimelimits $-1$ $-1$ 2197.03\,ns ) \\
If {\it muDecayTimeMax} is less or equal to zero, this command is ignored,
@ -617,16 +728,16 @@ Three special volumes ``Target, M0, M1 and M2''.
more events than saved in the TURTLE file), the TURTLE file be be rewind to its
beginning. Note that this does not mean that the same events will be simulated
after the rewind, because the random seed will be set differently than at the
beginning of the simulation. Note that the muons initialised
beginning of the simulation. The muons initialised
at the same position and with the same momentum will have completely different
(random) multiple scattering, penetration depths, decay times,
decay positron energies and angles, ..., and therefore will be (almost completely)
different events not affecting the statistical quality of the sample.
different events not affecting the statistical quality of the simulation.
\item{\bf /gun/turtleZ0position \emph{z0\_InitialTurtle} \emph{unit}}\\
Set the z-position which has been used to generate the TURTLE file.\\
If this value differs from the $z0$ value of the ``/gun/vertex'' command,
than the particle initial position is extrapolated from $z0\_InitialTurtle$
than the particle initial position is linearly extrapolated from $z0\_InitialTurtle$
to the point corresponding to $z0$, using the direction of its momenta.\\
MORE DETAILS:\\
When running TURTLE (e.g. when generating the TURTLE file using the TURTLE program),
@ -674,7 +785,8 @@ Three special volumes ``Target, M0, M1 and M2''.
the decays of radioactive atoms and for other purposes.
Whenever the /gps/ keyword is used, the ``G4ParticleGun'' is not initiated
(and all /gun/* commands are ignored).
The description of GPS can be found on the web, some of the useful commands are:\\
The description of GPS can be found at http://reat.space.qinetiq.com/gps/~,
some of the useful commands are:\\
/gps/particle ion\\
/gps/ion 38 90 0 0\\
/gps/position 0 0 0\\
@ -804,7 +916,7 @@ in order to simulate optical photons:
\end{description}
\subsection*{Tips and tricks}
\subsection*{Tips and tricks for optical photons}
\begin{itemize}
\item One has to assign a non-zero EFFICIENCY and a REFLECTIVITY smaller than 1 to a boundary surface
between the scintillator and sensitive device (e.g.\ an APD).
@ -816,6 +928,7 @@ in order to simulate optical photons:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Some other parameters}
\label{sec:otherParameters}
%
\begin{description}
\item{\bf /run/beamOn \emph{nrOfEvents}}\\

6
include/musrRootOutput.hh
View File

@ -80,7 +80,7 @@ class musrRootOutput {
G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex, G4int particleID) ;
void SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA,
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast,
G4double timeB, G4double timeC, G4double timeD, G4double timeMean, G4double timeLast,
G4double timeCFD, G4double amplCFD);
void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z, G4double time,
@ -221,7 +221,7 @@ class musrRootOutput {
static G4bool store_odet_timeB;
static G4bool store_odet_timeC;
static G4bool store_odet_timeD;
static G4bool store_odet_timeE;
static G4bool store_odet_timeMean;
static G4bool store_odet_timeLast;
static G4bool store_odet_timeCFD;
static G4bool store_odet_amplCFD;
@ -329,7 +329,7 @@ class musrRootOutput {
G4double odet_timeB[odet_nMax];
G4double odet_timeC[odet_nMax];
G4double odet_timeD[odet_nMax];
G4double odet_timeE[odet_nMax];
G4double odet_timeMean[odet_nMax];
G4double odet_timeLast[odet_nMax];
G4double odet_timeCFD[odet_nMax];
G4double odet_amplCFD[odet_nMax];

11
include/musrScintSD.hh
View File

@ -78,12 +78,14 @@ class musrScintSD : public G4VSensitiveDetector
// Optical Photon Signal Analysis (OPSA)
void Set_OPSA_minNrOfDetectedPhotons(G4int val) {OPSA_minNrOfDetectedPhotons=val;}
void Set_OPSA_SignalSeparationTime(G4double val) {OPSA_signalSeparationTime=val;}
void Set_OPSA_frac(G4double a, G4double b, G4double c, G4double d, G4double e)
{OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=d; OPSA_fracE=e;}
void Set_OPSA_frac(G4double a, G4double b, G4double c, G4double d)
{OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=d;}
void Set_OPSA_CFD(G4double a1, G4double delay, G4double timeShiftOffset)
{OPSA_CFD_a1=a1; OPSA_CFD_delay=delay; OPSA_CFD_timeShiftOffset = timeShiftOffset;}
void AddEventIDToMultimapOfEventIDsForOPSAhistos (G4int ev_ID, G4int detector_ID) {
bool_multimapOfEventIDsForOPSAhistosEXISTS=true;
bool_multimapOfEventIDsForOPSAhistosEXISTS=true;
if (ev_ID==-1) bool_StoreThisOPSAhistSUMMED = true;
if (ev_ID==-2) bool_StoreThisOPSAhistALL = true;
multimapOfEventIDsForOPSAhistos.insert(std::pair<G4int,G4int>(ev_ID,detector_ID));
}
void SetOPSAhistoBinning(Int_t nBins, Double_t min, Double_t max) {
@ -115,11 +117,12 @@ class musrScintSD : public G4VSensitiveDetector
G4double OPSA_fracB;
G4double OPSA_fracC;
G4double OPSA_fracD;
G4double OPSA_fracE;
typedef std::multimap<G4int,signalInfo*> OPSA_signal_MapType;
OPSA_signal_MapType OPSA_signal_Map;
G4bool bool_multimapOfEventIDsForOPSAhistosEXISTS;
G4bool bool_StoreThisOPSAhistSUMMED;
G4bool bool_StoreThisOPSAhistALL;
typedef std::multimap<G4int,G4int> multimapOfEventIDsForOPSAhistos_Type;
multimapOfEventIDsForOPSAhistos_Type multimapOfEventIDsForOPSAhistos;
TH1D* OPSAhisto;

8
src/musrDetectorConstruction.cc
View File

@ -711,9 +711,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
myMusrScintSD -> Set_OPSA_SignalSeparationTime(fVarValue*nanosecond);
}
else if (strcmp(varName,"photonFractions")==0) {
double a, b, c, d, e;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %lf %lf",&a, &b, &c, &d, &e);
myMusrScintSD -> Set_OPSA_frac(a,b,c,d,e);
double a, b, c, d;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %lf",&a, &b, &c, &d);
myMusrScintSD -> Set_OPSA_frac(a,b,c,d);
}
else if (strcmp(varName,"eventsForOPSAhistos")==0) {
int i_eventID, i_detectorID;
@ -1169,7 +1169,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"odet_timeB")==0) {musrRootOutput::store_odet_timeB = false;}
if (strcmp(tmpString2,"odet_timeC")==0) {musrRootOutput::store_odet_timeC = false;}
if (strcmp(tmpString2,"odet_timeD")==0) {musrRootOutput::store_odet_timeD = false;}
if (strcmp(tmpString2,"odet_timeE")==0) {musrRootOutput::store_odet_timeE = false;}
if (strcmp(tmpString2,"odet_timeMean")==0) {musrRootOutput::store_odet_timeMean = false;}
if (strcmp(tmpString2,"odet_timeLast")==0) {musrRootOutput::store_odet_timeLast = false;}
}
else if(strcmp(storeIt,"on")==0) {

15
src/musrRootOutput.cc
View File

@ -156,7 +156,7 @@ G4bool musrRootOutput::store_odet_timeA = true;
G4bool musrRootOutput::store_odet_timeB = true;
G4bool musrRootOutput::store_odet_timeC = true;
G4bool musrRootOutput::store_odet_timeD = true;
G4bool musrRootOutput::store_odet_timeE = true;
G4bool musrRootOutput::store_odet_timeMean = true;
G4bool musrRootOutput::store_odet_timeLast = true;
G4bool musrRootOutput::store_odet_timeCFD = true;
G4bool musrRootOutput::store_odet_amplCFD = true;
@ -173,6 +173,11 @@ void musrRootOutput::BeginOfRunAction() {
// sprintf(RootOutputFileName, "data/musr_%i.root", tmpRunNr);
sprintf(RootOutputFileName, "%s/musr_%i.root",rootOutputDirectoryName,tmpRunNr);
rootFile=new TFile(RootOutputFileName,"recreate");
if (rootFile->IsZombie()) {
char message[200];
sprintf(message,"musrRootOutput::BeginOfRunAction() Root output file %s can not be created",RootOutputFileName);
musrErrorMessage::GetInstance()->musrError(FATAL,message,false);
}
rootTree=new TTree("t1","a simple Tree with simple variables");
if (store_runID) {rootTree->Branch("runID",&runID,"runID/I");}
if (store_eventID) {rootTree->Branch("eventID",&eventID,"eventID/I");}
@ -283,7 +288,7 @@ void musrRootOutput::BeginOfRunAction() {
}
if (store_odet_ID || store_odet_nPhot || store_odet_timeFirst || store_odet_timeA || store_odet_timeB ||
store_odet_timeC || store_odet_timeD || store_odet_timeE || store_odet_timeLast || store_odet_timeCFD || store_odet_amplCFD)
store_odet_timeC || store_odet_timeD || store_odet_timeMean || store_odet_timeLast || store_odet_timeCFD || store_odet_amplCFD)
{rootTree->Branch("odet_n",&odet_n,"odet_n/I");}
if (store_odet_ID) {rootTree->Branch("odet_ID",&odet_ID,"odet_ID[odet_n]/I");}
if (store_odet_nPhot) {rootTree->Branch("odet_nPhot",&odet_nPhot,"odet_nPhot[odet_n]/I");}
@ -292,7 +297,7 @@ void musrRootOutput::BeginOfRunAction() {
if (store_odet_timeB) {rootTree->Branch("odet_timeB",&odet_timeB,"odet_timeB[odet_n]/D");}
if (store_odet_timeC) {rootTree->Branch("odet_timeC",&odet_timeC,"odet_timeC[odet_n]/D");}
if (store_odet_timeD) {rootTree->Branch("odet_timeD",&odet_timeD,"odet_timeD[odet_n]/D");}
if (store_odet_timeE) {rootTree->Branch("odet_timeE",&odet_timeE,"odet_timeE[odet_n]/D");}
if (store_odet_timeMean) {rootTree->Branch("odet_timeMean",&odet_timeMean,"odet_timeMean[odet_n]/D");}
if (store_odet_timeLast) {rootTree->Branch("odet_timeLast",&odet_timeLast,"odet_timeLast[odet_n]/D");}
if (store_odet_timeCFD) {rootTree->Branch("odet_timeCFD",&odet_timeCFD,"odet_timeCFD[odet_n]/D");}
if (store_odet_amplCFD) {rootTree->Branch("odet_amplCFD",&odet_amplCFD,"odet_amplCFD[odet_n]/D");}
@ -534,7 +539,7 @@ void musrRootOutput::SetDetectorInfoVvv (G4int nDetectors,
void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA,
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast, G4double timeCFD, G4double amplCFD)
G4double timeB, G4double timeC, G4double timeD, G4double timeMean, G4double timeLast, G4double timeCFD, G4double amplCFD)
{
if ((nDetectors<0)||(nDetectors>=(odet_nMax-1))) {
char message[200];
@ -551,7 +556,7 @@ void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4
odet_timeB[nDetectors]=timeB/microsecond;
odet_timeC[nDetectors]=timeC/microsecond;
odet_timeD[nDetectors]=timeD/microsecond;
odet_timeE[nDetectors]=timeE/microsecond;
odet_timeMean[nDetectors]=timeMean/microsecond;
odet_timeLast[nDetectors]=timeLast/microsecond;
odet_timeCFD[nDetectors]=timeCFD/microsecond;
odet_amplCFD[nDetectors]=amplCFD;

31
src/musrScintSD.cc
View File

@ -72,8 +72,9 @@ musrScintSD::musrScintSD(G4String name)
OPSA_fracB = 0.05;
OPSA_fracC = 0.10;
OPSA_fracD = 0.20;
OPSA_fracE = 0.5;
bool_multimapOfEventIDsForOPSAhistosEXISTS = false;
bool_StoreThisOPSAhistSUMMED = false;
bool_StoreThisOPSAhistALL = false;
OPSAhistoNbin = 100;
OPSAhistoMin =0;
OPSAhistoMax = 10.;
@ -530,13 +531,14 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
G4bool boolStoreThisOPSAhist = false;
G4bool boolStoreThisOPSAhistSUMMED = false;
// G4bool boolStoreThisOPSAhistSUMMED = false;
G4int OPSA_detID= it->first;
optHitDetectorMapType* optHitDetectorMap = it->second;
// Check whether OPSA histograming of times of optical photon detection is required for this eventID.
if (bool_multimapOfEventIDsForOPSAhistosEXISTS) {
if (multimapOfEventIDsForOPSAhistos.find(eeeventID)!=multimapOfEventIDsForOPSAhistos.end()) {
if (bool_StoreThisOPSAhistALL) boolStoreThisOPSAhist = true;
else if (multimapOfEventIDsForOPSAhistos.find(eeeventID)!=multimapOfEventIDsForOPSAhistos.end()) {
// Now check whether the histogramming is required for the currently analysed detector.
std::pair<multimapOfEventIDsForOPSAhistos_Type::iterator,multimapOfEventIDsForOPSAhistos_Type::iterator> retOPSAhist;
multimapOfEventIDsForOPSAhistos_Type::iterator itOPSAhist;
@ -547,10 +549,10 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
if ( (itOPSAhist->second == 0) || (itOPSAhist->second == OPSA_detID) ) boolStoreThisOPSAhist = true;
}
}
if (multimapOfEventIDsForOPSAhistos.find(-1)!=multimapOfEventIDsForOPSAhistos.end()) {
// The user requires to store OPSA timing histograms summed up for all events together
boolStoreThisOPSAhistSUMMED = true;
}
// if (multimapOfEventIDsForOPSAhistos.find(-1)!=multimapOfEventIDsForOPSAhistos.end()) {
// // The user requires to store OPSA timing histograms summed up for all events together
// boolStoreThisOPSAhistSUMMED = true;
// }
}
if (optHitDetectorMap->empty()) continue;
@ -564,7 +566,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
G4double OPSA_timeB = -1000000;
G4double OPSA_timeC = -1000000;
G4double OPSA_timeD = -1000000;
G4double OPSA_timeE = -1000000;
G4double OPSA_timeMean = -1000000;
G4double OPSA_timeLast = -1000000;
G4double OPSA_CFD_time = -1000000;
G4double OPSA_CFD_ampl = -1000;
@ -593,7 +595,6 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
G4int NB = int (OPSA_fracB * OPSA_f_nPhot + 0.5); if (NB<=0) NB=1;
G4int NC = int (OPSA_fracC * OPSA_f_nPhot + 0.5); if (NC<=0) NC=1;
G4int ND = int (OPSA_fracD * OPSA_f_nPhot + 0.5); if (ND<=0) ND=1;
G4int NE = int (OPSA_fracE * OPSA_f_nPhot + 0.5); if (NE<=0) NE=1;
Int_t nP=0;
// Define OPSA histograms if required for this event
@ -614,7 +615,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
OPSA_CFD = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
}
}
if (boolStoreThisOPSAhistSUMMED) {
if (bool_StoreThisOPSAhistSUMMED) {
iHistNrSUM++;
char nameHist[200]; sprintf(nameHist,"OPSAhistSUM_%d_%d",OPSA_detID,iHistNrSUM);
char nameHist0[200]; sprintf(nameHist0,"OPSAhistSUM0_%d_%d",OPSA_detID,iHistNrSUM);
@ -642,10 +643,9 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
if (nP==NB) OPSA_timeB = timePhot;
if (nP==NC) OPSA_timeC = timePhot;
if (nP==ND) OPSA_timeD = timePhot;
if (nP==NE) OPSA_timeE = timePhot;
if (nP==OPSA_nPhot) OPSA_timeLast = timePhot;
if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) OPSAhisto->Fill(timePhot-OPSA_timeFirst+0.00000000001);
if (boolStoreThisOPSAhistSUMMED) {
if (bool_StoreThisOPSAhistSUMMED) {
OPSAhistoSUM ->Fill(timePhot-OPSA_timeFirst+0.00000000001);
OPSAhistoSUM0->Fill(timePhot+0.00000000001);
}
@ -653,6 +653,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
// OPSAhisto->Fit(poiss,"Q");
if (boolStoreThisOPSAhist) OPSAhisto->Write();
OPSA_timeMean = OPSAhisto->GetMean();
// if required, convert the histogram with photon times to histogram of electronic pulse shape
if (bool_pulseShapeExists) {
@ -674,7 +675,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
Double_t xValue = (iBin-0.5)*OPSAhistoBinWidth;
OPSA_CFD->Fill(xValue,signalValue*OPSA_CFD_a1);
Double_t xValueShifted = xValue + OPSA_CFD_delay;
if (xValueShifted>OPSAhistoMax) break; //return if out of range;
if (xValueShifted>OPSAhistoMax) continue; //return if out of range;
OPSA_CFD->Fill(xValueShifted,signalValue);
}
if (boolStoreThisOPSAhist) OPSA_CFD -> Write();
@ -708,7 +709,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
}
signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_timeFirst,OPSA_timeA,OPSA_timeB,OPSA_timeC,
OPSA_timeD,OPSA_timeE,OPSA_timeLast,OPSA_CFD_time,OPSA_CFD_ampl);
OPSA_timeD,OPSA_timeMean,OPSA_timeLast,OPSA_CFD_time,OPSA_CFD_ampl);
OPSA_signal_Map.insert(std::pair<G4int,signalInfo*>(OPSA_nPhot,mySignalInfo) );
}
OPSA_nPhot = 0;
@ -717,7 +718,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
OPSA_timeB = -1000000;
OPSA_timeC = -1000000;
OPSA_timeD = -1000000;
OPSA_timeE = -1000000;
OPSA_timeMean = -1000000;
OPSA_timeLast = -1000000;
OPSA_CFD_time = -1000000;
OPSA_CFD_ampl = -1000;

69
src/musrSteppingAction.cc
View File

@ -89,14 +89,21 @@ void musrSteppingAction::DoAtTheBeginningOfEvent() {
void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
G4Track* aTrack = aStep->GetTrack();
G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
G4ThreeVector preStepPosition = preStepPoint->GetPosition();
G4ThreeVector postStepPosition = postStepPoint->GetPosition();
// suspend the track if too many steps has already happened (relevant at high field)
if (aTrack->GetCurrentStepNumber()>100000) {
musrErrorMessage::GetInstance()->musrError(WARNING,
"musrSteppingAction: Current number of steps for the track > 100000 ==> TRACK KILLED",true);
G4double x=aStep->GetPostStepPoint()->GetPosition().x()/mm;
G4double y=aStep->GetPostStepPoint()->GetPosition().y()/mm;
G4double z=aStep->GetPostStepPoint()->GetPosition().z()/mm;
// G4double x=aStep->GetPostStepPoint()->GetPosition().x()/mm;
// G4double y=aStep->GetPostStepPoint()->GetPosition().y()/mm;
// G4double z=aStep->GetPostStepPoint()->GetPosition().z()/mm;
G4double x = postStepPosition.x()/mm;
G4double y = postStepPosition.y()/mm;
G4double z = postStepPosition.z()/mm;
G4double E=aTrack->GetVertexKineticEnergy()/MeV;
myRootOutput->htest1->Fill(x,y);
myRootOutput->htest2->Fill(sqrt(x*x+y*y),z);
@ -118,13 +125,14 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
// Temporary fix to avoid crashes caused by particles with unphysically high energies
// (probably corrupted event?)
if ((aStep->GetPreStepPoint()->GetKineticEnergy()) > (1*GeV)) {
// if ((aStep->GetPreStepPoint()->GetKineticEnergy()) > (1*GeV)) {
if ((preStepPoint->GetKineticEnergy()) > (1*GeV)) {
musrErrorMessage::GetInstance()->musrError(SERIOUS,
"musrSteppingAction: kinetic energy of a particle larger than 1GeV! STRANGE FOR muSR!",false);
G4RunManager* fRunManager = G4RunManager::GetRunManager();
G4cout<<" Event nr.:"<<fRunManager->GetCurrentEvent()->GetEventID()
<<", the particle \""<< aTrack->GetDynamicParticle()->GetDefinition()->GetParticleName()
<<"\" has energy of "<<(aStep->GetPreStepPoint()->GetKineticEnergy())/GeV<<" GeV."<<G4endl;
<<"\" has energy of "<<(preStepPoint->GetKineticEnergy())/GeV<<" GeV."<<G4endl;
G4cout<<" Deleting the event!"<<G4endl;
G4cout.flush();
myRootOutput->SetEventWeight(0);
@ -188,7 +196,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
if (aTrack->GetCreatorProcess()->GetProcessName() == "RadioactiveDecay") {
if (aTrack->GetDefinition()->GetParticleName()=="e-") {
if (aTrack->GetCurrentStepNumber()==1) {
G4double electron_kinetic_energy=aStep->GetPreStepPoint()->GetKineticEnergy();
G4double electron_kinetic_energy=preStepPoint->GetKineticEnergy();
myRootOutput->htest4->Fill(electron_kinetic_energy);
}
}
@ -211,17 +219,17 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
G4int tmpVolumeID=saveVolumeMapping[actualVolume];
if (tmpVolumeID!=0) {
G4int particle_id_save=p_definition->GetPDGEncoding();
G4double ke_save=aStep->GetPreStepPoint()->GetKineticEnergy();
G4double x_save=aStep->GetPreStepPoint()->GetPosition().x();
G4double y_save=aStep->GetPreStepPoint()->GetPosition().y();
G4double z_save=aStep->GetPreStepPoint()->GetPosition().z();
G4double time_save=aStep->GetPreStepPoint()->GetGlobalTime();
G4double px_save=aStep->GetPreStepPoint()->GetMomentum().x();
G4double py_save=aStep->GetPreStepPoint()->GetMomentum().y();
G4double pz_save=aStep->GetPreStepPoint()->GetMomentum().z();
G4double polx_save=aStep->GetPreStepPoint()->GetPolarization().x();
G4double poly_save=aStep->GetPreStepPoint()->GetPolarization().y();
G4double polz_save=aStep->GetPreStepPoint()->GetPolarization().z();
G4double ke_save=preStepPoint->GetKineticEnergy();
G4double x_save=preStepPosition.x();
G4double y_save=preStepPosition.y();
G4double z_save=preStepPosition.z();
G4double time_save=preStepPoint->GetGlobalTime();
G4double px_save=preStepPoint->GetMomentum().x();
G4double py_save=preStepPoint->GetMomentum().y();
G4double pz_save=preStepPoint->GetMomentum().z();
G4double polx_save=preStepPoint->GetPolarization().x();
G4double poly_save=preStepPoint->GetPolarization().y();
G4double polz_save=preStepPoint->GetPolarization().z();
myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,time_save,px_save,py_save,pz_save,polx_save,poly_save,polz_save);
}
}
@ -236,10 +244,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
muAlreadyWasInTargetInThisEvent=true;
myRootOutput->SetPolInTarget(aTrack->GetPolarization());
myRootOutput->SetTimeInTarget(aTrack->GetGlobalTime());
G4ThreeVector p_beforeEnteringTarget = aTrack->GetMomentum() - aStep->GetDeltaMomentum(); // must be - sign because DeltaMomentum is negative
// std::cout<<"(aStep->GetDeltaMomentum()).z() = "<<(aStep->GetDeltaMomentum()).z()<<std::endl;
myRootOutput->SetMomentumInTarget(p_beforeEnteringTarget);
// myRootOutput->SetMomentumInTarget(aTrack->GetMomentum());
myRootOutput->SetMomentumInTarget(preStepPoint->GetMomentum());
}
}
else if ((actualVolume=="log_M0")||(actualVolume=="log_m0")) {
@ -271,10 +276,9 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
// force Geant to use very small step size.
if (boolCalculateFieldIntegral) {
if (F04GlobalField::Exists()) {
G4ThreeVector position_tmp = aStep->GetPostStepPoint()->GetPosition();
CoordinateForFieldIntegral[0] = position_tmp.x();
CoordinateForFieldIntegral[1] = position_tmp.y();
CoordinateForFieldIntegral[2] = position_tmp.z();
CoordinateForFieldIntegral[0] = postStepPosition.x();
CoordinateForFieldIntegral[1] = postStepPosition.y();
CoordinateForFieldIntegral[2] = postStepPosition.z();
CoordinateForFieldIntegral[3] = aTrack->GetGlobalTime();
F04GlobalField::getObject() -> GetFieldValue(CoordinateForFieldIntegral,FieldForFieldIntegral);
// G4cout<<"B=("<<FieldForFieldIntegral[0]/tesla<<","<<FieldForFieldIntegral[1]/tesla<<","<<FieldForFieldIntegral[2]/tesla<<")"<<G4endl;
@ -285,15 +289,15 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
BxIntegral += stepLength * FieldForFieldIntegral[0];
ByIntegral += stepLength * FieldForFieldIntegral[1];
BzIntegral += stepLength * FieldForFieldIntegral[2];
BzIntegral1 += ( position_tmp.z() - aStep->GetPreStepPoint()->GetPosition().z() ) * FieldForFieldIntegral[2];
BzIntegral1 += ( postStepPosition.z() - preStepPoint->GetPosition().z() ) * FieldForFieldIntegral[2];
BzIntegral2 += stepLength;
BzIntegral3 += ( position_tmp.z() - aStep->GetPreStepPoint()->GetPosition().z() );
BzIntegral3 += ( postStepPosition.z() - preStepPoint->GetPosition().z() );
// G4cout<<"BzIntegral="<<BzIntegral<<" stepLength="<<stepLength<<"FieldForFieldIntegral[2]="<<FieldForFieldIntegral[2]<<G4endl;
}
}
// Pick up process "DecayWithSpin":
const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
const G4VProcess* process = postStepPoint->GetProcessDefinedStep();
if (process!=NULL) {
G4String processName = process->GetProcessName();
if (processName=="DecayWithSpin") {
@ -324,9 +328,8 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
// Store the information about the decaying muon and store it in the Root tree
G4double timeOfDecay_tmp=aTrack->GetGlobalTime();
G4ThreeVector positionOfDecay_tmp = aStep->GetPostStepPoint()->GetPosition();
G4double BFieldAtOrigin[6] = {0.,0.,0.,0.,0.,0.};
G4double PointOfDecay[4] ={positionOfDecay_tmp.x(),positionOfDecay_tmp.y(),positionOfDecay_tmp.z(),timeOfDecay_tmp};
G4double PointOfDecay[4] ={postStepPosition.x(),postStepPosition.y(),postStepPosition.z(),timeOfDecay_tmp};
if (F04GlobalField::Exists()) {
F04GlobalField* myGlobalField = F04GlobalField::getObject();
myGlobalField->GetFieldValue(PointOfDecay,BFieldAtOrigin);
@ -344,7 +347,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
}
myRootOutput->SetDecayTime(timeOfDecay_tmp);
myRootOutput->SetDecayPolarisation(aTrack->GetPolarization());
myRootOutput->SetDecayPosition(positionOfDecay_tmp);
myRootOutput->SetDecayPosition(postStepPosition);
myRootOutput->SetDecayDetectorID(actualVolume);
myRootOutput->SetBField(BFieldAtOrigin);
if (boolCalculateFieldIntegral) {
@ -362,10 +365,6 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
}
}
}
// G4ThreeVector position = aStep->GetPostStepPoint()->GetPosition();
// G4ThreeVector polarization=aTrack->GetDynamicParticle()->GetPolarization();
}
else { // particle is not muon