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

Kamil Sedlak 7.2.2013

Browse Source 
Implemented changes of James Lord.  Compiling was possible, but
no testing was done so far.
This commit is contained in:
sedlak 2013-02-07 14:54:04 +00:00
parent db7b6a80f4
commit c5fb2ca46f
20 changed files with 3967 additions and 2894 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
CMakeMacroParseArguments.cmake Normal file
View File

@ -0,0 +1,30 @@
MACRO(PARSE_ARGUMENTS prefix arg_names option_names)
SET(DEFAULT_ARGS)
FOREACH(arg_name ${arg_names})
SET(${prefix}_${arg_name})
ENDFOREACH(arg_name)
FOREACH(option ${option_names})
SET(${prefix}_${option} FALSE)
ENDFOREACH(option)
SET(current_arg_name DEFAULT_ARGS)
SET(current_arg_list)
FOREACH(arg ${ARGN})
SET(larg_names ${arg_names})
LIST(FIND larg_names "${arg}" is_arg_name)
IF (is_arg_name GREATER -1)
SET(${prefix}_${current_arg_name} ${current_arg_list})
SET(current_arg_name ${arg})
SET(current_arg_list)
ELSE (is_arg_name GREATER -1)
SET(loption_names ${option_names})
LIST(FIND loption_names "${arg}" is_option)
IF (is_option GREATER -1)
SET(${prefix}_${arg} TRUE)
ELSE (is_option GREATER -1)
SET(current_arg_list ${current_arg_list} ${arg})
ENDIF (is_option GREATER -1)
ENDIF (is_arg_name GREATER -1)
ENDFOREACH(arg)
SET(${prefix}_${current_arg_name} ${current_arg_list})
ENDMACRO(PARSE_ARGUMENTS)

138
CMakeParseArguments.cmake Normal file
View File

@ -0,0 +1,138 @@
# CMAKE_PARSE_ARGUMENTS(<prefix> <options> <one_value_keywords> <multi_value_keywords> args...)
#
# CMAKE_PARSE_ARGUMENTS() is intended to be used in macros or functions for
# parsing the arguments given to that macro or function.
# It processes the arguments and defines a set of variables which hold the
# values of the respective options.
#
# The <options> argument contains all options for the respective macro,
# i.e. keywords which can be used when calling the macro without any value
# following, like e.g. the OPTIONAL keyword of the install() command.
#
# The <one_value_keywords> argument contains all keywords for this macro
# which are followed by one value, like e.g. DESTINATION keyword of the
# install() command.
#
# The <multi_value_keywords> argument contains all keywords for this macro
# which can be followed by more than one value, like e.g. the TARGETS or
# FILES keywords of the install() command.
#
# When done, CMAKE_PARSE_ARGUMENTS() will have defined for each of the
# keywords listed in <options>, <one_value_keywords> and
# <multi_value_keywords> a variable composed of the given <prefix>
# followed by "_" and the name of the respective keyword.
# These variables will then hold the respective value from the argument list.
# For the <options> keywords this will be TRUE or FALSE.
#
# All remaining arguments are collected in a variable
# <prefix>_UNPARSED_ARGUMENTS, this can be checked afterwards to see whether
# your macro was called with unrecognized parameters.
#
# As an example here a my_install() macro, which takes similar arguments as the
# real install() command:
#
# function(MY_INSTALL)
# set(options OPTIONAL FAST)
# set(oneValueArgs DESTINATION RENAME)
# set(multiValueArgs TARGETS CONFIGURATIONS)
# cmake_parse_arguments(MY_INSTALL "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} )
# ...
#
# Assume my_install() has been called like this:
# my_install(TARGETS foo bar DESTINATION bin OPTIONAL blub)
#
# After the cmake_parse_arguments() call the macro will have set the following
# variables:
# MY_INSTALL_OPTIONAL = TRUE
# MY_INSTALL_FAST = FALSE (this option was not used when calling my_install()
# MY_INSTALL_DESTINATION = "bin"
# MY_INSTALL_RENAME = "" (was not used)
# MY_INSTALL_TARGETS = "foo;bar"
# MY_INSTALL_CONFIGURATIONS = "" (was not used)
# MY_INSTALL_UNPARSED_ARGUMENTS = "blub" (no value expected after "OPTIONAL"
#
# You can the continue and process these variables.
#
# Keywords terminate lists of values, e.g. if directly after a one_value_keyword
# another recognized keyword follows, this is interpreted as the beginning of
# the new option.
# E.g. my_install(TARGETS foo DESTINATION OPTIONAL) would result in
# MY_INSTALL_DESTINATION set to "OPTIONAL", but MY_INSTALL_DESTINATION would
# be empty and MY_INSTALL_OPTIONAL would be set to TRUE therefor.
#=============================================================================
# Copyright 2010 Alexander Neundorf <neundorf@kde.org>
#
# Distributed under the OSI-approved BSD License (the "License");
# see accompanying file Copyright.txt for details.
#
# This software is distributed WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the License for more information.
#=============================================================================
# (To distribute this file outside of CMake, substitute the full
# License text for the above reference.)
if(__CMAKE_PARSE_ARGUMENTS_INCLUDED)
return()
endif()
set(__CMAKE_PARSE_ARGUMENTS_INCLUDED TRUE)
function(CMAKE_PARSE_ARGUMENTS prefix _optionNames _singleArgNames _multiArgNames)
# first set all result variables to empty/FALSE
foreach(arg_name ${_singleArgNames} ${_multiArgNames})
set(${prefix}_${arg_name})
endforeach(arg_name)
foreach(option ${_optionNames})
set(${prefix}_${option} FALSE)
endforeach(option)
set(${prefix}_UNPARSED_ARGUMENTS)
set(insideValues FALSE)
set(currentArgName)
# now iterate over all arguments and fill the result variables
foreach(currentArg ${ARGN})
list(FIND _optionNames "${currentArg}" optionIndex) # ... then this marks the end of the arguments belonging to this keyword
list(FIND _singleArgNames "${currentArg}" singleArgIndex) # ... then this marks the end of the arguments belonging to this keyword
list(FIND _multiArgNames "${currentArg}" multiArgIndex) # ... then this marks the end of the arguments belonging to this keyword
if(${optionIndex} EQUAL -1 AND ${singleArgIndex} EQUAL -1 AND ${multiArgIndex} EQUAL -1)
if(insideValues)
if("${insideValues}" STREQUAL "SINGLE")
set(${prefix}_${currentArgName} ${currentArg})
set(insideValues FALSE)
elseif("${insideValues}" STREQUAL "MULTI")
list(APPEND ${prefix}_${currentArgName} ${currentArg})
endif()
else(insideValues)
list(APPEND ${prefix}_UNPARSED_ARGUMENTS ${currentArg})
endif(insideValues)
else()
if(NOT ${optionIndex} EQUAL -1)
set(${prefix}_${currentArg} TRUE)
set(insideValues FALSE)
elseif(NOT ${singleArgIndex} EQUAL -1)
set(currentArgName ${currentArg})
set(${prefix}_${currentArgName})
set(insideValues "SINGLE")
elseif(NOT ${multiArgIndex} EQUAL -1)
set(currentArgName ${currentArg})
set(${prefix}_${currentArgName})
set(insideValues "MULTI")
endif()
endif()
endforeach(currentArg)
# propagate the result variables to the caller:
foreach(arg_name ${_singleArgNames} ${_multiArgNames} ${_optionNames})
set(${prefix}_${arg_name} ${${prefix}_${arg_name}} PARENT_SCOPE)
endforeach(arg_name)
set(${prefix}_UNPARSED_ARGUMENTS ${${prefix}_UNPARSED_ARGUMENTS} PARENT_SCOPE)
endfunction(CMAKE_PARSE_ARGUMENTS _options _singleArgs _multiArgs)

BIN
doc/musrSim.pdf
View File

Binary file not shown.

5
doc/musrSim.tex
View File

@ -23,10 +23,11 @@
%\title{GEANT Simulation Program for the %\title{GEANT Simulation Program for the
%$\mathbf{\mu}$SR Instruments} %$\mathbf{\mu}$SR Instruments}
\author{Kamil Sedl\'ak$^1$, Toni Shiroka$^1$, Zaher Salman$^1$, Jose Rodriguez$^1$, Tom Lancaster$^2$, Thomas Prokscha$^1$, Taofiq Para\"{\i}so$^1$, Robert Scheuermann$^1$} \author{Kamil Sedl\'ak$^1$, Toni Shiroka$^1$, Zaher Salman$^1$, Jose Rodriguez$^1$, Tom Lancaster$^2$, Thomas Prokscha$^1$, Taofiq Para\"{\i}so$^1$, Robert Scheuermann$^1$, James S. Lord$^3$}
\address{{$^1$ Laboratory for Muon Spin Spectroscopy, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland}\\ \address{{$^1$ Laboratory for Muon Spin Spectroscopy, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland}\\
$^2$ Clarendon Laboratory, Department of Physics, Oxford University, Parks Road, Oxford OX1 3PU, UK} $^2$ Clarendon Laboratory, Department of Physics, Oxford University, Parks Road, Oxford OX1 3PU, UK\\
$^3$ ISIS Facility, Rutherford Appleton Laboratory, Chilton, Oxon OX11 0QX, U.K.}
%\address{{\rm (on behalf of the H1 collaboration)}\\ %\address{{\rm (on behalf of the H1 collaboration)}\\
%Institute of Physics AS\,CR\\ %Institute of Physics AS\,CR\\
%%Academy of Sciences of the Czech Republic %%Academy of Sciences of the Czech Republic

BIN
doc/musrSimAna.pdf
View File

Binary file not shown.

64
doc/musrSimAna.tex
View File

@ -9,8 +9,10 @@
\newcommand{\musrSimAna}{\emph{musrSimAna}} \newcommand{\musrSimAna}{\emph{musrSimAna}}
\title{Manual of musrSimAna} \title{Manual of musrSimAna}
\author{Kamil Sedl\'ak$^1$} \author{Kamil Sedl\'ak$^1$}
\author{James S. Lord$^2$}
\address{{$^1$ Laboratory for Muon Spin Spectroscopy, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland}} \address{{$^1$ Laboratory for Muon Spin Spectroscopy, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland}}
\address{{$^2$ ISIS Facility, Rutherford Appleton Laboratory, Chilton, Oxon OX11 0QX, U.K.}}
\maketitle \maketitle
@ -41,8 +43,8 @@ hand, it requires special program for each detector set-up.
A second possibility is to use \musrSimAna, which is intended to be a general analysis A second possibility is to use \musrSimAna, which is intended to be a general analysis
program for \musr\ instruments. program for \musr\ instruments.
At the moment, however, \musrSimAna\ is tuned to the needs of continuous muon beam facilities. \musrSimAna\ was originally written for continuous muon beam facilities.
Some modifications might be necessary for the pulsed muon beams. Some extensions have been added for pulsed muon beams, see section \ref{sect_pulsed}.
As in the case of \musrSim, the user of \musrSimAna\ specifies all parameters As in the case of \musrSim, the user of \musrSimAna\ specifies all parameters
required by the analysis in an input text file (steering file). The initial structure required by the analysis in an input text file (steering file). The initial structure
@ -276,12 +278,24 @@ All events should/have to be (?) saved in the Root tree
two subsequently generated muons, and therefore decreases the incoming muon rate two subsequently generated muons, and therefore decreases the incoming muon rate
(number of muons per second). This parameter should be changed in order to set the (number of muons per second). This parameter should be changed in order to set the
incoming muon rate to a given required value, typically to 30\,000\,$\mu/$s.\\ incoming muon rate to a given required value, typically to 30\,000\,$\mu/$s.\\
For Pulsed beams this adjusts the mean number of muons per frame and can again be used to match the experimental count rate.\\
See also variable ``INFINITELYLOWMUONRATE''. See also variable ``INFINITELYLOWMUONRATE''.
\item{\bf INFINITELYLOWMUONRATE} \\ \item{\bf INFINITELYLOWMUONRATE} \\
If INFINITELYLOWMUONRATE is specified, each event is treated independently of any other If INFINITELYLOWMUONRATE is specified, each event is treated independently of any other
event. This corresponds to a situation of infinitely low rate of incoming muons, and event. This corresponds to a situation of infinitely low rate of incoming muons, and
no pileup can be observed. The variable ``MUONRATEFACTOR'' becomes irrelevant when no pileup can be observed. The variable ``MUONRATEFACTOR'' becomes irrelevant when
INFINITELYLOWMUONRATE is specified. INFINITELYLOWMUONRATE is specified.
\item{\bf FRAMEINTERVAL=\emph{value}} \\
For pulsed data analysis, the interval between muon pulses in ms. Affects the number of muons/frame (as does MUONRATEFACTOR) and the calculated run duration and count rate.
\item{\bf PARTIALFRAMEATEND} \\
For pulsed data analysis. If the events in the input file do not fill an exact number of frames, this requests that the partially filled frame left over at the end will be processed too and included in the results.
This will happen anyway if it would be the only frame (low statistics input file or very high requested flux). Normally such a frame would be discarded as it would have the wrong pile-up fraction.
\item{\bf DEADTIME=\emph{value}} \\
Set the ``dead time'' recovery time in ns for pulsed data analysis, for all the detectors defined after this point in the file.
\item{\bf MUONPULSEWIDTHFACTOR=\emph{value}} \\
Scale the random start times (and therefore the muon pulse width) in the input file by this value.
\item{\bf COMMONTHRESHOLD=\emph{value}} \\
A short-cut to adjust the threshold settings for all detectors (after this point in the file). Overrides the individual values.
\item{\bf DATAWINDOWMIN=\emph{value}} \\ \item{\bf DATAWINDOWMIN=\emph{value}} \\
Beginning of the data interval for the positron counters in $\mu$s. Beginning of the data interval for the positron counters in $\mu$s.
\item{\bf DATAWINDOWMAX=\emph{value}} \\ \item{\bf DATAWINDOWMAX=\emph{value}} \\
@ -300,7 +314,7 @@ All events should/have to be (?) saved in the Root tree
for {\tt PROMPTPEAKMIN}.) for {\tt PROMPTPEAKMIN}.)
\item{\bf MUSRMODE=\emph{string}} \\ \item{\bf MUSRMODE=\emph{string}} \\
Defines the mode of \musr\ experiment -- presently only ``D'', corresponding to Defines the mode of \musr\ experiment -- presently only ``D'', corresponding to
the time differential mode is implemented. the time differential mode, and ``P'', for pulsed mode, are implemented.
\item{\bf REWINDTIMEBINS=\emph{value}} \\ \item{\bf REWINDTIMEBINS=\emph{value}} \\
A technical parameter specifying when a roll-over of all hits has to be done. A technical parameter specifying when a roll-over of all hits has to be done.
It is specified in TDC bins, and normally there should be no need to change this parameter. It is specified in TDC bins, and normally there should be no need to change this parameter.
@ -381,6 +395,14 @@ All events should/have to be (?) saved in the Root tree
\item[pileup\_muDecayDetID] \ldots detector ID, in which $\mu_2$ stopped and decayed. \item[pileup\_muDecayDetID] \ldots detector ID, in which $\mu_2$ stopped and decayed.
\item[pileup\_muDecayPosZ] \ldots $z$-coordinate of where the $\mu_2$ stopped and decayed. \item[pileup\_muDecayPosZ] \ldots $z$-coordinate of where the $\mu_2$ stopped and decayed.
\item[pileup\_muDecayPosR] \ldots radius of where the $\mu_2$ stopped and decayed. \item[pileup\_muDecayPosR] \ldots radius of where the $\mu_2$ stopped and decayed.
\item[BFieldAtDecay\_Bx] \ldots The field at the location the muon decayed (X component).
\item[BFieldAtDecay\_By] \ldots The field at the location the muon decayed (Y).
\item[BFieldAtDecay\_Bz] \ldots The field at the location the muon decayed (Z).
\item[fieldNomVal0] \ldots The nominal magnitude of the first field map used. Useful if the field is scanned in \musrSim \ and some background muons stop away from the central uniform field region.
\item[fieldNomVal1] \ldots The nominal magnitude of the second mapped field.
\item[pos\_detID] \ldots The detector number of the positron detector recording this signal.
\item[pos\_detID\_doublehit] \ldots For pulsed mode, if this is not the first positron signal, the detector number of that first signal.
\item[multiHitInterval] \ldots (Pulsed) If a double count, the time interval (in $\mu$s) between the first positron signal and this one.
\end{description} \end{description}
Variables are usually set to -1000 if they can not be calculated (e.g.\ {\tt det\_posEdep} = -1000 Variables are usually set to -1000 if they can not be calculated (e.g.\ {\tt det\_posEdep} = -1000
if there was no hit in any positron counter). if there was no hit in any positron counter).
@ -401,22 +423,23 @@ All events should/have to be (?) saved in the Root tree
The \emph{conditionName} is one of the following: The \emph{conditionName} is one of the following:
\begin{description} \begin{description}
\item[alwaysTrue] \ldots true for every hit in the m-counter (there can be more than one M-counter hit per event). \item[alwaysTrue] \ldots true for every hit in the m-counter (there can be more than one M-counter hit per event).
\item[oncePerEvent] \ldots true once for every event (the first hit in M-counter, if any, is considered). \item[oncePerEvent] \ldots true once for every event (the first hit in M-counter, if any, is considered). For Pulsed mode, this is the first hit in the lowest numbered positron detector, if any.
\item[muonDecayedInSample\_gen] \ldots true if muon stopped and decayed in the sample. \item[muonDecayedInSample\_gen] \ldots true if muon stopped and decayed in the sample (one count for each positron hit, if double counting in pulsed mode).
\item[muonDecayedInSampleOnce\_gen] \ldots true if muon stopped and decayed in the sample (true only once mer muon, use for sample/background fractions).
\item[muonTriggered\_gen] \ldots true if muon passed through the M-counter (irrespective of the deposited energy) \item[muonTriggered\_gen] \ldots true if muon passed through the M-counter (irrespective of the deposited energy)
-- not a measurable variable. -- not a measurable variable.
\item[muonTriggered\_det] \ldots true if a good muon candidate was found in the M-counter (using coincidences, vetoes, ...). \item[muonTriggered\_det] \ldots true if a good muon candidate was found in the M-counter (using coincidences, vetoes, ...).
Double hits within the pileup window are excluded. Double hits within the pileup window are excluded.
\item[positronHit\_det] \ldots true if a good positron candidate was found in the positron counter. \item[positronHit\_det] \ldots true if a good positron candidate was found in the positron counter.
Double hits within the data window are excluded. Double hits within the data window are excluded.
\item[goodEvent\_det] \ldots true if {\tt muonTriggered\_det} and {\tt positronHit\_det}. \item[goodEvent\_det] \ldots true if {\tt muonTriggered\_det} and {\tt positronHit\_det}. In Pulsed mode, any positron event over threshold and not killed by pileup.
\item[goodEvent\_gen] \ldots true if muon passed through the M-counter, and the muon stopped anywhere \item[goodEvent\_gen] \ldots true if muon passed through the M-counter, and the muon stopped anywhere
(i.e.\ did not leave the World volume of the simulation). No requirement (i.e.\ did not leave the World volume of the simulation). No requirement
on the positron is implied, i.e.\ the positron may or may not be detected. on the positron is implied, i.e.\ the positron may or may not be detected.
Not a measurable variable. Not a measurable variable. For pulsed mode, the first (or only) positron event over threshold (even if piled up).
\item[goodEvent\_det\_AND\_goodEvent\_gen] \item[goodEvent\_det\_AND\_goodEvent\_gen]
\item[pileupEventCandidate] \ldots M-counter hit and positron counter hit both come from two different events. \item[pileupEventCandidate] \ldots M-counter hit and positron counter hit both come from two different events.
\item[pileupEvent] \ldots {\tt pileupEventCandidate} and {\tt goodEvent\_det}. \item[pileupEvent] \ldots {\tt pileupEventCandidate} and {\tt goodEvent\_det}. Or for pulsed mode, an event lost due to ``dead time''.
\item[goodEvent\_det\_AND\_muonDecayedInSample\_gen] \item[goodEvent\_det\_AND\_muonDecayedInSample\_gen]
\item[goodEvent\_F\_det] \ldots {\tt goodEvent\_det}, where the positron was detected in the forward detectors \item[goodEvent\_F\_det] \ldots {\tt goodEvent\_det}, where the positron was detected in the forward detectors
defined by the command {\tt counterGrouping}. defined by the command {\tt counterGrouping}.
@ -425,7 +448,8 @@ All events should/have to be (?) saved in the Root tree
\item[goodEvent\_D\_det] \ldots like {\tt goodEvent\_F\_det} but for lower positron counters. \item[goodEvent\_D\_det] \ldots like {\tt goodEvent\_F\_det} but for lower positron counters.
\item[goodEvent\_L\_det] \ldots like {\tt goodEvent\_F\_det} but for left positron counters. \item[goodEvent\_L\_det] \ldots like {\tt goodEvent\_F\_det} but for left positron counters.
\item[goodEvent\_R\_det] \ldots like {\tt goodEvent\_F\_det} but for right positron counters. \item[goodEvent\_R\_det] \ldots like {\tt goodEvent\_F\_det} but for right positron counters.
\item[goodEvent\_F\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_F\_det} and {\tt pileupEvent}. \item[goodEvent\_F\_det\_AND\_muonDecayedInSample\_gen] \ldots the count in a Forward counter came from a muon in the sample. Equivalent conditions for the other groups.
\item[goodEvent\_F\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_F\_det} and {\tt pileupEvent}. For pulsed mode, an event lost from the F counters.
\item[goodEvent\_B\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_B\_det} and {\tt pileupEvent}. \item[goodEvent\_B\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_B\_det} and {\tt pileupEvent}.
\item[goodEvent\_U\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_U\_det} and {\tt pileupEvent}. \item[goodEvent\_U\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_U\_det} and {\tt pileupEvent}.
\item[goodEvent\_D\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_D\_det} and {\tt pileupEvent}. \item[goodEvent\_D\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_D\_det} and {\tt pileupEvent}.
@ -438,6 +462,11 @@ All events should/have to be (?) saved in the Root tree
\item[promptPeakD] \ldots like {\tt goodEvent\_D\_det} and {\tt promptPeak}. \item[promptPeakD] \ldots like {\tt goodEvent\_D\_det} and {\tt promptPeak}.
\item[promptPeakL] \ldots like {\tt goodEvent\_L\_det} and {\tt promptPeak}. \item[promptPeakL] \ldots like {\tt goodEvent\_L\_det} and {\tt promptPeak}.
\item[promptPeakR] \ldots like {\tt goodEvent\_R\_det} and {\tt promptPeak}. \item[promptPeakR] \ldots like {\tt goodEvent\_R\_det} and {\tt promptPeak}.
\item[doubleHitEvent\_gen] \ldots The second (or subsequent) signal in a multiple hit event event regardless of pileup
\item[doubleHit] \ldots A second (or subsequent) signal that is actually counted
\item[goodEvent\_F\_det\_AND\_doubleHit] \ldots A double event counted in the Forward counters. For pulsed analysis. Note that the first hit might not have been in the same bank. (Equivalent conditions for the other banks)
\item[singleHitEvent\_gen] \ldots A single hit, only one positron signal from this muon was over threshold. (Pulsed analysis)
\item[stackedEvent\_gen] \ldots The opposite of pileup, two or more small events (or noise) within a small time interval adding up to the threshold level. (Pulsed)
\end{description} \end{description}
Additional conditions may be implemented on request. Additional conditions may be implemented on request.
\item{\bf draw \emph{histogramName} \emph{conditionID} } \\ \item{\bf draw \emph{histogramName} \emph{conditionID} } \\
@ -794,6 +823,23 @@ Many different ``*.v1190'' files are stored in the file:
Note that the syntax of the ``fit'' command was changed Note that the syntax of the ``fit'' command was changed
at some point, and therefore the ``fit'' command might cause problems at some point, and therefore the ``fit'' command might cause problems
(the {\tt ``option''} has to be added in the old ``*.v1190'' files). (the {\tt ``option''} has to be added in the old ``*.v1190'' files).
%========================================================================================================
\section{Pulsed data collection}
\label{sect_pulsed}
The differences for a pulsed muon beam and instrument are:
\begin{itemize}
\item No ``M'' counter needs to be defined (but if present in \musrSim`s .mac file it must also be defined in the .v1190 file and will be ignored)
\item The \musrSim \ command ``/gun/starttimesigma \emph{sigma}'' defines the pulse width of the generated muon beam (in ns). A positive value is a Gaussian distribution, negative gives a rectangular pulse extending to $\pm$sigma.
% \item When varying the magnetic field or any physical arrangement of the instrument, \musrSim \ command ``/musr/command storeOnlyEventsWithHits false'' must be given so that the number of muons per frame, indicated count rate, and pileup, scale correctly with field.
\item It is advisable to set RESOLUTION to a very small value (1\,ps) so that pileup events do not coincide exactly.
\item \musrSimAna \ outputs the indicated count rate (in Mevents per hour) and this should be adjusted to match the experimental rate using the parameter MUONRATEFACTOR.
\item For ISIS operation with 4 pulses to TS1 (and muons) and one pulse to TS2, set FRAMEINTERVAL to the mean value of 25ms = 1/(40Hz).
\end{itemize}
\subsection{Dead time and pileup}
The model used is that each event ``charges up'' the discriminator by an amount equal to the energy deposited, and it then decays exponentially back towards zero with a time constant set by DEADTIME=\emph{value}. If the event causes the discriminator to cross the specified threshold in the positive direction, it will be recorded. At present there is no hysteresis and no minimum time interval between hits recorded in the TDC.
There are two likely artefacts. The first is the expected pileup or dead time distortion: a second event (itself over threshold) comes before the discriminator has recovered from the first, and will be missed. The second, referred to here as ``stacked'' events, are when a first small event charges the discriminator to half way and a second one soon after it can then reach the threshold and trigger a count. This is more common if the threshold is set too high and there are many normal positron hits not reaching the threshold on their own, and results in an opposite distortion to dead time, though with the same time dependence (worse for early times).
%======================================================================================================== %========================================================================================================
\begin{thebibliography}{0} \begin{thebibliography}{0}

6
musrSim.cc
View File

@ -10,7 +10,9 @@
#include "G4RunManager.hh" #include "G4RunManager.hh"
#include "G4UImanager.hh" #include "G4UImanager.hh"
#include "G4UIterminal.hh" #include "G4UIterminal.hh"
#include "G4UItcsh.hh" #ifdef G4UI_USE_TCSH
#include "G4UItcsh.hh" //JSL: should be commented on windows ?
#endif
//#include <TApplication.h> //#include <TApplication.h>
//#include <TSystem.h> //#include <TSystem.h>
@ -24,7 +26,7 @@
#include "Randomize.hh" #include "Randomize.hh"
#include <X11/Xlib.h> // #include <X11/Xlib.h> //JSL
#ifdef G4VIS_USE #ifdef G4VIS_USE
// #include "musrVisManager.hh" // #include "musrVisManager.hh"

612
musrSimAna/musrAnalysis.cxx
View File

@ -53,7 +53,7 @@ void musrAnalysis::Loop(char* runChar, char* v1190FileName, Int_t nrEvents)
if (ientry < 0) break; if (ientry < 0) break;
nb = fChain->GetEntry(jentry); InitialiseEvent(); nbytes += nb; nb = fChain->GetEntry(jentry); InitialiseEvent(); nbytes += nb;
// if (Cut(ientry) < 0) continue; // if (Cut(ientry) < 0) continue;
if (jentry%100000==0) std::cout << "Analysing event nr. jentry=" << jentry << std::endl; if (jentry%10000==0) std::cout << "Analysing event nr. jentry=" << jentry << " and found " << numberOfGoodMuons << " muons so far" << std::endl;
AnalyseEvent(jentry); AnalyseEvent(jentry);
} }
SaveHistograms(runChar,v1190FileName); SaveHistograms(runChar,v1190FileName);
@ -103,6 +103,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
pcoincwin = 40; pcoincwin = 40;
vcoincwin = 80; vcoincwin = 80;
muonRateFactor = 1.; muonRateFactor = 1.;
muonPulseWidthFactor = 1.;
// rewindTimeBins = 1000000000; // rewindTimeBins = 1000000000;
dataWindowMin = -0.2; dataWindowMin = -0.2;
dataWindowMax = 10.0; dataWindowMax = 10.0;
@ -113,6 +114,11 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
musrMode = 'D'; musrMode = 'D';
overallBinDelay = 0; overallBinDelay = 0;
boolInfinitelyLowMuonRate = false; boolInfinitelyLowMuonRate = false;
frameInterval = 25000.0 * microsecond; // JSL: ISIS "40Hz" average with TS1 & TS2 running
detRecoveryTime=10.0 * nanosecond; // JSL: typical ISIS PMT detector
doPartialFrameAtEnd = false; // tidy end of run
// processInBulk = true;
commonThreshold = -1000. ; // allow individual settings
// globTime = 0.; // globTime = 0.;
int tmpivar; int tmpivar;
@ -174,9 +180,31 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
sscanf(&line[strlen("MUONRATEFACTOR=")],"%g",&tmpfvar); sscanf(&line[strlen("MUONRATEFACTOR=")],"%g",&tmpfvar);
muonRateFactor = tmpfvar; muonRateFactor = tmpfvar;
} }
else if (strncmp(tmpString0,"MUONPULSEWIDTHFACTOR=",strlen("MUONPULSEWIDTHFACTOR="))==0) {
sscanf(&line[strlen("MUONPULSEWIDTHFACTOR=")],"%g",&tmpfvar);
muonPulseWidthFactor = tmpfvar;
}
else if (strncmp(tmpString0,"INFINITELYLOWMUONRATE=",strlen("INFINITELYLOWMUONRATE"))==0) { else if (strncmp(tmpString0,"INFINITELYLOWMUONRATE=",strlen("INFINITELYLOWMUONRATE"))==0) {
boolInfinitelyLowMuonRate=true; boolInfinitelyLowMuonRate=true;
} }
else if (strncmp(tmpString0,"FRAMEINTERVAL=",strlen("FRAMEINTERVAL="))==0) { // JSL: pulsed mode params. Frame interval in ms
sscanf(&line[strlen("FRAMEINTERVAL=")],"%g",&tmpfvar);
frameInterval = tmpfvar * microsecond * 1000.0;
}
else if (strncmp(tmpString0,"DEADTIME=",strlen("DEADTIME="))==0) { // JSL: pulsed mode params, dead time in ns. Redefine part way between detectors if not all the same
sscanf(&line[strlen("DEADTIME=")],"%g",&tmpfvar);
detRecoveryTime = tmpfvar * nanosecond;
}
else if (strncmp(tmpString0,"COMMONTHRESHOLD=",strlen("COMMONTHRESHOLD="))==0) { // JSL: pulsed mode params, dead time in ns. Redefine part way between detectors if not all the same
sscanf(&line[strlen("COMMONTHRESHOLD=")],"%g",&tmpfvar);
commonThreshold = tmpfvar;
}
else if (strncmp(tmpString0,"PARTIALFRAMEATEND=",strlen("PARTIALFRAMEATEND"))==0) {
doPartialFrameAtEnd=true;
}
// else if (strncmp(tmpString0,"PROCESSINBULK=",strlen("PROCESSINBULK"))==0) {
//processInBulk=true;
// }
else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) { else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) {
sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar); sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar);
if (tmpivar<0) {rewindTimeBins = Long64_t(tmpivar)* Long64_t(tmpivar);} if (tmpivar<0) {rewindTimeBins = Long64_t(tmpivar)* Long64_t(tmpivar);}
@ -411,6 +439,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
if (strcmp(conditionNameTMP,"alwaysTrue")==0) conditionMap[iConditionTMP]=&alwaysTrue; if (strcmp(conditionNameTMP,"alwaysTrue")==0) conditionMap[iConditionTMP]=&alwaysTrue;
else if (strcmp(conditionNameTMP,"oncePerEvent")==0) conditionMap[iConditionTMP]=&oncePerEvent; else if (strcmp(conditionNameTMP,"oncePerEvent")==0) conditionMap[iConditionTMP]=&oncePerEvent;
else if (strcmp(conditionNameTMP,"muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&muonDecayedInSample_gen; else if (strcmp(conditionNameTMP,"muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"muonDecayedInSampleOnce_gen")==0) conditionMap[iConditionTMP]=&muonDecayedInSampleOnce_gen;
else if (strcmp(conditionNameTMP,"muonTriggered_gen")==0) conditionMap[iConditionTMP]=&muonTriggered_gen; else if (strcmp(conditionNameTMP,"muonTriggered_gen")==0) conditionMap[iConditionTMP]=&muonTriggered_gen;
else if (strcmp(conditionNameTMP,"muonTriggered_gen_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&muonTriggered_gen_AND_muonDecayedInSample_gen; else if (strcmp(conditionNameTMP,"muonTriggered_gen_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&muonTriggered_gen_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"muonTriggered_det")==0) conditionMap[iConditionTMP]=&muonTriggered_det; else if (strcmp(conditionNameTMP,"muonTriggered_det")==0) conditionMap[iConditionTMP]=&muonTriggered_det;
@ -427,6 +456,12 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
else if (strcmp(conditionNameTMP,"goodEvent_D_det")==0) conditionMap[iConditionTMP]=&goodEvent_D_det; else if (strcmp(conditionNameTMP,"goodEvent_D_det")==0) conditionMap[iConditionTMP]=&goodEvent_D_det;
else if (strcmp(conditionNameTMP,"goodEvent_L_det")==0) conditionMap[iConditionTMP]=&goodEvent_L_det; else if (strcmp(conditionNameTMP,"goodEvent_L_det")==0) conditionMap[iConditionTMP]=&goodEvent_L_det;
else if (strcmp(conditionNameTMP,"goodEvent_R_det")==0) conditionMap[iConditionTMP]=&goodEvent_R_det; else if (strcmp(conditionNameTMP,"goodEvent_R_det")==0) conditionMap[iConditionTMP]=&goodEvent_R_det;
else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_D_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_D_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_L_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_L_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_R_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_R_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_pileupEvent; else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_pileupEvent; else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_pileupEvent; else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_pileupEvent;
@ -441,6 +476,15 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
else if (strcmp(conditionNameTMP,"promptPeakL")==0) conditionMap[iConditionTMP]=&promptPeakL; else if (strcmp(conditionNameTMP,"promptPeakL")==0) conditionMap[iConditionTMP]=&promptPeakL;
else if (strcmp(conditionNameTMP,"promptPeakR")==0) conditionMap[iConditionTMP]=&promptPeakR; else if (strcmp(conditionNameTMP,"promptPeakR")==0) conditionMap[iConditionTMP]=&promptPeakR;
else if (strcmp(conditionNameTMP,"doubleHit")==0) conditionMap[iConditionTMP]=&doubleHit; else if (strcmp(conditionNameTMP,"doubleHit")==0) conditionMap[iConditionTMP]=&doubleHit;
else if (strcmp(conditionNameTMP,"doubleHitEvent_gen")==0) conditionMap[iConditionTMP]=&doubleHitEvent_gen;
else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_doubleHit")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_doubleHit;
else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_doubleHit")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_doubleHit;
else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_doubleHit")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_doubleHit;
else if (strcmp(conditionNameTMP,"goodEvent_D_det_AND_doubleHit")==0) conditionMap[iConditionTMP]=&goodEvent_D_det_AND_doubleHit;
else if (strcmp(conditionNameTMP,"goodEvent_L_det_AND_doubleHit")==0) conditionMap[iConditionTMP]=&goodEvent_L_det_AND_doubleHit;
else if (strcmp(conditionNameTMP,"goodEvent_R_det_AND_doubleHit")==0) conditionMap[iConditionTMP]=&goodEvent_R_det_AND_doubleHit;
else if (strcmp(conditionNameTMP,"singleHitEvent_gen")==0) conditionMap[iConditionTMP]=&singleHitEvent_gen;
else if (strcmp(conditionNameTMP,"stackedEvent_gen")==0) conditionMap[iConditionTMP]=&stackedEvent_gen;
else { else {
std::cout<<" !!! ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx S T O P !!!"<<std::endl; std::cout<<" !!! ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx S T O P !!!"<<std::endl;
exit(1); exit(1);
@ -710,8 +754,10 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
int tmpChannel=-1, tmpTimeShift=0; int tmpChannel=-1, tmpTimeShift=0;
char tmpName[200], tmpType[200]; char tmpName[200], tmpType[200];
float tmpThreshold =0.; float tmpThreshold =0.;
Bool_t tmp_keep=(musrMode=='P'); // JSL: keep the small pulses in case different events add up to a big one
sscanf(&line[0],"%d %s %s %g %d",&tmpChannel,tmpName,tmpType,&tmpThreshold,&tmpTimeShift); sscanf(&line[0],"%d %s %s %g %d",&tmpChannel,tmpName,tmpType,&tmpThreshold,&tmpTimeShift);
musrCounter* counter = new musrCounter(tmpChannel,tmpName,tmpType[0],tmpThreshold,tmpTimeShift); if(commonThreshold>0.0) { tmpThreshold=commonThreshold; }
musrCounter* counter = new musrCounter(tmpChannel,tmpName,tmpType[0],tmpThreshold,tmpTimeShift, detRecoveryTime / tdcresolution, tmp_keep);
if (tmpType[0]=='M') {mCounter =counter; allCounterMap[tmpChannel]=counter; overallBinDelay=tmpTimeShift;} if (tmpType[0]=='M') {mCounter =counter; allCounterMap[tmpChannel]=counter; overallBinDelay=tmpTimeShift;}
else if (tmpType[0]=='P') {pCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} else if (tmpType[0]=='P') {pCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;}
else if (tmpType[0]=='K') {kCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} else if (tmpType[0]=='K') {kCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;}
@ -852,8 +898,9 @@ void musrAnalysis::CreateHistograms() {
numberOfRewinds=0; numberOfRewinds=0;
numberOfGoodMuons=0; numberOfGoodMuons=0;
lastPreprocessedEntry=-1; lastPreprocessedEntry=-1;
firstPreprocessedEntry=0x7FFFFFFFFFFFFFFFLL;
// currentEventID =0; // currentEventID =0;
fChain->GetEntry(1); Double_t fieldValue = fieldNomVal[0]; if (nFieldNomVal>1) fieldValue += fieldNomVal[1]; fChain->GetEntry(1); InitialiseEvent(); Double_t fieldValue = fieldNomVal[0]; if (nFieldNomVal>1) fieldValue += fieldNomVal[1];
// omega = 851.372*fabs(fieldValue); // omega = 851.372*fabs(fieldValue);
// omega = 851.372*fieldValue; // value set originally in this program // omega = 851.372*fieldValue; // value set originally in this program
// omega = 850.62*fieldValue; // value used in Geant ? // omega = 850.62*fieldValue; // value used in Geant ?
@ -863,6 +910,8 @@ void musrAnalysis::CreateHistograms() {
hInfo->Fill(6, runID); hInfo->Fill(6, runID);
hInfo->Fill(7, hGeantParameters->GetBinContent(7)); hInfo->Fill(7, hGeantParameters->GetBinContent(7));
std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<<fieldValue<<"T, omega="<<omega<<std::endl; std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<<fieldValue<<"T, omega="<<omega<<std::endl;
// zero the condition counters - could use a histogram for this as well?
for (Int_t i=0; i<nrConditions; i++) { conditionCounter[i]=0; }
} }
//================================================================ //================================================================
@ -878,8 +927,10 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
for(std::list<musrTH*>::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { for(std::list<musrTH*>::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) {
(*it)->DrawTH2DdrawList("hist"); (*it)->DrawTH2DdrawList("hist");
} }
for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) { if(musrMode == 'D') {
(it->second)->DrawTDChistogram(); for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
(it->second)->DrawTDChistogram();
}
} }
if (humanDecayHistograms) { if (humanDecayHistograms) {
// humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,indexHuman,nHumanDecayArray); // humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,indexHuman,nHumanDecayArray);
@ -889,14 +940,21 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
humanDecayPileupHistograms->DrawTH1DdrawList("text"); humanDecayPileupHistograms->DrawTH1DdrawList("text");
} }
Long64_t numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates(); // need to re-define these for Pulsed acquisition (and anything else)
Long64_t numberOfMuonCandidatesAfterVK = mCounter->GetNumberOfMuonCandidatesAfterVK(); Long64_t numberOfMuonCandidates;
Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval = mCounter->GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval(); Long64_t numberOfMuonCandidatesAfterVK;
Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000; Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval;
Double_t durationOfExperiment;
if(musrMode=='D') {
numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates();
numberOfMuonCandidatesAfterVK = mCounter->GetNumberOfMuonCandidatesAfterVK();
numberOfMuonCandidatesAfterVKandDoubleHitRemoval = mCounter->GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval();
durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
hInfo->Fill(5,(Double_t) numberOfMuonCandidates); // number of "triggers" hInfo->Fill(5,(Double_t) numberOfMuonCandidates); // number of "triggers"
hInfo->Fill(4,(Double_t) numberOfMuonCandidatesAfterVK); // number of "triggers" hInfo->Fill(4,(Double_t) numberOfMuonCandidatesAfterVK); // number of "triggers"
hInfo->Fill(3,durationOfExperiment); hInfo->Fill(3,durationOfExperiment);
hInfo->Fill(8,(Double_t) numberOfMuonCandidatesAfterVKandDoubleHitRemoval); hInfo->Fill(8,(Double_t) numberOfMuonCandidatesAfterVKandDoubleHitRemoval);
}
//============================== //==============================
// Write out the histograms // Write out the histograms
std::cout<<"musrAnalysis::SaveHistograms()"<<std::endl; std::cout<<"musrAnalysis::SaveHistograms()"<<std::endl;
@ -921,9 +979,15 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
} }
} }
//============================== // condition counts
std::cout << "Number of times each condition was recorded" << std::endl;
for(int i=0; i<nrConditions; i++) {
std::cout << i << " " << conditionCounter[i] << " " << std::endl;
}
//==============================
std::cout<<"==================================================================="<<std::endl; std::cout<<"==================================================================="<<std::endl;
if(musrMode=='D') {
std::cout<<"Number of raw trigger counts (ignoring pileup gate, vetos and coincidences): "<<numberOfMuonCandidates<<std::endl; std::cout<<"Number of raw trigger counts (ignoring pileup gate, vetos and coincidences): "<<numberOfMuonCandidates<<std::endl;
std::cout<<"Number of trigger counts (after vetos and coincidences but not pile-up rejected): "<<numberOfMuonCandidatesAfterVK<<std::endl; std::cout<<"Number of trigger counts (after vetos and coincidences but not pile-up rejected): "<<numberOfMuonCandidatesAfterVK<<std::endl;
std::cout<<"Number of triggered events (i.e. only good \"muons\"): "<<numberOfGoodMuons<<std::endl; std::cout<<"Number of triggered events (i.e. only good \"muons\"): "<<numberOfGoodMuons<<std::endl;
@ -936,13 +1000,29 @@ void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
std::cout<<" Trigger rate after V & K = "<<numberOfMuonCandidatesAfterVK/durationOfExperiment<<" muons/second"; std::cout<<" Trigger rate after V & K = "<<numberOfMuonCandidatesAfterVK/durationOfExperiment<<" muons/second";
std::cout<<" (to get 30000 events/second, set MUONRATEFACTOR = " std::cout<<" (to get 30000 events/second, set MUONRATEFACTOR = "
<< muonRateFactor*( (numberOfMuonCandidatesAfterVK/durationOfExperiment)/30000 )<<")"<<std::endl; << muonRateFactor*( (numberOfMuonCandidatesAfterVK/durationOfExperiment)/30000 )<<")"<<std::endl;
} else if (musrMode=='P') {
Double_t duration_tmp;
if(numberOfRewinds==0) {
duration_tmp = nextUnfilledEventTime;
std::cout <<"One frame only, partial filling = " <<nextUnfilledEventTime/frameInterval << std::endl;
} else if(nextUnfilledEventTime>0 && (doPartialFrameAtEnd)) {
duration_tmp = numberOfRewinds*frameInterval + nextUnfilledEventTime;
std::cout <<"Number of frames: "<<numberOfRewinds<<" whole frames and a " <<nextUnfilledEventTime/frameInterval << "partial one" << std::endl;
} else {
duration_tmp = (numberOfRewinds)*frameInterval;
std::cout << "Number of whole frames: "<<numberOfRewinds<<std::endl;
}
std::cout <<"Number of muon counts (inc. doubles, minus dead): " << numberOfGoodMuons << std::endl;
std::cout << "Duration of the \"experiment\": " << duration_tmp/1000000.0 <<" seconds"<<std::endl;
std::cout << "Displayed count rate: " << numberOfGoodMuons / (duration_tmp/1000000.0) * 3600/1.E6 <<" Mevents/hour"<<std::endl;
}
std::cout<<"========================== E N D =================================="<<std::endl; std::cout<<"========================== E N D =================================="<<std::endl;
} }
//================================================================ //================================================================
void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) { void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
Bool_t pulsedFinishing = false;
// std::cout<<"\n______________________________________________________________________________________________\n"; // std::cout<<"\n______________________________________________________________________________________________\n";
// std::cout<<"musrAnalysis::AnalyseEvent("<<iiiEntry<<") currentTime="<<currentTime<<" (bin="<<Long64_t(currentTime/tdcresolution)<<")"<<std::endl; // std::cout<<"musrAnalysis::AnalyseEvent("<<iiiEntry<<") currentTime="<<currentTime<<" (bin="<<Long64_t(currentTime/tdcresolution)<<")"<<std::endl;
@ -951,13 +1031,27 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
// Loop over several next event and preprocess them (i.e. fill // Loop over several next event and preprocess them (i.e. fill
// them into the lists/maps of the class musrCounter). // them into the lists/maps of the class musrCounter).
// JSL: pulsed mode, loop over and preprocess a whole frame worth of events
if (bool_debugingRequired) { if (bool_debugingRequired) {
if (debugEventMap[eventID]>0) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________ \"Preprocessing Event\"_________"<<std::endl;} if (debugEventMap[eventID]>0) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________ \"Preprocessing Event\"_________"<<std::endl;}
} }
if(musrMode=='P') {
while (((iiiEntry>lastPreprocessedEntry)||((nextUnfilledEventTime<frameInterval)&&(!boolInfinitelyLowMuonRate))) && (lastPreprocessedEntry+1<nentries)) {
// std::cout << "preprocessing event no. " << lastPreprocessedEntry+1 << std::endl;
Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
nextUnfilledEventTime+=deltaT;
};
// avoid a partial frame at the end with lower than usual dead time effects
// but process everything if less than one frame worth of muons to work with
pulsedFinishing = ((lastPreprocessedEntry+1>=nentries) && (nextUnfilledEventTime<frameInterval) ); // flag it for later
// if(pulsedFinishing) { std::cout << "set Finishing flag" << std::endl; }
} else { // Continuous 'D'
while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) { while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) {
Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1); Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
nextUnfilledEventTime+=deltaT; nextUnfilledEventTime+=deltaT;
}; };
};
fChain->GetEntry(iiiEntry); InitialiseEvent(); fChain->GetEntry(iiiEntry); InitialiseEvent();
if (bool_debugingRequired) { if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) PrintHitsInAllCounters(); if (debugEventMap[eventID]>2) PrintHitsInAllCounters();
@ -980,7 +1074,15 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
if (bool_debugingRequired) { if (bool_debugingRequired) {
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(\"Filling Histograms\"_________"<<std::endl;} if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(\"Filling Histograms\"_________"<<std::endl;}
} }
FillHistograms(iiiEntry); if(musrMode=='P') {
// no longer use Not-process-in-bulk option...
//if(!processInBulk && (!pulsedFinishing || boolInfinitelyLowMuonRate || doPartialFrameAtEnd || (numberOfRewinds==0))) {
// std::cout<<"FillHistogramsPulsed("<<iiiEntry<<")"<<std::endl;
// FillHistogramsPulsed(iiiEntry,iiiEntry);
// }
} else {
FillHistograms(iiiEntry);
}
if (bool_debugingRequired) { if (bool_debugingRequired) {
// if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"FillHistograms\"_________"<<std::endl;} // if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"FillHistograms\"_________"<<std::endl;}
if (debugEventMap[eventID]>1) {std::cout<<"____________________________________________________________________"<<std::endl;} if (debugEventMap[eventID]>1) {std::cout<<"____________________________________________________________________"<<std::endl;}
@ -988,14 +1090,36 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
if(musrMode=='P') {
if(iiiEntry==lastPreprocessedEntry) {
if(!pulsedFinishing || doPartialFrameAtEnd || (numberOfRewinds==0) || boolInfinitelyLowMuonRate) {
// std::cout<<"FillHistogramsPulsed(-1000)"<<std::endl;
FillHistogramsPulsed(firstPreprocessedEntry,lastPreprocessedEntry);
}
if(pulsedFinishing && !doPartialFrameAtEnd && (numberOfRewinds>0)) {
nextUnfilledEventTime=0.0;
// std::cout << "resetting NextUnfilledEventTime"<<std::endl;
} else {
if(!pulsedFinishing) {
// std::cout << "Rewinding/clearing Frame" << std::endl;
Long64_t currentTimeBin = Long64_t(frameInterval/tdcresolution);
RemoveOldHitsFromCounters(currentTimeBin-1); // Remove all hits as they can not influence the next frame
firstPreprocessedEntry=0x7FFFFFFFFFFFFFFFLL; // reset
nextUnfilledEventTime -= frameInterval; // only time variable to rewind, otherwise working from start of "this" frame
numberOfRewinds++;
}
};
};
} else { // 'D' - continuous
currentTime+=timeToNextEvent*muonRateFactor; currentTime+=timeToNextEvent*muonRateFactor;
Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution); Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution);
if (!boolInfinitelyLowMuonRate) RemoveOldHitsFromCounters(currentTimeBin-1); // Remove all hits that can not influence the next event if (!boolInfinitelyLowMuonRate ) RemoveOldHitsFromCounters(currentTimeBin-1); // Remove all hits that can not influence the next event
else RemoveOldHitsFromCounters(std::numeric_limits<Long64_t>::max()/2); else RemoveOldHitsFromCounters(std::numeric_limits<Long64_t>::max()/2);
if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere; if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere;
// RewindAllTimeInfo(rewindTimeBins); // RewindAllTimeInfo(rewindTimeBins);
RewindAllTimeInfo(); RewindAllTimeInfo();
} }
}
} }
//================================================================ //================================================================
@ -1044,7 +1168,6 @@ void musrAnalysis::PrintHitsInAllCounters() {
void musrAnalysis::FillHistograms(Int_t iiiEntry) { void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl; // std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl;
oncePerEvent = true; oncePerEvent = true;
Bool_t mCounterHitExistsForThisEventID = false; Bool_t mCounterHitExistsForThisEventID = false;
Bool_t pCounterHitExistsForThisEventID = false; Bool_t pCounterHitExistsForThisEventID = false;
@ -1126,7 +1249,7 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// //
if (pCounterHitExistsForThisEventID&&(posEntry>0)) { if (pCounterHitExistsForThisEventID&&(posEntry>0)) {
// Get variables of the detector hit corresponding to the positron candidate: // Get variables of the detector hit corresponding to the positron candidate:
fChain->GetEntry(posEntry); fChain->GetEntry(posEntry); InitialiseEvent();
Double_t detP_edep = det_edep[idetP]; Double_t detP_edep = det_edep[idetP];
if (detP_edep!=idetP_edep) {std::cout<<"KIKS: detP_edep="<<detP_edep<<"; idetP_edep="<<idetP_edep<<"; idetP= "<<idetP<<std::endl; exit(1);} if (detP_edep!=idetP_edep) {std::cout<<"KIKS: detP_edep="<<detP_edep<<"; idetP_edep="<<idetP_edep<<"; idetP= "<<idetP<<std::endl; exit(1);}
detP_x = det_x[idetP]; detP_x = det_x[idetP];
@ -1148,7 +1271,7 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
pileup_muDecayPosR = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY); pileup_muDecayPosR = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
// fChain->GetEntry(iiiEntry); // fChain->GetEntry(iiiEntry);
} }
fChain->GetEntry(iiiEntry); fChain->GetEntry(iiiEntry); InitialiseEvent();
} }
} }
// if (doubleHitP) { // if (doubleHitP) {
@ -1184,7 +1307,6 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// det_time0 = timeBin0*tdcresolution; // det_time0 = timeBin0*tdcresolution;
// det_time1 = timeBin1*tdcresolution; // det_time1 = timeBin1*tdcresolution;
det_time10 = ((timeBin0==-100000000)|| (timeBin1==-100000000)) ? -100000000 : (timeBin1-timeBin0)*tdcresolution; det_time10 = ((timeBin0==-100000000)|| (timeBin1==-100000000)) ? -100000000 : (timeBin1-timeBin0)*tdcresolution;
//PROMPT PEAK INVESTIGATION: if (timeBin1-timeBin0==7) std::cout<<"HHHHHHHHHHHH eventID="<<eventID<<", t0="<<timeBin0*tdcresolution<<std::endl;
det_time20 = ((timeBin0==-100000000)|| (timeBin2==-100000000)) ? -100000000 : (timeBin2-timeBin0)*tdcresolution; det_time20 = ((timeBin0==-100000000)|| (timeBin2==-100000000)) ? -100000000 : (timeBin2-timeBin0)*tdcresolution;
//cks if ((det_time10>-1)&&(det_time10<0.01)) {std::cout<<" eventID="<<eventID<<" det_time10="<< det_time10<<std::endl;} //cks if ((det_time10>-1)&&(det_time10<0.01)) {std::cout<<" eventID="<<eventID<<" det_time10="<< det_time10<<std::endl;}
gen_time10 = muDecayTime-muM0Time; gen_time10 = muDecayTime-muM0Time;
@ -1204,6 +1326,7 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// Find whether the muon decayed in the sample; the sample may consist of more than just one volume. The ID // Find whether the muon decayed in the sample; the sample may consist of more than just one volume. The ID
// of the sample volume (or volumes) is defined in the *.v1190 file after the keyword "sampleID": // of the sample volume (or volumes) is defined in the *.v1190 file after the keyword "sampleID":
muonDecayedInSample_gen = (find(SampleDetIDList.begin(), SampleDetIDList.end(), muDecayDetID)) != SampleDetIDList.end() ; muonDecayedInSample_gen = (find(SampleDetIDList.begin(), SampleDetIDList.end(), muDecayDetID)) != SampleDetIDList.end() ;
muonDecayedInSampleOnce_gen = muonDecayedInSample_gen && oncePerEvent;
muonTriggered_gen = muM0Time > -1000; muonTriggered_gen = muM0Time > -1000;
muonTriggered_gen_AND_muonDecayedInSample_gen = muonTriggered_gen && muonDecayedInSample_gen; muonTriggered_gen_AND_muonDecayedInSample_gen = muonTriggered_gen && muonDecayedInSample_gen;
muonTriggered_det = mCounterHitExistsForThisEventID; muonTriggered_det = mCounterHitExistsForThisEventID;
@ -1230,6 +1353,14 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// debugEventMap.insert(std::pair<int,int>(eventID,10)); // debugEventMap.insert(std::pair<int,int>(eventID,10));
} }
// for better assignment of backgrounds versus useful signal
goodEvent_F_det_AND_muonDecayedInSample_gen = goodEvent_F_det && muonDecayedInSample_gen;
goodEvent_B_det_AND_muonDecayedInSample_gen = goodEvent_B_det && muonDecayedInSample_gen;
goodEvent_U_det_AND_muonDecayedInSample_gen = goodEvent_U_det && muonDecayedInSample_gen;
goodEvent_D_det_AND_muonDecayedInSample_gen = goodEvent_D_det && muonDecayedInSample_gen;
goodEvent_L_det_AND_muonDecayedInSample_gen = goodEvent_L_det && muonDecayedInSample_gen;
goodEvent_R_det_AND_muonDecayedInSample_gen = goodEvent_R_det && muonDecayedInSample_gen;
goodEvent_F_det_AND_pileupEvent = goodEvent_F_det && pileupEvent; goodEvent_F_det_AND_pileupEvent = goodEvent_F_det && pileupEvent;
goodEvent_B_det_AND_pileupEvent = goodEvent_B_det && pileupEvent; goodEvent_B_det_AND_pileupEvent = goodEvent_B_det && pileupEvent;
goodEvent_U_det_AND_pileupEvent = goodEvent_U_det && pileupEvent; goodEvent_U_det_AND_pileupEvent = goodEvent_U_det && pileupEvent;
@ -1269,6 +1400,7 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) { for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
int i = it->first; int i = it->first;
condition[i]=*(it->second); condition[i]=*(it->second);
if(condition[i]) { conditionCounter[i]++; }
} }
// A special condition can me added here: // A special condition can me added here:
@ -1306,6 +1438,429 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
} }
//================================================================
void musrAnalysis::FillHistogramsPulsed(Int_t iiiEntry1,Int_t iiiEntry2) {
// std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl;
// edited for Pulsed Mode acquisition rather than fill "FillHistograms" with lots of if() statements
// may also be useful for Integral Mode acquisition?
// search for events iiiEntry1 to iiiEntry2 (maybe whole frame)
// process "once per event" for all, regardless of whether there are hits
//
// det_time_10: time difference between "mean t0" and the positron counter hit (=measured), will include the muon's time of flight to the sample
// gen_time_20: time from the muon entering the target to the positron event (-1000 for decay in flight, background, etc)
// det_time_10_MINUS_gen_time_10 = position in pulse, after initial arrival transport
// det_time_20 and gen_time_20 adjusted for phase shifts
// pileup* not applicable
oncePerEvent = true; // each incoming muon event, one randomly chosen detector hit if there are any at all but probably from a low numbered counter
// Bool_t mCounterHitExistsForThisEventID = false;
Int_t positronQuality = -1000;
Long64_t timeBinOfThePreviouslyProcessedHit = -100000000;
Long64_t timeBin0 = -100000000;
Long64_t timeBin1 = -100000000;
Long64_t timeBin2 = -100000000;
Long64_t timeBinDoubleHit;
Int_t kEntry = 0;
Int_t posEntry = 0;
Int_t idetM = 0;
Int_t idetM_ID = 0;
Double_t idetM_edep = 0.;
Int_t idetP = 0;
Int_t idetP_ID = 0;
Double_t idetP_edep = 0.;
Int_t idetFirstMulti = -1000;
Int_t bank_this = 0;
Int_t bank_multi = 0;
//cks 4.10.2011 Bool_t doubleHitM = false;
//cks 4.10.2011 Bool_t doubleHitP = false;
// typedef std::map<Int_t,Bool_t> unrecordedEventMap;
// unrecordedEventMap unrecordedEvents;
// unrecordedEventMap::const_iterator it2;
typedef std::multimap<Long64_t,hitInfo*> sortedByEventMap;
sortedByEventMap sortedByEvent;
sortedByEventMap::const_iterator it3;
hitInfo* theEvent;
// new improved
// start with empty multimap of hits (key=entry num)
// loop through counters
// scan time order
// "copy" events to map
// initialise iterator into sorted map
// loop through event numbers in range
// if events to be processed then do them (1st=OnceOnly)
// else do OnceOnly with no hit
//
// hitInfo needs PositronQuality field (filled in during GetNextGoodPositronPulsed() and timeBin1 )
// GetNextGoodPositronPulsed() now returns ptr to the hitInfo (or NULL)..?
// std::cout << "starting FHP; map contents=" << sortedByEvent.size() << std::endl;
for (counterMapType::const_iterator it = pCounterMap.begin(); it != pCounterMap.end(); ++it) {
Long64_t dataBinMin = dataWindowBinMin;
Long64_t dataBinMax = dataWindowBinMax;
timeBin1 = -100000000;
timeBin2 = -100000000;
do {
positronQuality = (it->second)->GetNextGoodPositronPulsed(iiiEntry1,iiiEntry2,dataBinMin,dataBinMax,timeBin1,timeBin2,posEntry,idetP,idetP_ID,idetP_edep,idetFirstMulti,theEvent);
if(positronQuality != -1000) {
// std::cout << "pushing event " << posEntry << " det " << idetP_ID << std::endl;
sortedByEvent.insert(std::pair<Long64_t,hitInfo*>(posEntry,theEvent));
}
} while (positronQuality != -1000);
}
it3=sortedByEvent.begin();
// std::cout << "now processing sorted list" << std:: endl;
//for(Int_t i=iiiEntry1; i<=iiiEntry2; i++) {
// unrecordedEvents.insert(std::pair<Int_t,Bool_t>(i,true));
//}
// std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
// mCounterHitExistsForThisEventID = MuonCounterHit(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep);
// not needed (JSL) mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep);
// timeBinOfThePreviouslyProcessedHit = timeBin0;
//___________________________________________________________
// loop through all good P counters and then one extra go with it==pCounterMap.end() to tidy up unrecorded events
// for (counterMapType::const_iterator it = pCounterMap.begin(); ; ++it) {
for(Int_t iec=iiiEntry1; iec<=iiiEntry2; iec++) {
oncePerEvent=true;
fChain->GetEntry(iec); InitialiseEvent(); // now opening each event once only and in order! (as well as once for inserting hits in counters)
// positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep);
// std::cout<<" timeBin0 ="<<timeBin0<<std::endl;
// Check whether there was good hit in the Positron counter
// Long64_t dataBinMin = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMin : timeBinOfThePreviouslyProcessedHit-100000000;
// Long64_t dataBinMax = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMax : timeBinOfThePreviouslyProcessedHit+100000000;
detP_x = -1001.;
detP_y = -1001.;
detP_z = -1001.;
detP_time_start = -1001.;
detP_time_end = -1001.;
detP_theta = -1001.;
detP_phi = -1001.;
detP_phi_MINUS_pos_Phi = -1001.;
detP_phi_MINUS_pos_Phi360 = -1001.;
detP_theta_MINUS_pos_Theta = -1001.;
detP_theta_MINUS_pos_Theta360 = -1001.;
detP_time_start_MINUS_muDecayTime = -1001.;
pileup_eventID = -1001;
pileup_muDecayDetID_double = -1001;
pileup_muDecayPosX = -1000000000;
pileup_muDecayPosY = -1000000000;
pileup_muDecayPosZ = -1000000000;
pileup_muDecayPosR = -1000000000;
det_time31 = -1001.;
timeBinDoubleHit = -100000000;
doubleHit = false;
pos_detID = -1.;
pos_detID_doubleHit = -1.;
Int_t pos_detID_doubleHit_INT = -1;
timeBin1 = -100000000;
timeBin2 = -100000000;
Long64_t dataBinMin = dataWindowBinMin;
Long64_t dataBinMax = dataWindowBinMax;
// Long64_t positronBinMax = timeBin0+pileupWindowBinMax+dataWindowBinMin;
Long64_t positronBinMax = dataBinMax-dataWindowBinMin; // note that "dataWindowBinMin" is normally negative, i.e. positronBinMax > dataBinMax !!!
// GetNextGoodPositronPulsed(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep) {
while (oncePerEvent || positronQuality != -1000) {
// does the next hit in sortedByEvent belong to this?
if(it3!=sortedByEvent.end() && it3->first == iec) {
positronQuality=it3->second->posQual;
timeBin1=it3->second->timeBin1;
timeBin2=it3->second->timeBin2;
posEntry=it3->second->eventEntry;
idetP=it3->second->det_i;
idetP_ID=it3->second->det_id;
idetP_edep=it3->second->det_edep;
if(((it3->second)->multiHitCtr > 1) && (it3->second)->firstMulti) {
idetFirstMulti = (it3->second)->firstMulti->det_i;
} else {
idetFirstMulti = -1000; // not a second or subsequent multi hit, or not recorded for some reason
}
// std::cout << "processing hit in det " << idetP_ID << " for event " << iec << std::endl;
it3++;
} else {
if(!oncePerEvent) {
// std::cout << "run out of hits for event " << iec << std::endl;
break;
}
positronQuality=-1000;
posEntry=iec;
timeBin1=-100000000;
timeBin2=-100000000;
idetP=-1000;
idetP_ID=-1000;
idetP_edep=-1000;
idetFirstMulti=-1000;
// std::cout << "processing non-hit for event " << iec << std::endl;
}
// if(it == pCounterMap.end()) {
// while((it2 != unrecordedEvents.end()) && (!it2->second)) { it2++; }
// if(it2 == unrecordedEvents.end()) { break; }
// positronQuality=-1000;
// posEntry=it2->first;
// timeBin1=-1000;
// timeBin2=-1000;
// idetP=-1000;
// idetP_ID=-1000;
// idetP_edep=-1000;
// idetFirstMulti=-1000;
// // std::cout << "found entry " << posEntry << " with no hits" << std::endl;
// } else {
// if((positronQuality = (it->second)->GetNextGoodPositronPulsed(iiiEntry1,iiiEntry2,dataBinMin,dataBinMax,timeBin1,timeBin2,posEntry,idetP,idetP_ID,idetP_edep,idetFirstMulti)) == -1000 ) { break; }
// // std::cout << "found a hit in detector " << idetP_ID << " from event " << posEntry << " with quality=" << positronQuality << std::endl;
// }
// // std::cout << "searching entry " << iiiEntry << " and found a hit in detector " << idetP_ID << " with quality=" << positronQuality << std::endl;
// fChain->GetEntry(posEntry); InitialiseEvent();
// oncePerEvent=unrecordedEvents[posEntry];
// unrecordedEvents[posEntry] = false;
if (muTargetTime != -1000.0) {
timeBin0 = Long64_t(muTargetTime / tdcresolution);
} else {
timeBin0=-100000000;
}
pos_detID = Double_t(idetP_ID);
if (debugEventMap[eventID]>2) {
std::cout<<"FillHistograms: GOOD positron candidate found: timeBin1="<<timeBin1
<<", deltat+800="<<timeBin1-timeBin0+800 <<std::endl;}
// Get variables of the detector hit corresponding to the positron candidate:
if(positronQuality != -1000) {
Double_t detP_edep = det_edep[idetP];
if (detP_edep!=idetP_edep) {std::cout<<"KIKS: detP_edep="<<detP_edep<<"; idetP_edep="<<idetP_edep<<"; idetP= "<<idetP<<std::endl; exit(1);}
detP_x = det_x[idetP];
detP_y = det_y[idetP];
detP_z = det_z[idetP];
detP_time_start = det_time_start[idetP];
detP_time_end = det_time_end[idetP];
detP_phi = myAtan2(detP_y,detP_x) * 180./pi;
detP_theta = myAtan2(sqrt(detP_y*detP_y+detP_x+detP_x),detP_z) * 180./pi;
} else {
detP_x = -1000;
detP_y = -1000;
detP_z = -1000;
detP_time_start = -1000;
detP_time_end = -1000;
detP_phi = -1000;
detP_theta = -1000;
}
//
// if (doubleHitP) {
// std::cout<<"yyyyyyyyyyyyyyyyy Total double hit nr="<<testIVar1++<<std::endl;
// }
// CALCULATE VARIABLES
// Initial muon
muIniPosR = sqrt(muIniPosX*muIniPosX+muIniPosY*muIniPosY);
muIniMomTrans = sqrt(muIniMomX*muIniMomX+muIniMomY*muIniMomY);
muTargetPol_Theta = myAtan2(sqrt(muTargetPolX*muTargetPolX+muTargetPolY*muTargetPolY),muTargetPolZ) * 180./pi;
muTargetPol_Theta360 = (muTargetPol_Theta<0) ? muTargetPol_Theta+360. : muTargetPol_Theta;
muTargetPol_Phi = myAtan2(muTargetPolY,muTargetPolX) * 180./pi;
muTargetPol_Phi360= (muTargetPol_Phi<0) ? muTargetPol_Phi+360. : muTargetPol_Phi;
muDecayPol_Theta = myAtan2(sqrt(muDecayPolX*muDecayPolX+muDecayPolY*muDecayPolY),muDecayPolZ) * 180./pi;
muDecayPol_Theta360 = (muDecayPol_Theta<0) ? muDecayPol_Theta+360. : muDecayPol_Theta;
muDecayPol_Phi = myAtan2(muDecayPolY,muDecayPolX) * 180./pi;
muDecayPol_Phi360= (muDecayPol_Phi<0) ? muDecayPol_Phi+360. : muDecayPol_Phi;
// Initial positron
pos_Trans_Momentum = sqrt(posIniMomX*posIniMomX+posIniMomY*posIniMomY);
pos_Momentum = sqrt(pos_Trans_Momentum*pos_Trans_Momentum+posIniMomZ*posIniMomZ);
pos_Radius = pos_Trans_Momentum/(-BFieldAtDecay_Bz)/0.3;
pos_Theta = myAtan2(pos_Trans_Momentum,posIniMomZ) * 180./pi;
pos_Theta360 = (pos_Theta<0) ? pos_Theta+360. : pos_Theta;
pos_Phi = myAtan2(posIniMomY,posIniMomX) * 180./pi;
pos_Phi360 = (pos_Phi<0) ? pos_Phi+360. : pos_Phi;
pos_Theta_MINUS_muDecayPol_Theta = deltaAngle(pos_Theta,muDecayPol_Theta);
pos_Theta_MINUS_muDecayPol_Theta360 = (pos_Theta_MINUS_muDecayPol_Theta<0) ? pos_Theta_MINUS_muDecayPol_Theta+360 : pos_Theta_MINUS_muDecayPol_Theta;
pos_Phi_MINUS_muDecayPol_Phi = deltaAngle(pos_Phi,muDecayPol_Phi);
pos_Phi_MINUS_muDecayPol_Phi360 = (pos_Phi_MINUS_muDecayPol_Phi<0) ? pos_Phi_MINUS_muDecayPol_Phi+360 : pos_Phi_MINUS_muDecayPol_Phi;
// Detector info
// det_m0edep = (mCounterHitExistsForThisEventID) ? idetM_edep : -1000.;
// det_time0 = timeBin0*tdcresolution;
// det_time1 = timeBin1*tdcresolution;
det_time10 = (timeBin1==-100000000) ? -1000.0 : timeBin1*tdcresolution;
det_time20 = (timeBin2==-100000000) ? -1000.0 : timeBin2*tdcresolution;
//cks if ((det_time10>-1)&&(det_time10<0.01)) {std::cout<<" eventID="<<eventID<<" det_time10="<< det_time10<<std::endl;}
if (muTargetTime != -1000.0) {
gen_time10 = muDecayTime-muTargetTime;
det_time10_MINUS_gen_time10 = (det_time10 - gen_time10)/picosecond;
} else {
gen_time10 = -1000.0;
det_time10_MINUS_gen_time10 = -1000.0;
}
det_posEdep = idetP_edep ;
det_time1_MINUS_muDecayTime = (timeBin1*tdcresolution-muDecayTime)/picosecond;
if ((detP_time_start>-998.)&&(muDecayTime>-998.)) detP_time_start_MINUS_muDecayTime = (detP_time_start - muDecayTime)/picosecond;
detP_phi_MINUS_pos_Phi = deltaAngle(detP_phi,pos_Phi);
// std::cout<<"detP_phi_MINUS_pos_Phi="<<detP_phi_MINUS_pos_Phi<<std::endl;
detP_phi_MINUS_pos_Phi360 = (detP_phi_MINUS_pos_Phi<0) ? detP_phi_MINUS_pos_Phi+360 : detP_phi_MINUS_pos_Phi;
detP_theta_MINUS_pos_Theta = deltaAngle(detP_theta,pos_Theta);
detP_theta_MINUS_pos_Theta360 = (detP_theta_MINUS_pos_Theta<0) ? detP_theta_MINUS_pos_Theta+360 : detP_theta_MINUS_pos_Theta;
// std::cout<<eventID<<" det_time10="<<det_time10<<" t1="<<timeBin1*tdcresolution<<" t0="<<timeBin0*tdcresolution<<" gen_time10="<<gen_time10<<std::endl;
// CALCULATE CONDITIONS
alwaysTrue = true;
// Find whether the muon decayed in the sample; the sample may consist of more than just one volume. The ID
// of the sample volume (or volumes) is defined in the *.v1190 file after the keyword "sampleID":
muonDecayedInSample_gen = (find(SampleDetIDList.begin(), SampleDetIDList.end(), muDecayDetID)) != SampleDetIDList.end() ;
muonDecayedInSampleOnce_gen = muonDecayedInSample_gen && oncePerEvent;
// muonTriggered_gen = muM0Time > -1000;
// muonTriggered_gen_AND_muonDecayedInSample_gen = muonTriggered_gen && muonDecayedInSample_gen;
// muonTriggered_det = mCounterHitExistsForThisEventID;
// positronHit_det = pCounterHitExistsForThisEventID;
goodEvent_det = (positronQuality >= -1);
if(goodEvent_det) numberOfGoodMuons++;
goodEvent_gen = (positronQuality== -3 || positronQuality== -2 || positronQuality==0 || positronQuality==1);
goodEvent_det_AND_goodEvent_gen = goodEvent_det && goodEvent_gen;
// pileupEventCandidate = ((kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) ? true:false;
// pileupEvent = pileupEventCandidate&&goodEvent_det;
goodEvent_det_AND_muonDecayedInSample_gen = goodEvent_det && muonDecayedInSample_gen;
// posCounterList_Iterator = find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID);
// goodEvent_F_det = posCounterList_Iterator != F_posCounterList.end()
Bool_t event_F = ( (find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID)) != F_posCounterList.end() );
Bool_t event_B = ( (find(B_posCounterList.begin(), B_posCounterList.end(), idetP_ID)) != B_posCounterList.end() );
Bool_t event_U = ( (find(U_posCounterList.begin(), U_posCounterList.end(), idetP_ID)) != U_posCounterList.end() );
Bool_t event_D = ( (find(D_posCounterList.begin(), D_posCounterList.end(), idetP_ID)) != D_posCounterList.end() );
Bool_t event_L = ( (find(L_posCounterList.begin(), L_posCounterList.end(), idetP_ID)) != L_posCounterList.end() );
Bool_t event_R = ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() );
bank_this = (event_F ? 1:0) + (event_B ? 2:0) + (event_U ? 3:0) + (event_D ? 4:0) + (event_L ? 5:0) + (event_R ? 6:0) ;
if(idetFirstMulti >= 0) {
det_multi_interval = det_time_start[idetP]-det_time_start[idetFirstMulti];
pos_detID_doubleHit = det_ID[idetFirstMulti];
// pos_doubleHit_dPhi: difference in detector phases between this and the first multi hit, plus 10*first-bank + 60*second-bank
bank_multi = (( (find(F_posCounterList.begin(), F_posCounterList.end(), det_ID[idetFirstMulti])) != F_posCounterList.end() ) ? 1:0) +
(( (find(B_posCounterList.begin(), B_posCounterList.end(), det_ID[idetFirstMulti])) != B_posCounterList.end() ) ? 2:0) +
(( (find(U_posCounterList.begin(), U_posCounterList.end(), det_ID[idetFirstMulti])) != U_posCounterList.end() ) ? 3:0) +
(( (find(D_posCounterList.begin(), D_posCounterList.end(), det_ID[idetFirstMulti])) != D_posCounterList.end() ) ? 4:0) +
(( (find(L_posCounterList.begin(), L_posCounterList.end(), det_ID[idetFirstMulti])) != L_posCounterList.end() ) ? 5:0) +
(( (find(R_posCounterList.begin(), R_posCounterList.end(), det_ID[idetFirstMulti])) != R_posCounterList.end() ) ? 6:0);
pos_doubleHit_dPhi= deltaAngle(phaseShiftMap[det_ID[idetFirstMulti]]*180/pi,phaseShiftMap[idetP_ID]*180/pi) + bank_this*1000 + bank_multi*6000;
} else {
det_multi_interval = -1000;
pos_detID_doubleHit = -1000;
pos_doubleHit_dPhi= -1000;
}
goodEvent_F_det = goodEvent_det && event_F;
goodEvent_B_det = goodEvent_det && event_B;
goodEvent_U_det = goodEvent_det && event_U;
goodEvent_D_det = goodEvent_det && event_D;
goodEvent_L_det = goodEvent_det && event_L;
goodEvent_R_det = goodEvent_det && event_R;
// for better assignment of backgrounds versus useful signal
goodEvent_F_det_AND_muonDecayedInSample_gen = goodEvent_F_det && muonDecayedInSample_gen;
goodEvent_B_det_AND_muonDecayedInSample_gen = goodEvent_B_det && muonDecayedInSample_gen;
goodEvent_U_det_AND_muonDecayedInSample_gen = goodEvent_U_det && muonDecayedInSample_gen;
goodEvent_D_det_AND_muonDecayedInSample_gen = goodEvent_D_det && muonDecayedInSample_gen;
goodEvent_L_det_AND_muonDecayedInSample_gen = goodEvent_L_det && muonDecayedInSample_gen;
goodEvent_R_det_AND_muonDecayedInSample_gen = goodEvent_R_det && muonDecayedInSample_gen;
// std::cout<<"goodEvent_F_det="<<goodEvent_F_det<<std::endl;
promptPeak = goodEvent_det && (det_time10>promptPeakWindowMin) && (det_time10<promptPeakWindowMax);
promptPeakF = promptPeak && goodEvent_F_det;
promptPeakB = promptPeak && goodEvent_B_det;
promptPeakU = promptPeak && goodEvent_U_det;
promptPeakD = promptPeak && goodEvent_D_det;
promptPeakL = promptPeak && goodEvent_L_det;
promptPeakR = promptPeak && goodEvent_R_det;
pileupEvent = (positronQuality > -1000 && positronQuality<=-2); // a good positron missed because of dead time
goodEvent_F_det_AND_pileupEvent = pileupEvent && event_F;
goodEvent_B_det_AND_pileupEvent = pileupEvent && event_B;
goodEvent_U_det_AND_pileupEvent = pileupEvent && event_U;
goodEvent_D_det_AND_pileupEvent = pileupEvent && event_D;
goodEvent_L_det_AND_pileupEvent = pileupEvent && event_L;
goodEvent_R_det_AND_pileupEvent = pileupEvent && event_R;
doubleHitEvent_gen = (positronQuality>1 || (positronQuality>-1000 && positronQuality<-3)); // any double hits
doubleHit = (positronQuality>1); // counted double hits only
goodEvent_F_det_AND_doubleHit = doubleHit && event_F;
goodEvent_B_det_AND_doubleHit = doubleHit && event_B;
goodEvent_U_det_AND_doubleHit = doubleHit && event_U;
goodEvent_D_det_AND_doubleHit = doubleHit && event_D;
goodEvent_L_det_AND_doubleHit = doubleHit && event_L;
goodEvent_R_det_AND_doubleHit = doubleHit && event_R;
singleHitEvent_gen = (positronQuality==0);
stackedEvent_gen = (positronQuality == -1);
// if (bool_debugingRequired && muonTriggered_det) {
// std::cout<<"DEBUG: goodEvent_det: eventID="<<eventID<<std::endl;
// if (goodEvent_det) std::cout<<" ___DETECTED___"<<std::endl;
// MyPrintTree();
// MyPrintConditions();
// }
// if (pileupEvent) {
// std::cout<<"\n NEW: pileupEvent: eventID="<<eventID<<", kEntry="<<kEntry<<", posEntry="<<posEntry<< std::endl;
// std::cout<<"det_time10 = "<<det_time10<<std::endl;
// }
// Fill pileup-variables, but only if positron comes from different muon than the trigger signal
// if (mCounterHitExistsForThisEventID&&pCounterHitExistsForThisEventID)
// std::cout<<eventID<<" "<<mCounterHitExistsForThisEventID<<": iiiEntry (trigger) = "<<iiiEntry<<" kEntry (muon hit)="<<kEntry<<" "<<pCounterHitExistsForThisEventID<<" posEntry (positron hit)="<<posEntry<<std::endl;
// if (goodEvent_det) {
// std::cout<<"wwwwwwwwwwwwww \nEvent Nr: "<<eventID<<" muDecayTime="<< muDecayTime <<" muM0Time="<<muM0Time<<" det_time_start[idetP]="<<det_time_start[idetP]<<" t1-t0="<<(timeBin1-timeBin0)*tdcresolution<<std::endl;
// }
for (Int_t i=0; i<nrConditions; i++) { condition[i]=false; }
// MyPrintConditions();
for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
int i = it->first;
condition[i]=*(it->second);
if(condition[i]) { conditionCounter[i]++; }
}
// A special condition can me added here:
// condition[29] = muIniMomZ>27;
// Fill all "Detector Histograms"
for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
(*it)->FillTH1D(wght,&condition[0]);
}
for(std::list<musrTH*>::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) {
(*it)->FillTH2D(wght,&condition[0]);
}
// JSL don't plot the individual histograms - too many of them!
// if ((goodEvent_det)&&(pCounterHitExistsForThisEventID)) {
// musrCounter* pCounterHitInThisEvent
// counterMapType::const_iterator itPcounterHitInThisEvent = pCounterMap.find(idetP_ID);
// if (itPcounterHitInThisEvent==pCounterMap.end()) {
// std::cout<<" ERROR! musrAnalysis: itPcounterHitInThisEvent not found ==> This should never happen!!!"<<std::endl;
// std::cout<<" idetP="<<idetP<<" idetP_ID="<<idetP_ID<< std::endl;
// }
// else{
// (itPcounterHitInThisEvent->second)->FillTDChistogram(float(timeBin1-timeBin0+overallBinDelay)+0.5,wght);
// }
// }
// Check whether there is a good trigger in the next event
// std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
// mCounterHitExistsForThisEventID = MuonCounterHit(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep);
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep);
timeBinOfThePreviouslyProcessedHit = timeBin1;
// if (mCounterHitExistsForThisEventID) std::cout<<" YYYYYYYYYYYYYYYYYYY check this : LOOOPING AGAIN"<<std::endl;
oncePerEvent=false;
}; // loop through possible hits
};// loop through events
}
//================================================================ //================================================================
void musrAnalysis::InitialiseEvent() { void musrAnalysis::InitialiseEvent() {
runID_double = runID; runID_double = runID;
@ -1320,6 +1875,7 @@ void musrAnalysis::InitialiseEvent() {
//================================================================ //================================================================
Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) { Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
// std::cout<<"musrAnalysis::PreprocessEvent()"<<std::endl; // std::cout<<"musrAnalysis::PreprocessEvent()"<<std::endl;
Int_t multiCtr;
fChain->GetEntry(iEn); InitialiseEvent(); fChain->GetEntry(iEn); InitialiseEvent();
// Clone some channels into different one, if requested by user // Clone some channels into different one, if requested by user
@ -1384,13 +1940,27 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
// std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<<eventID<<std::endl; // std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<<eventID<<std::endl;
// MyPrintTree(); // MyPrintTree();
Double_t globTime = nextUnfilledEventTime; Double_t globTime;
if(musrMode=='P') {
globTime=muIniTime * muonPulseWidthFactor; // JSL: use time within pulse
} else { // continuous
globTime= nextUnfilledEventTime;
}
multiCtr=0; // JSL: number the above-threshold Good counts as 0 (single) or 1,2 etc (multi); small ones all -1
// must insert events in time order for correct labelling of multiple hits as "first", "second", etc
// use multimap to sort them. Sorting by GEANT raw time is OK here
doubleHitSorterType doubleHitSorter;
for (Int_t i=0; i<det_n; i++) { for (Int_t i=0; i<det_n; i++) {
doubleHitSorter.insert(std::pair<double,int>(det_time_start[i],i));
}
for(doubleHitSorterType::const_iterator dhi = doubleHitSorter.begin(); dhi != doubleHitSorter.end(); dhi++) {
Int_t i=dhi->second;
// // Int_t detID=det_ID[i]; // // Int_t detID=det_ID[i];
std::map<int,musrCounter*>::iterator it; std::map<int,musrCounter*>::iterator it;
it = allCounterMap.find(det_ID[i]); it = allCounterMap.find(det_ID[i]);
if (it==allCounterMap.end()) { if (it==allCounterMap.end()) {
// std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl; std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
} }
else { else {
// Double_t omega = 851.372*fieldNomVal[0]; // Double_t omega = 851.372*fieldNomVal[0];
@ -1402,13 +1972,15 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
Long64_t timeBin = Long64_t(t1); Long64_t timeBin = Long64_t(t1);
Long64_t timeBin2 = Long64_t(t2); Long64_t timeBin2 = Long64_t(t2);
(*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i],eventID); (*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i],eventID, multiCtr);
// std::cout << "pushed event into counter "<<det_ID[i]<<" resulting in multiCtr="<<multiCtr<<std::endl;
} }
} }
// std::cout<<"lastPreprocessedEntry+1="<<lastPreprocessedEntry+1<<" iEn="<<iEn<<" eventID="<<eventID<<std::endl; // std::cout<<"lastPreprocessedEntry+1="<<lastPreprocessedEntry+1<<" iEn="<<iEn<<" eventID="<<eventID<<std::endl;
if ((lastPreprocessedEntry+1)!=iEn) {std::cout<<"ERROR PreprocessEvent - should never happen!!!"<<std::endl; } if ((lastPreprocessedEntry+1)!=iEn) {std::cout<<"ERROR PreprocessEvent - should never happen!!!"<<std::endl; }
lastPreprocessedEntry = iEn; lastPreprocessedEntry = iEn;
if(iEn<firstPreprocessedEntry) { firstPreprocessedEntry=iEn; }
return timeToNextEvent*muonRateFactor; return timeToNextEvent*muonRateFactor;
} }

67
musrSimAna/musrAnalysis.hh
View File

@ -23,6 +23,12 @@
//#include "musrTH.hh" //#include "musrTH.hh"
#include <math.h> #include <math.h>
const double pi=3.14159265358979324; // initialise here instead JSL
const double microsecond=1.;
const double nanosecond=0.001;
const double picosecond=0.000001;
class musrTH; class musrTH;
//#include "musrSimGlobal.hh" //#include "musrSimGlobal.hh"
@ -192,6 +198,7 @@ public :
Double_t eventID_double; Double_t eventID_double;
Double_t muDecayDetID_double; Double_t muDecayDetID_double;
Double_t det_n_double; Double_t det_n_double;
Double_t det_multifirst_double;
Double_t muDecayPosR; Double_t muDecayPosR;
Double_t wght; Double_t wght;
Double_t det_m0edep; Double_t det_m0edep;
@ -219,6 +226,7 @@ public :
Double_t pos_Phi_MINUS_muDecayPol_Phi360; Double_t pos_Phi_MINUS_muDecayPol_Phi360;
Double_t pos_detID; Double_t pos_detID;
Double_t pos_detID_doubleHit; Double_t pos_detID_doubleHit;
Double_t pos_doubleHit_dPhi;
// Double_t det_time0; // Double_t det_time0;
// Double_t get_time0; // Double_t get_time0;
// Double_t det_time1; // Double_t det_time1;
@ -229,6 +237,7 @@ public :
Double_t det_time20; Double_t det_time20;
Double_t det_time31; Double_t det_time31;
Double_t det_time1_MINUS_muDecayTime; Double_t det_time1_MINUS_muDecayTime;
Double_t det_multi_interval;
Double_t detP_x; Double_t detP_x;
Double_t detP_y; Double_t detP_y;
Double_t detP_z; Double_t detP_z;
@ -254,6 +263,7 @@ public :
Bool_t alwaysTrue; Bool_t alwaysTrue;
Bool_t oncePerEvent; Bool_t oncePerEvent;
Bool_t muonDecayedInSample_gen; Bool_t muonDecayedInSample_gen;
Bool_t muonDecayedInSampleOnce_gen;
Bool_t muonTriggered_gen; Bool_t muonTriggered_gen;
Bool_t muonTriggered_gen_AND_muonDecayedInSample_gen; Bool_t muonTriggered_gen_AND_muonDecayedInSample_gen;
Bool_t muonTriggered_det; Bool_t muonTriggered_det;
@ -264,6 +274,12 @@ public :
Bool_t pileupEventCandidate; Bool_t pileupEventCandidate;
Bool_t pileupEvent; Bool_t pileupEvent;
Bool_t goodEvent_det_AND_muonDecayedInSample_gen; Bool_t goodEvent_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_F_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_B_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_U_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_D_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_L_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_R_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_F_det; Bool_t goodEvent_F_det;
Bool_t goodEvent_B_det; Bool_t goodEvent_B_det;
Bool_t goodEvent_U_det; Bool_t goodEvent_U_det;
@ -284,6 +300,25 @@ public :
Bool_t promptPeakL; Bool_t promptPeakL;
Bool_t promptPeakR; Bool_t promptPeakR;
Bool_t doubleHit; Bool_t doubleHit;
// pulsed:
// goodEvent_det = (discrim triggered, however happened)
// goodEvent_gen = (1st or only over-threshold event, even if killed by dead time)
// doubleHitEvent_gen = (2nd onwards of a double hit sequence, whether detected or dead)
Bool_t doubleHitEvent_gen;
// doubleHit = (2nd and subsequent hits, counted)
// to subdivide double counting among banks (2nd and subsequent)
Bool_t goodEvent_F_det_AND_doubleHit;
Bool_t goodEvent_B_det_AND_doubleHit;
Bool_t goodEvent_U_det_AND_doubleHit;
Bool_t goodEvent_D_det_AND_doubleHit;
Bool_t goodEvent_L_det_AND_doubleHit;
Bool_t goodEvent_R_det_AND_doubleHit;
// singleHitEvent_gen = (single hits only)
Bool_t singleHitEvent_gen;
// stackedEvent_gen = (multiple small hits adding to threshold, muon details refer to last one)
Bool_t stackedEvent_gen;
// pileupEvent = (all events killed by dead time)
// ..._AND_pileupEvent = subdivided by bank (note the _det isn't appropriate here)
musrAnalysis(TTree *tree=0); musrAnalysis(TTree *tree=0);
virtual ~musrAnalysis(); virtual ~musrAnalysis();
@ -298,6 +333,7 @@ public :
virtual void CreateHistograms(); virtual void CreateHistograms();
virtual void AnalyseEvent(Long64_t iiiEntry); virtual void AnalyseEvent(Long64_t iiiEntry);
virtual void FillHistograms(Int_t iiiEntry); virtual void FillHistograms(Int_t iiiEntry);
virtual void FillHistogramsPulsed(Int_t iiiEntry1, Int_t iiiEntry2);
virtual void SaveHistograms(char* runChar,char* v1190FileName); virtual void SaveHistograms(char* runChar,char* v1190FileName);
virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit); virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit);
// virtual void RewindAllTimeInfo(Double_t timeToRewind); // virtual void RewindAllTimeInfo(Double_t timeToRewind);
@ -319,11 +355,13 @@ public :
TH1D* hGeantParameters; TH1D* hGeantParameters;
TH1D* hInfo; TH1D* hInfo;
typedef std::multimap<double,int> doubleHitSorterType;
// typedef std::map<int,int> debugEventMapType; // typedef std::map<int,int> debugEventMapType;
// debugEventMapType debugEventMap; // debugEventMapType debugEventMap;
// Bool_t bool_debugingRequired; // Bool_t bool_debugingRequired;
static const Double_t pi=3.14159265358979324; //JSL static const Double_t pi=3.14159265358979324;
static musrWriteDump* myWriteDump; static musrWriteDump* myWriteDump;
private: private:
@ -335,9 +373,10 @@ public :
Int_t mdelay; Int_t mdelay;
Int_t pdelay; Int_t pdelay;
Int_t mcoincwin; Int_t mcoincwin;
Int_t pcoincwin; Int_t pcoincwin; // pulsed: for dead time or add-up overlap. Model as exponential decay of this time constant?
Int_t vcoincwin; Int_t vcoincwin;
Double_t muonRateFactor; Double_t muonRateFactor; // also mean rate for pulsed mode
Double_t muonPulseWidthFactor; // pulsed
Double_t dataWindowMin; Double_t dataWindowMin;
Double_t dataWindowMax; Double_t dataWindowMax;
Double_t pileupWindowMin; Double_t pileupWindowMin;
@ -346,16 +385,22 @@ public :
Double_t promptPeakWindowMax; Double_t promptPeakWindowMax;
Int_t overallBinDelay; Int_t overallBinDelay;
Bool_t boolInfinitelyLowMuonRate; Bool_t boolInfinitelyLowMuonRate;
Double_t frameInterval; // JSL: pulsed: to divide muons into frames
Double_t detRecoveryTime; // JSL: pulsed: for dead time effects
Bool_t doPartialFrameAtEnd; // process the left over muons if the beam goes off mid frame having done some complete frames
// Bool_t processInBulk; // do all the muons in one go rather than re-filtering at end of run
Double_t commonThreshold; // set all thresholds to this value regardless of individual settings
char musrMode; // D = time diferential; I = time integral; char musrMode; // D = time diferential; I = time integral; P = pulsed time differential
Double_t safeTimeWindow; Double_t safeTimeWindow;
Double_t currentTime; Double_t currentTime;
// Double_t nextEventTime; // Double_t nextEventTime;
Double_t nextUnfilledEventTime; Double_t nextUnfilledEventTime;
Long64_t numberOfRewinds; Long64_t numberOfRewinds; // JSL: number of frames for Pulsed Mode
Long64_t numberOfGoodMuons; Long64_t numberOfGoodMuons;
// Int_t currentEventID; // Int_t currentEventID;
Long64_t lastPreprocessedEntry; Long64_t lastPreprocessedEntry;
Long64_t firstPreprocessedEntry; // JSL: first in this frame
musrCounter* mCounter; musrCounter* mCounter;
typedef std::map<int,musrCounter*> counterMapType; typedef std::map<int,musrCounter*> counterMapType;
counterMapType pCounterMap; counterMapType pCounterMap;
@ -367,9 +412,9 @@ public :
Bool_t bool_muDecayTimeTransformation; Bool_t bool_muDecayTimeTransformation;
Double_t muDecayTime_Transformation_min, muDecayTime_Transformation_max, muDecayTime_t_min, muDecayTime_t_max; Double_t muDecayTime_Transformation_min, muDecayTime_Transformation_max, muDecayTime_t_min, muDecayTime_t_max;
static const Double_t microsecond=1.; // static const Double_t microsecond=1.;
static const Double_t nanosecond=0.001; // static const Double_t nanosecond=0.001;
static const Double_t picosecond=0.000001; // static const Double_t picosecond=0.000001;
// static const Double_t rewindTime=1000000.; // static const Double_t rewindTime=1000000.;
// static const Double_t rewindTime=1000.; // static const Double_t rewindTime=1000.;
// static const Long64_t rewindTimeBins=1000000000000000; // Max Long64_t can be +9,223,372,036,854,775,807 // static const Long64_t rewindTimeBins=1000000000000000; // Max Long64_t can be +9,223,372,036,854,775,807
@ -383,6 +428,7 @@ public :
public: public:
static const Int_t nrConditions = 31; static const Int_t nrConditions = 31;
Bool_t condition[nrConditions]; Bool_t condition[nrConditions];
Long64_t conditionCounter[nrConditions];
static Long64_t rewindTimeBins; static Long64_t rewindTimeBins;
// static Int_t clock_channelID; // static Int_t clock_channelID;
// static Long64_t clock_interval; // static Long64_t clock_interval;
@ -488,7 +534,10 @@ musrAnalysis::musrAnalysis(TTree *tree)
variableMap["posIniMomZ"]=&posIniMomZ; variableMap["posIniMomZ"]=&posIniMomZ;
// variableMap["nFieldNomVal"]=&nFieldNomVal_double; // variableMap["nFieldNomVal"]=&nFieldNomVal_double;
// variableMap["fieldNomVal0"]=...; //[nFieldNomVal] // variableMap["fieldNomVal0"]=...; //[nFieldNomVal]
variableMap["fieldNomVal0"]=&fieldNomVal[0]; // JSL: allow nominal field especially for background muons stopping elsewhere in map
variableMap["fieldNomVal1"]=&fieldNomVal[1]; // a second field if defined. Could add more?
variableMap["det_n"]=&det_n_double; variableMap["det_n"]=&det_n_double;
variableMap["det_multifirst"]=&det_multifirst_double;
// //
variableMap["muDecayPosR"]=&muDecayPosR; variableMap["muDecayPosR"]=&muDecayPosR;
variableMap["wght"]=&wght; variableMap["wght"]=&wght;
@ -517,6 +566,7 @@ musrAnalysis::musrAnalysis(TTree *tree)
variableMap["pos_Phi_MINUS_muDecayPol_Phi360"]=&pos_Phi_MINUS_muDecayPol_Phi360; variableMap["pos_Phi_MINUS_muDecayPol_Phi360"]=&pos_Phi_MINUS_muDecayPol_Phi360;
variableMap["pos_detID"]=&pos_detID; variableMap["pos_detID"]=&pos_detID;
variableMap["pos_detID_doubleHit"]=&pos_detID_doubleHit; variableMap["pos_detID_doubleHit"]=&pos_detID_doubleHit;
variableMap["pos_doubleHit_dPhi"]=&pos_doubleHit_dPhi;
// variableMap["det_time0"]=&det_time0; // variableMap["det_time0"]=&det_time0;
// variableMap["gen_time0"]=&gen_time0; // variableMap["gen_time0"]=&gen_time0;
// variableMap["det_time1"]=&det_time1; // variableMap["det_time1"]=&det_time1;
@ -525,6 +575,7 @@ musrAnalysis::musrAnalysis(TTree *tree)
variableMap["gen_time10"]=&gen_time10; variableMap["gen_time10"]=&gen_time10;
variableMap["det_time10_MINUS_gen_time10"]=&det_time10_MINUS_gen_time10; variableMap["det_time10_MINUS_gen_time10"]=&det_time10_MINUS_gen_time10;
variableMap["det_time1_MINUS_muDecayTime"]=&det_time1_MINUS_muDecayTime; variableMap["det_time1_MINUS_muDecayTime"]=&det_time1_MINUS_muDecayTime;
variableMap["multiHitInterval"]=&det_multi_interval;
variableMap["detP_x"]=&detP_x; variableMap["detP_x"]=&detP_x;
variableMap["detP_y"]=&detP_x; variableMap["detP_y"]=&detP_x;
variableMap["detP_z"]=&detP_x; variableMap["detP_z"]=&detP_x;

122
musrSimAna/musrCounter.cxx
View File

@ -16,13 +16,15 @@ Bool_t musrCounter::bool_WriteDataToDumpFile = false;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT) { musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT, double DEADTIME, Bool_t KEEP_BELOW_THRESH) {
// pointerToAnalysis=this; // pointerToAnalysis=this;
counterNr = CHANNEL_NR; counterNr = CHANNEL_NR;
strcpy(counterName,CHANNEL_NAME); strcpy(counterName,CHANNEL_NAME);
counterType = CHANNEL_TYPE; counterType = CHANNEL_TYPE;
couterEThreshold = E_THRESH; couterEThreshold = E_THRESH;
counterTimeShift = (Long64_t) TIME_SHIFT; counterTimeShift = (Long64_t) TIME_SHIFT;
keepBelowThresholdHits = KEEP_BELOW_THRESH;
discriminatorRecoveryTime = DEADTIME;
std::cout<<"musrCounter::musrCounter: Creating counter "<<counterNr<<" "<<counterName<<" "<<counterType<<" "<<couterEThreshold<<" "<<counterTimeShift<<std::endl; std::cout<<"musrCounter::musrCounter: Creating counter "<<counterNr<<" "<<counterName<<" "<<counterType<<" "<<couterEThreshold<<" "<<counterTimeShift<<std::endl;
strcpy(TDC_histoName,"Unset"); strcpy(TDC_histoName,"Unset");
TDC_t0=0; TDC_t0=0;
@ -73,13 +75,43 @@ void musrCounter::DrawTDChistogram() {
} }
//================================================================ //================================================================
void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID, Int_t eventNum){ void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID, Int_t eventNum, Int_t& multiCtr){
static hitInfo *hFlagged; // common to all histograms
//cDEL std::cout<<"FillHitInCounter: timeBin="<<timeBin<<" timeBin2="<<timeBin2<<" counterTimeShift="<< counterTimeShift<<std::endl; //cDEL std::cout<<"FillHitInCounter: timeBin="<<timeBin<<" timeBin2="<<timeBin2<<" counterTimeShift="<< counterTimeShift<<std::endl;
//cDEL std::cout<<" FillHitInCounter I: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl; //cDEL std::cout<<" FillHitInCounter I: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl;
// std::cout<<"musrCounter::FillHitInCounter:"<<counterNr<<std::endl; // std::cout<<"musrCounter::FillHitInCounter:"<<counterNr<<std::endl;
if (edep>=couterEThreshold) { // JSL: allow keeping of small hits (but label them)
if ((edep>=couterEThreshold) || ((counterType=='P') && keepBelowThresholdHits)) {
hitInfo* hInfo = new hitInfo(kEntry,eveID,iDet,detectorID,edep,timeBin2-counterTimeShift); // multiple count labelling (JSL)
// counter provided by caller (per event across all counters), reset for new muon
// first big hit numbered 0
// first hit relabelled 1 when second significant hit comes in
// second and subsequent big hits labelled 2,3,etc
// small ones labelled -1 regardless of sequence
// Only for P counters; V,K,M all labelled 0.
Int_t multiCtrTmp;
if(counterType=='P') {
if(edep>=couterEThreshold) {
if((multiCtr==1) && hFlagged) {
// a second good hit, relabel the first hit
// std::cout << "flagged double count in event "<<kEntry<<" hist "<<detectorID<<" time "<<timeBin<<std::endl;
hFlagged->multiHitCtr = 1;
// hFlagged=NULL;
multiCtr++;
}
multiCtrTmp=multiCtr;
multiCtr++;
} else {
multiCtrTmp=-1;
}
} else {
multiCtrTmp=0; // M, V or K counter
}
hitInfo* hInfo = new hitInfo(kEntry,eveID,iDet,detectorID,edep,timeBin-counterTimeShift,timeBin2-counterTimeShift, multiCtrTmp, NULL,0);
if((multiCtrTmp==0) && (counterType=='P')) { hFlagged=hInfo; } // save pointer to the first hit
if(multiCtrTmp>0) {
hInfo->firstMulti = hFlagged;
}
//cDEL std::cout<<detectorID<<" FillHitInCounter II: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl; //cDEL std::cout<<detectorID<<" FillHitInCounter II: timeBin-counterTimeShift="<<timeBin-counterTimeShift<<" timeBin2-counterTimeShift="<<timeBin2-counterTimeShift<<std::endl;
hitMap.insert( std::pair<Long64_t,hitInfo*>(timeBin-counterTimeShift,hInfo) ); hitMap.insert( std::pair<Long64_t,hitInfo*>(timeBin-counterTimeShift,hInfo) );
if (bool_WriteDataToDumpFile) { // write data into an output file if (bool_WriteDataToDumpFile) { // write data into an output file
@ -93,16 +125,20 @@ void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t tim
//================================================================ //================================================================
void musrCounter::RemoveHitsInCounter(Long64_t timeBinLimit) { void musrCounter::RemoveHitsInCounter(Long64_t timeBinLimit) {
// Remove the obsolete hits (i.e. hits that happaned well before t0) from the Counter class // Remove the obsolete hits (i.e. hits that happaned well before t0) from the Counter class
hitMap_TYPE::iterator it2; // JSL
if (hitMap.empty()) return; if (hitMap.empty()) return;
// myPrintThisCounter(); // myPrintThisCounter();
// if (counterNr==1) {std::cout<<"ooooo1 timeBinLimit="<<timeBinLimit<<std::endl; myPrintThisCounter();} // if (counterNr==1) {std::cout<<"ooooo1 timeBinLimit="<<timeBinLimit<<std::endl; myPrintThisCounter();}
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end();) { // JSL: move ++it into body
// std::cout<<" musrCounter::RemoveHitsInCounter: counterNr="<<counterNr<<" timeBinLimit="<<timeBinLimit<<" maxCoincidenceTimeWindow="<<maxCoincidenceTimeWindow<<" counterTimeShift="<<counterTimeShift<<std::endl; // std::cout<<" musrCounter::RemoveHitsInCounter: counterNr="<<counterNr<<" timeBinLimit="<<timeBinLimit<<" maxCoincidenceTimeWindow="<<maxCoincidenceTimeWindow<<" counterTimeShift="<<counterTimeShift<<std::endl;
if ((it->first)>(timeBinLimit+maxCoincidenceTimeWindow-counterTimeShift)) return; //note that maxCoincidenceTimeWindow is usually negative number if ((it->first)>(timeBinLimit+maxCoincidenceTimeWindow-counterTimeShift)) return; //note that maxCoincidenceTimeWindow is usually negative number
else { else {
// std::cout<<" Deleting hit from counter "<<counterNr<<", time bin = "<<(it->first)<<std::endl; // std::cout<<" Deleting hit from counter "<<counterNr<<", time bin = "<<(it->first)<<std::endl;
it2=it;
it2++;
delete (it->second); delete (it->second);
hitMap.erase(it); hitMap.erase(it);
it=it2;
} }
} }
// if (counterNr==1) {std::cout<<"ooooo2"<<std::endl; myPrintThisCounter();} // if (counterNr==1) {std::cout<<"ooooo2"<<std::endl; myPrintThisCounter();}
@ -344,6 +380,80 @@ Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_
return positronQuality; return positronQuality;
} }
//================================================================
Int_t musrCounter::GetNextGoodPositronPulsed(Int_t entWanted1, Int_t entWanted2, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Int_t& idetFirstMulti, hitInfo*& hitStruct) {
// pulsed. Sort through one detector
// return positron quality: (>= -1 are "counted" events, 0,1,-2,-3 are "ideal")
// -1000: no more positron candidates
// 0: a normal "good" positron, the only above threshold hit for that muon
// 1..n: "good" positrons, but double count events. 1 - 1st such hit 2=2nd, etc
// -1: "stacked" event, this positron plus remaining "charge" adds to threshold (double hit status of such events not recorded)
// -2..-999: "dead" event, good positron but will not be counted. -2=only hit -3=1st of double -4=2nd of double, etc
Int_t positronQuality=-1000;
Int_t idMin=+1000000000;
Int_t idMax=-1000000000;
//
// entWanted: must be in range
// timeBinMin, timeBinMax: range to search
// timeBinOfNextGoodHit: starting point on entry, filled with time. If < timeBinMin, reset discriminator to beginning and run through pre-time-window events
// timeBinOfNextGoodHit_phaseShifted: adjusted (out only)
// kEntry: entry number for this event
// idet: detector number
// idetID:
// idetEdep: energy deposited by this positron event (may be < threshold!)
//
// may no longer need to unpackage all the vars?
if(timeBinOfNextGoodHit < timeBinMin || timeBinOfNextGoodHit < discriminatorLastTime) {
// rewind discriminator
discriminatorLastVolts = 0.0;
discriminatorLastTime = timeBinOfNextGoodHit;
discriminatorIterator = hitMap.begin();
}
// process events until a good hit found, or end
while((discriminatorIterator != hitMap.end()) && (discriminatorIterator->first <= timeBinMax)) {
if(discriminatorIterator->second->eventEntry > idMax) { idMax=discriminatorIterator->second->eventEntry; }
if(discriminatorIterator->second->eventEntry < idMin) { idMin=discriminatorIterator->second->eventEntry; }
Long_t dT=((discriminatorIterator->first)-discriminatorLastTime);
double discriminatorDropVolts=discriminatorLastVolts * exp(-dT/discriminatorRecoveryTime);
discriminatorLastVolts=discriminatorDropVolts + discriminatorIterator->second->det_edep;
discriminatorLastTime=discriminatorIterator->first;
// discriminate it!
Bool_t discriminatorTriggered = (discriminatorDropVolts < couterEThreshold && discriminatorLastVolts > couterEThreshold);
// have we found the wanted event?
if ((((discriminatorIterator->second->eventEntry >= entWanted1) && (discriminatorIterator->second->eventEntry <= entWanted2)) && discriminatorIterator->first > timeBinMin) && (discriminatorTriggered || (discriminatorIterator->second->det_edep > couterEThreshold))) {
//yes
if(discriminatorTriggered) {
positronQuality=discriminatorIterator->second->multiHitCtr; // already checked to see if it was above threshold on its own and/or double
} else { // a missed good count, record these too
positronQuality= -2 - discriminatorIterator->second->multiHitCtr;
}
kEntry = (discriminatorIterator->second)->eventEntry;
idet = (discriminatorIterator->second)->det_i;
idetID = (discriminatorIterator->second)->det_id;
idetEdep = (discriminatorIterator->second)->det_edep;
timeBinOfNextGoodHit = discriminatorIterator->first;
timeBinOfNextGoodHit_phaseShifted = (discriminatorIterator->second) -> timeBin2;
if(((discriminatorIterator->second)->multiHitCtr > 1) && (discriminatorIterator->second)->firstMulti) {
idetFirstMulti = (discriminatorIterator->second)->firstMulti->det_i;
} else {
idetFirstMulti = -1000; // not a second or subsequent multi hit, or not recorded for some reason
}
hitStruct = discriminatorIterator->second;
discriminatorIterator->second->posQual = positronQuality;
discriminatorIterator++;
break;
} else {
discriminatorIterator++;
hitStruct = NULL;
}
}
if(positronQuality==-1000 && ((idMax>idMin) && timeBinOfNextGoodHit<0)) {
// std::cout << "Didn't find entry " << entWanted << " but only range " << idMin << " to " << idMax << " from " << timeBinOfNextGoodHit << std::endl;
}
return positronQuality;
}
//================================================================ //================================================================
void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max) { void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max) {
// std::cout<<"QQQQQQQQQQQQQQQQQQQQQ koincidenceCounterMap.size()="<<koincidenceCounterMap.size()<<std::endl; // std::cout<<"QQQQQQQQQQQQQQQQQQQQQ koincidenceCounterMap.size()="<<koincidenceCounterMap.size()<<std::endl;

17
musrSimAna/musrCounter.hh
View File

@ -11,7 +11,7 @@
class hitInfo { class hitInfo {
public: public:
hitInfo(Int_t kEntry, Int_t evID, Int_t deI, Int_t detectorID, Double_t deEDEP, Long64_t timeBIN2) {eventEntry=kEntry; eventIDnumber=evID; det_i=deI; det_id = detectorID; det_edep=deEDEP; timeBin2=timeBIN2;} hitInfo(Int_t kEntry, Int_t evID, Int_t deI, Int_t detectorID, Double_t deEDEP, Long64_t timeBIN1, Long64_t timeBIN2, Int_t multiCtr, hitInfo* firstMulti, Int_t posQuality) {eventEntry=kEntry; eventIDnumber=evID; det_i=deI; det_id = detectorID; det_edep=deEDEP; timeBin1=timeBIN1; timeBin2=timeBIN2; multiHitCtr=multiCtr; firstMulti=NULL;posQual=posQuality;}
~hitInfo() {} ~hitInfo() {}
void RewindTimeBin2(Long64_t timeBinsToRewind) {timeBin2-=timeBinsToRewind;} void RewindTimeBin2(Long64_t timeBinsToRewind) {timeBin2-=timeBinsToRewind;}
Int_t eventEntry; Int_t eventEntry;
@ -19,7 +19,11 @@ public:
Int_t det_i; Int_t det_i;
Int_t det_id; Int_t det_id;
Double_t det_edep; Double_t det_edep;
Long64_t timeBin1; // repeated here in case this is not stored in a Map with First=timeBin1
Long64_t timeBin2; Long64_t timeBin2;
Int_t multiHitCtr;
hitInfo* firstMulti;
Int_t posQual;
// extern double GlobalKamil; // extern double GlobalKamil;
// typedef std::map<int,int> debugEventMapType; // typedef std::map<int,int> debugEventMapType;
@ -29,7 +33,7 @@ public:
class musrCounter { class musrCounter {
public: public:
musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT); musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT, double DEADTIME, Bool_t KEEP_BELOW_THRESH);
~musrCounter(); ~musrCounter();
int GetCounterNr() {return counterNr;} int GetCounterNr() {return counterNr;}
char GetCounterType() {return counterType;} char GetCounterType() {return counterType;}
@ -57,7 +61,7 @@ class musrCounter {
void SetTDChistogram(char hName[200],int t0,int t1,int t2,int hNr,char hNameAdd[200]); void SetTDChistogram(char hName[200],int t0,int t1,int t2,int hNr,char hNameAdd[200]);
void FillTDChistogram(Double_t variable, Double_t vaha); void FillTDChistogram(Double_t variable, Double_t vaha);
void DrawTDChistogram(); void DrawTDChistogram();
void FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID, Int_t eventNum); void FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID, Int_t eventNum, Int_t& multiCtr);
void RemoveHitsInCounter(Long64_t timeBinLimit); void RemoveHitsInCounter(Long64_t timeBinLimit);
// void RewindHitsInCounter(Long64_t timeBinsToRewind); // void RewindHitsInCounter(Long64_t timeBinsToRewind);
void RewindHitsInCounter(); void RewindHitsInCounter();
@ -65,6 +69,7 @@ class musrCounter {
Bool_t GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep); Bool_t GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep);
Bool_t CheckForPileupMuons(Long64_t timeBin0); Bool_t CheckForPileupMuons(Long64_t timeBin0);
Int_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep); Int_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep);
Int_t GetNextGoodPositronPulsed(Int_t entWanted1, Int_t entWanted2, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Int_t& idetFirstMulti, hitInfo*& hitStruct);
void myPrintThisCounter(Int_t evtID, Int_t detail=2); void myPrintThisCounter(Int_t evtID, Int_t detail=2);
// void CheckClockInfo(Long64_t timeBin); // void CheckClockInfo(Long64_t timeBin);
Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;} Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;}
@ -79,6 +84,7 @@ class musrCounter {
char counterName[200]; char counterName[200];
char counterType; char counterType;
double couterEThreshold; double couterEThreshold;
Bool_t keepBelowThresholdHits; // JSL: for pulsed mode overlap
Long64_t counterTimeShift; Long64_t counterTimeShift;
typedef std::map<int,musrCounter*> counterMapType; typedef std::map<int,musrCounter*> counterMapType;
counterMapType koincidenceCounterMap; counterMapType koincidenceCounterMap;
@ -94,6 +100,11 @@ class musrCounter {
typedef std::map<Long64_t,hitInfo*> hitMap_TYPE; // Long64_t = timeBin, hitInfo = eventID and det_i typedef std::map<Long64_t,hitInfo*> hitMap_TYPE; // Long64_t = timeBin, hitInfo = eventID and det_i
hitMap_TYPE hitMap; hitMap_TYPE hitMap;
// std::list<Double_t> timeOfHitsList; // std::list<Double_t> timeOfHitsList;
// for pulsed double hit processing
double discriminatorLastVolts;
Long64_t discriminatorLastTime;
hitMap_TYPE::iterator discriminatorIterator;
Double_t discriminatorRecoveryTime; // in TDC bins!
Int_t doubleHitN; Int_t doubleHitN;
Long64_t numberOfMuonCandidates; Long64_t numberOfMuonCandidates;

4
musrSimAna/musrTH.cxx
View File

@ -49,7 +49,7 @@ void musrTH::FillTH1D(Double_t vaha, Bool_t* cond){
for (Int_t i=0; i<musrAnalysis::nrConditions; i++) { for (Int_t i=0; i<musrAnalysis::nrConditions; i++) {
if (bool_rotating_reference_frame) { if (bool_rotating_reference_frame) {
// Double_t var = *variableToBeFilled_X; // Double_t var = *variableToBeFilled_X;
Double_t waha = vaha*exp((*variableToBeFilled_X)/2.19703)*cos(2*musrAnalysis::pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase); Double_t waha = vaha*exp((*variableToBeFilled_X)/2.19703)*cos(2*pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase);
// std::cout<<"rot_ref_frequency="<<rot_ref_frequency<<std::endl; // std::cout<<"rot_ref_frequency="<<rot_ref_frequency<<std::endl;
if (cond[i]) histArray1D[i]->Fill(*variableToBeFilled_X,waha); if (cond[i]) histArray1D[i]->Fill(*variableToBeFilled_X,waha);
} }
@ -225,7 +225,7 @@ void musrTH::ListHistograms() {
//============================================================================================== //==============================================================================================
void musrTH::FitHistogramsIfRequired(Double_t omega) { void musrTH::FitHistogramsIfRequired(Double_t omega) {
if (funct==NULL) return; if (funct==NULL) return;
if (bool_rotating_reference_frame) omega = fabs(omega) - 2*musrAnalysis::pi*fabs(rot_ref_frequency); if (bool_rotating_reference_frame) omega = fabs(omega) - 2*pi*fabs(rot_ref_frequency);
std::cout<<"============================================================================================================"<<std::endl; std::cout<<"============================================================================================================"<<std::endl;
std::cout<<"Fitting \""<<funct->GetName()<<"\", funct_xMin="<<funct_xMin<<" funct_xMax="<<funct_xMax<<" omega="<<omega<<std::endl; std::cout<<"Fitting \""<<funct->GetName()<<"\", funct_xMin="<<funct_xMin<<" funct_xMax="<<funct_xMax<<" omega="<<omega<<std::endl;
if (strcmp(funct->GetName(),"funct1")==0) {funct->FixParameter(0,omega);} if (strcmp(funct->GetName(),"funct1")==0) {funct->FixParameter(0,omega);}

20
src/meyer.cc
View File

@ -60,6 +60,10 @@
#include <iostream> #include <iostream>
#include <stdlib.h> #include <stdlib.h>
#include <ios> #include <ios>
// James Lord:
// already defines its own constant Pi so no need to look for M_PI anywhere
// #define _USE_MATH_DEFINES 1
// #include <math.h>
using namespace std; using namespace std;
meyer::meyer() meyer::meyer()
@ -258,7 +262,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2,
} }
// std::cout<< "M4: "<<Meyer_faktor4<<std::endl; // std::cout<< "M4: "<<Meyer_faktor4<<std::endl;
// std::cout<< "Ekin: "<<Ekin <<std::endl; // std::cout<< "Ekin: "<<Ekin <<std::endl;
thetaMax = thetaSchlangeMax / Meyer_faktor4 / Ekin/M_PI*180; thetaMax = thetaSchlangeMax / Meyer_faktor4 / Ekin/Pi*180;
if (thetaMax>50.) if (thetaMax>50.)
{ {
thetaStep = .5; thetaStep = .5;
@ -292,7 +296,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2,
// ! Berechne aus theta das 'reduzierte' thetaSchlange (dabei gleich // ! Berechne aus theta das 'reduzierte' thetaSchlange (dabei gleich
// ! noch von degree bei theta in Radiant bei thetaSchlange umrechnen): // ! noch von degree bei theta in Radiant bei thetaSchlange umrechnen):
// //
thetaSchlange = Meyer_faktor4 * Ekin * theta *M_PI/180; thetaSchlange = Meyer_faktor4 * Ekin * theta *Pi/180;
// ! Auslesen der Tabellenwerte fuer die f-Funktionen: // ! Auslesen der Tabellenwerte fuer die f-Funktionen:
@ -305,7 +309,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2,
} }
// ! Berechnen der Streuintensitaet: // ! Berechnen der Streuintensitaet:
F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /M_PI * (f1 - Meyer_faktor3*f2);// TAO, Anselm was: Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2); F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /Pi * (f1 - Meyer_faktor3*f2);// TAO, Anselm was: Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2);
nBin = nBin + 1; nBin = nBin + 1;
if (nBin>nBinMax) if (nBin>nBinMax)
@ -317,7 +321,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2,
value[nBin] = sin(theta)*F; value[nBin] = sin(theta)*F;
fValues[nBin+1] = F; // ! fuer Testzwecke fValues[nBin+1] = F; // ! fuer Testzwecke
fValuesFolded[nBin+1] = sin(theta/180*M_PI)*F;// ! fuer Testzwecke fValuesFolded[nBin+1] = sin(theta/180*Pi)*F;// ! fuer Testzwecke
}// end of do loop }// end of do loop
@ -348,7 +352,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2,
//! Berechne die Werte fuer theta=0: //! Berechne die Werte fuer theta=0:
F_Functions_Meyer(tau,0.,&f1,&f2); F_Functions_Meyer(tau,0.,&f1,&f2);
F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /M_PI * (f1 - Meyer_faktor3*f2);// TAO, Anselm was: Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2); F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /Pi * (f1 - Meyer_faktor3*f2);// TAO, Anselm was: Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2);
fValues[1] = F; fValues[1] = F;
fValuesFolded[1] = 0.; fValuesFolded[1] = 0.;
@ -463,7 +467,7 @@ void meyer:: Get_F_Function(double tau,double theta, double Ekin, double Z1, dou
nBin = 0; nBin = 0;
thetaSchlange = Meyer_faktor4 * Ekin * theta *M_PI/180; thetaSchlange = Meyer_faktor4 * Ekin * theta *Pi/180;
F_Functions_Meyer(tau,thetaSchlange,&f1,&f2); F_Functions_Meyer(tau,thetaSchlange,&f1,&f2);
@ -472,7 +476,7 @@ void meyer:: Get_F_Function(double tau,double theta, double Ekin, double Z1, dou
goto bigtheta; goto bigtheta;
} }
*F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /M_PI * (f1 - Meyer_faktor3*f2); *F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /Pi * (f1 - Meyer_faktor3*f2);
bigtheta:for( i = 1;i<= nBin; i++) bigtheta:for( i = 1;i<= nBin; i++)
{ {
@ -486,7 +490,7 @@ void meyer:: Get_F_Function(double tau,double theta, double Ekin, double Z1, dou
//! Berechne die Werte fuer theta=0: //! Berechne die Werte fuer theta=0:
double norm; double norm;
F_Functions_Meyer(tau,0.,&f1,&f2); F_Functions_Meyer(tau,0.,&f1,&f2);
norm = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /M_PI * (f1 - Meyer_faktor3*f2); norm = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /Pi * (f1 - Meyer_faktor3*f2);
*F=*F/norm; *F=*F/norm;
} }

70
src/musrDetectorConstruction.cc
View File

@ -44,6 +44,7 @@
#include "G4Box.hh" #include "G4Box.hh"
#include "G4Cons.hh" #include "G4Cons.hh"
#include "G4Polycone.hh" #include "G4Polycone.hh"
#include "G4Torus.hh"
#include "G4LogicalVolume.hh" #include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh" #include "G4PVPlacement.hh"
#include "G4Tubs.hh" #include "G4Tubs.hh"
@ -329,6 +330,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
solidName+=name; solidName+=name;
solid = new G4Para(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg,x6*deg); solid = new G4Para(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg,x6*deg);
} }
else if (strcmp(tmpString2,"torus")==0){ // added JSL. Torus piece, Rmin (bore), Rmax (outer), Rtorus (swept), pSPhi (start angle), pDPhi (swept angle)
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
solid = new G4Torus(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
}
else if (strcmp(tmpString2,"cylpart")==0){ // Volume introduced by Pavel Bakule on 12 May 2009 else if (strcmp(tmpString2,"cylpart")==0){ // Volume introduced by Pavel Bakule on 12 May 2009
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %s %lf %lf %lf %s %s", sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix); name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix);
@ -665,7 +673,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID); sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID);
solidName+=name; solidName+=name;
G4Box* solidGPDcutOut = new G4Box ("solidGPDcutOut",x1*mm,x1*mm,(x3+0.01)*mm); G4Box* solidGPDcutOut = new G4Box ("solidGPDcutOut",x1*mm,x1*mm,(x3+0.01)*mm);
G4Box* solidGPDmHolder = new G4Box ("solidGPDmHolder",sqrt(2)*x1*mm,x2*mm,x3*mm); G4Box* solidGPDmHolder = new G4Box ("solidGPDmHolder",sqrt((double)2)*x1*mm,x2*mm,x3*mm);
// G4RotationMatrix* rot = new G4RotationMatrix(45*deg,0,0); // G4RotationMatrix* rot = new G4RotationMatrix(45*deg,0,0);
G4RotationMatrix* rot = new G4RotationMatrix(G4ThreeVector(0,0,1),45*deg); G4RotationMatrix* rot = new G4RotationMatrix(G4ThreeVector(0,0,1),45*deg);
// G4RotationMatrix* rot = new G4RotationMatrix(); // G4RotationMatrix* rot = new G4RotationMatrix();
@ -673,6 +681,28 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
solid = new G4SubtractionSolid(solidName,solidGPDmHolder,solidGPDcutOut,rot,trans); solid = new G4SubtractionSolid(solidName,solidGPDmHolder,solidGPDcutOut,rot,trans);
// solid = new G4SubtractionSolid(solidName,solidGPDcutOut,solidGPDmHolder,rot,trans); // solid = new G4SubtractionSolid(solidName,solidGPDcutOut,solidGPDmHolder,rot,trans);
} }
else if (strcmp(tmpString2,"cruciform")==0){
// JSL: Create a solid cruciform (union of 3 mutually perpendicular cylinders). Top port can be different to bottom.
// x1=radius of main (+-z) x2=radius of +-x and -y, x3=radius of +y x4=z extent x5=x,y extent
// could enclose another of these inside, filled with vacuum
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
//G4double roundingErr=0.01*mm; // to avoid some displaying problems of the subtracted volumes
G4Tubs* solidMainBore = new G4Tubs("solidMainBore" ,0.*mm,x1*mm,x4*mm,0*deg,360*deg);
G4Tubs* solidHorizCrossBore = new G4Tubs("solidHorizCrossBore",0.*mm,x2*mm,x5*mm,0*deg,360*deg);
G4Tubs* solidBottomBore = new G4Tubs("solidBottomBore" ,0.*mm,x2*mm,x5/2*mm,0*deg,360*deg);
G4Tubs* solidTopBore = new G4Tubs("solidTopBore" ,0.*mm,x3*mm,x5/2*mm,0*deg,360*deg);
G4RotationMatrix* rotx = new G4RotationMatrix(G4ThreeVector(0.,1.,0.),90.*deg);
G4RotationMatrix* roty = new G4RotationMatrix(G4ThreeVector(1.,0.,0.),90.*deg);
G4ThreeVector nowhere( 0,0,0);
G4ThreeVector UpOffset( 0., x5/2.*mm, 0.);
G4ThreeVector DownOffset( 0., -x5/2.*mm, 0.);
G4UnionSolid* partXZonly = new G4UnionSolid("partXZonly", solidMainBore, solidHorizCrossBore, rotx, nowhere);
G4UnionSolid* partXZBonly = new G4UnionSolid("partXZBonly", partXZonly, solidBottomBore, roty, DownOffset);
solid = new G4UnionSolid(solidName, partXZBonly, solidTopBore, roty, UpOffset);
}
else ReportGeometryProblem(line); else ReportGeometryProblem(line);
G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm); G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm);
@ -1631,6 +1661,41 @@ void musrDetectorConstruction::DefineMaterials()
Air->AddElement(N, fractionmass=0.7); Air->AddElement(N, fractionmass=0.7);
Air->AddElement(O, fractionmass=0.3); Air->AddElement(O, fractionmass=0.3);
G4Material* Air1mbar =
new G4Material("Air1mbar" , density= 1.290e-3*mg/cm3, ncomponents=2);
Air1mbar->AddElement(N, fractionmass=0.7);
Air1mbar->AddElement(O, fractionmass=0.3);
G4Material* AirE1mbar =
new G4Material("AirE1mbar" , density= 1.290e-4*mg/cm3, ncomponents=2);
AirE1mbar->AddElement(N, fractionmass=0.7);
AirE1mbar->AddElement(O, fractionmass=0.3);
G4Material* AirE2mbar =
new G4Material("AirE2mbar" , density= 1.290e-5*mg/cm3, ncomponents=2);
AirE2mbar->AddElement(N, fractionmass=0.7);
AirE2mbar->AddElement(O, fractionmass=0.3);
G4Material* AirE3mbar =
new G4Material("AirE3mbar" , density= 1.290e-6*mg/cm3, ncomponents=2);
AirE3mbar->AddElement(N, fractionmass=0.7);
AirE3mbar->AddElement(O, fractionmass=0.3);
G4Material* AirE4mbar =
new G4Material("AirE4mbar" , density= 1.290e-7*mg/cm3, ncomponents=2);
AirE4mbar->AddElement(N, fractionmass=0.7);
AirE4mbar->AddElement(O, fractionmass=0.3);
G4Material* AirE5mbar =
new G4Material("AirE5mbar" , density= 1.290e-8*mg/cm3, ncomponents=2);
AirE5mbar->AddElement(N, fractionmass=0.7);
AirE5mbar->AddElement(O, fractionmass=0.3);
G4Material* AirE6mbar =
new G4Material("AirE6mbar" , density= 1.290e-9*mg/cm3, ncomponents=2);
AirE6mbar->AddElement(N, fractionmass=0.7);
AirE6mbar->AddElement(O, fractionmass=0.3);
// //
// examples of vacuum // examples of vacuum
// //
@ -1661,6 +1726,9 @@ void musrDetectorConstruction::DefineMaterials()
new G4Material("HeliumGas80mbar",z=2., a=4.002602*g/mole, density= 0.00001410112*g/cm3); new G4Material("HeliumGas80mbar",z=2., a=4.002602*g/mole, density= 0.00001410112*g/cm3);
new G4Material("HeliumGas90mbar",z=2., a=4.002602*g/mole, density= 0.00001586376*g/cm3); new G4Material("HeliumGas90mbar",z=2., a=4.002602*g/mole, density= 0.00001586376*g/cm3);
new G4Material("HeliumGas100mbar",z=2.,a=4.002602*g/mole, density= 0.00001762640*g/cm3); new G4Material("HeliumGas100mbar",z=2.,a=4.002602*g/mole, density= 0.00001762640*g/cm3);
new G4Material("HeliumGasSat4K",z=2., a=4.002602*g/mole, density= 0.016891*g/cm3); // saturated vapour above liquid at 4.2K (JSL)
new G4Material("HeliumGas5mbar4K",z=2.,a=4.002602*g/mole, density= 0.016891*5.0/1013.0*g/cm3); // typical cold exchange gas, 4.2K and 5 mbar
new G4Material("HeliumGas2mbar4K",z=2.,a=4.002602*g/mole, density= 0.016891*2.0/1013.0*g/cm3); // typical cold exchange gas, 4.2K and 5 mbar
if (musrParameters::boolG4OpticalPhotons) { if (musrParameters::boolG4OpticalPhotons) {

1
src/musrErrorMessage.cc
View File

@ -47,6 +47,7 @@ void musrErrorMessage::musrError(SEVERITY severity, G4String message, G4bool sil
actualErrorMessage.mesSeverity = severity; actualErrorMessage.mesSeverity = severity;
actualErrorMessage.nTimes = 1; actualErrorMessage.nTimes = 1;
ErrorMapping[message]=actualErrorMessage; ErrorMapping[message]=actualErrorMessage;
it = ErrorMapping.find(message);
G4cout<<"!!!"<<severityWord[severity]<<"!!! "<<message<<" (First time occurrence)"<<G4endl; G4cout<<"!!!"<<severityWord[severity]<<"!!! "<<message<<" (First time occurrence)"<<G4endl;
nErrors++; nErrors++;
} }

4
src/musrEventAction.cc
View File

@ -120,7 +120,7 @@ void musrEventAction::EndOfEventAction(const G4Event* evt) {
fRunManager->AbortRun(true); fRunManager->AbortRun(true);
} }
if (thisEventNr != 0 and thisEventNr%nHowOftenToPrintEvent == 0) { if ((thisEventNr != 0) && (thisEventNr%nHowOftenToPrintEvent == 0)) {
time_t curr=time(0); time_t curr=time(0);
//char * ctime(const time_t * tp); //char * ctime(const time_t * tp);
G4cout << ">>> Event " << evt->GetEventID() <<". Running already for "<<curr-timeOfRunStart<<" seconds. Present time: "<< ctime(&curr); G4cout << ">>> Event " << evt->GetEventID() <<". Running already for "<<curr-timeOfRunStart<<" seconds. Present time: "<< ctime(&curr);
@ -133,7 +133,7 @@ void musrEventAction::EndOfEventAction(const G4Event* evt) {
for (G4int i=0; i<n_trajectories; i++) for (G4int i=0; i<n_trajectories; i++)
{ G4Trajectory* trj = (G4Trajectory*) { G4Trajectory* trj = (G4Trajectory*)
((*(evt->GetTrajectoryContainer()))[i]); ((*(evt->GetTrajectoryContainer()))[i]);
trj->DrawTrajectory(1000); trj->DrawTrajectory(); //removed argument 1000 (JSL)
} }
} }
} }

16
src/musrPhysicsList.cc
View File

@ -264,6 +264,7 @@ void musrPhysicsList::ConstructProcess()
// For optical photons // For optical photons
#include "G4Scintillation.hh" #include "G4Scintillation.hh"
#include "G4Cerenkov.hh"
// For models // For models
#include "G4ProcessTable.hh" #include "G4ProcessTable.hh"
@ -336,6 +337,13 @@ void musrPhysicsList::ConstructEM()
else if (stringProcessName=="G4OpRayleigh") pManager->AddDiscreteProcess(new G4OpRayleigh); else if (stringProcessName=="G4OpRayleigh") pManager->AddDiscreteProcess(new G4OpRayleigh);
else if (stringProcessName=="G4OpBoundaryProcess") pManager->AddDiscreteProcess(new G4OpBoundaryProcess); else if (stringProcessName=="G4OpBoundaryProcess") pManager->AddDiscreteProcess(new G4OpBoundaryProcess);
else if (stringProcessName=="G4OpWLS") pManager->AddDiscreteProcess(new G4OpWLS); else if (stringProcessName=="G4OpWLS") pManager->AddDiscreteProcess(new G4OpWLS);
else if (stringProcessName=="G4Cerenkov") {
G4Cerenkov* myCerenkov = new G4Cerenkov();
myCerenkov->SetMaxNumPhotonsPerStep(20); // settings from GEANT example N06
myCerenkov->SetMaxBetaChangePerStep(10.0);
myCerenkov->SetTrackSecondariesFirst(false); // as for Scintillation, do main high energy particles first
pManager->AddDiscreteProcess(myCerenkov); // added JSL
}
else { else {
sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess. It can be easily added.", sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess. It can be easily added.",
charProcessName); charProcessName);
@ -373,6 +381,14 @@ void musrPhysicsList::ConstructEM()
if (musrParameters::finiteRiseTimeInScintillator) myScintillation->SetFiniteRiseTime(true); if (musrParameters::finiteRiseTimeInScintillator) myScintillation->SetFiniteRiseTime(true);
pManager->AddProcess(myScintillation,nr1,nr2,nr3); pManager->AddProcess(myScintillation,nr1,nr2,nr3);
} }
else if (stringProcessName=="G4Cerenkov") {
G4Cerenkov* myCerenkov = new G4Cerenkov();
myCerenkov->SetMaxNumPhotonsPerStep(20); // settings from GEANT example N06
myCerenkov->SetMaxBetaChangePerStep(10.0);
myCerenkov->SetTrackSecondariesFirst(false); // as for Scintillation, do main high energy particles first
pManager->AddProcess(myCerenkov,nr1,nr2,nr3); // added JSL
pManager->SetProcessOrdering(myCerenkov,idxPostStep); // from Example N06
}
// cks: musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments // cks: musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments
// in his original "musrPhysicsList.cc about implementing musrMuScatter. // in his original "musrPhysicsList.cc about implementing musrMuScatter.
// else if (stringProcessName=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3); // else if (stringProcessName=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3);

3
src/musrRootOutput.cc
View File

@ -833,8 +833,9 @@ void musrRootOutput::SetPhotDetTime(G4double time) {
nOptPhotDet++; nOptPhotDet++;
} }
else { else {
nOptPhotDet++; // still want to know how many in total (JSL)
char message[200]; char message[200];
sprintf(message,"musrRootOutput.cc::SetPhotDetTime: number of individual photons larger than maxNOptPhotDet (=%d)",maxNOptPhotDet); sprintf(message,"musrRootOutput.cc::SetPhotDetTime: number of individual photons larger than maxNOptPhotDet (=%d)",maxNOptPhotDet);
musrErrorMessage::GetInstance()->musrError(WARNING,message,false); musrErrorMessage::GetInstance()->musrError(WARNING,message,true); // had silent=false and printed all messages (JSL)
} }
} }

21
src/musrScintSD.cc
View File

@ -35,6 +35,7 @@
#include "musrErrorMessage.hh" #include "musrErrorMessage.hh"
#include "musrSteppingAction.hh" #include "musrSteppingAction.hh"
//#include "TCanvas.h" //#include "TCanvas.h"
#include "TH1D.h"
#include "TMath.h" #include "TMath.h"
#include "TF1.h" #include "TF1.h"
#include <vector> #include <vector>
@ -131,6 +132,7 @@ musrScintSD::musrScintSD(G4String name)
APDcellsTimeVariationRequested=false; APDcellsTimeVariationRequested=false;
APDcrossTalkRequested=false; APDcrossTalkRequested=false;
APDcrossTalkProb=0.; APDcrossTalkProb=0.;
// verboseLevel = 9; // JSL
} }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -167,8 +169,10 @@ void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
if (musrParameters::boolG4OpticalPhotons) { if (musrParameters::boolG4OpticalPhotons) {
if (!musrParameters::boolG4OpticalPhotonsUnprocess) { if (!musrParameters::boolG4OpticalPhotonsUnprocess) {
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) { for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting optHitDetectorMap" <<"\n";
delete (it->second); delete (it->second);
} }
if (verboseLevel>1) G4cout<<"VERBOSE 2 : clearing optHitMap" <<"\n";
optHitMap.clear(); optHitMap.clear();
} }
} }
@ -208,6 +212,7 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new musrScintHit" <<"\n";
musrScintHit* newHit = new musrScintHit(); musrScintHit* newHit = new musrScintHit();
// newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName()); // newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
newHit->SetParticleName(particleName); newHit->SetParticleName(particleName);
@ -286,6 +291,7 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
// G4cout<<" detID ="<<detID<<" actualVolume="<<actualVolume<<G4endl; // G4cout<<" detID ="<<detID<<" actualVolume="<<actualVolume<<G4endl;
optHitMapType::iterator iter = optHitMap.find(detID); optHitMapType::iterator iter = optHitMap.find(detID);
if (iter==optHitMap.end()) { // optHitDetectorMapType does not exist ==> create it if (iter==optHitMap.end()) { // optHitDetectorMapType does not exist ==> create it
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new OptHitDetectorMap" <<"\n";
optHitDetectorMapType* optHitDetectorMapTmp = new optHitDetectorMapType; optHitDetectorMapType* optHitDetectorMapTmp = new optHitDetectorMapType;
optHitMap.insert(std::pair<Int_t,optHitDetectorMapType*>(detID,optHitDetectorMapTmp)); optHitMap.insert(std::pair<Int_t,optHitDetectorMapType*>(detID,optHitDetectorMapTmp));
iter = optHitMap.find(detID); iter = optHitMap.find(detID);
@ -380,7 +386,7 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
// in which case times are huge and due to the limited rounding // in which case times are huge and due to the limited rounding
// precision become equal --> map ignores the same "keys", // precision become equal --> map ignores the same "keys",
// multimap does not. // multimap does not.
std::map<G4double,G4int>::iterator it; std::multimap<G4double,G4int>::iterator it;
for (G4int i=0; i<NbHits; i++) { for (G4int i=0; i<NbHits; i++) {
musrScintHit* aHit = (*scintCollection)[i]; musrScintHit* aHit = (*scintCollection)[i];
G4double tmptime=aHit->GetGlobalTime(); G4double tmptime=aHit->GetGlobalTime();
@ -670,6 +676,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
iHistNr++; iHistNr++;
char nameHist[200]; sprintf(nameHist,"OPSAhist_%d_%d_%d",eeeventID,OPSA_detID,iHistNr); char nameHist[200]; sprintf(nameHist,"OPSAhist_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
char nameHistTitle[200]; sprintf(nameHistTitle,"OPSAhist_%d_%d_%d;time (ns);Nr of photons",eeeventID,OPSA_detID,iHistNr); char nameHistTitle[200]; sprintf(nameHistTitle,"OPSAhist_%d_%d_%d;time (ns);Nr of photons",eeeventID,OPSA_detID,iHistNr);
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new TH1D for OPSAhisto" <<"\n";
OPSAhisto = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); OPSAhisto = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
// poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2 // poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2
// poiss->SetParameter(0,1); // poiss->SetParameter(0,1);
@ -677,9 +684,11 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
if (bool_pulseShapeExists) { if (bool_pulseShapeExists) {
sprintf(nameHist,"OPSAshape_%d_%d_%d",eeeventID,OPSA_detID,iHistNr); sprintf(nameHist,"OPSAshape_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
sprintf(nameHistTitle,"OPSAshape_%d_%d_%d;time (ns);Pulse signal",eeeventID,OPSA_detID,iHistNr); sprintf(nameHistTitle,"OPSAshape_%d_%d_%d;time (ns);Pulse signal",eeeventID,OPSA_detID,iHistNr);
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new TH1D for OPSAshape" <<"\n";
OPSAshape = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); OPSAshape = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
sprintf(nameHist,"OPSA_CFD_%d_%d_%d",eeeventID,OPSA_detID,iHistNr); sprintf(nameHist,"OPSA_CFD_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
sprintf(nameHistTitle,"OPSA_CFD_%d_%d_%d;time (ns);CFD signal",eeeventID,OPSA_detID,iHistNr); sprintf(nameHistTitle,"OPSA_CFD_%d_%d_%d;time (ns);CFD signal",eeeventID,OPSA_detID,iHistNr);
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new TH1D for OPSA_CFD" <<"\n";
OPSA_CFD = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); OPSA_CFD = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
} }
} }
@ -695,6 +704,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
else { // create histogram because it does not exist else { // create histogram because it does not exist
char nameHistTitle[200]; sprintf(nameHistTitle,"OPSAhistSUM_%d_%d;time (ns);Nr of photons",OPSA_detID,iHistNrSUM); char nameHistTitle[200]; sprintf(nameHistTitle,"OPSAhistSUM_%d_%d;time (ns);Nr of photons",OPSA_detID,iHistNrSUM);
char nameHistTitle0[200]; sprintf(nameHistTitle0,"OPSAhistSUM0_%d_%d;time (ns);Nr of photons",OPSA_detID,iHistNrSUM); char nameHistTitle0[200]; sprintf(nameHistTitle0,"OPSAhistSUM0_%d_%d;time (ns);Nr of photons",OPSA_detID,iHistNrSUM);
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating new TH1Ds for OPSAhistoSUM,0" <<"\n";
OPSAhistoSUM = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); OPSAhistoSUM = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
OPSAhistoSUM0= new TH1D(nameHist0,nameHistTitle0,OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); OPSAhistoSUM0= new TH1D(nameHist0,nameHistTitle0,OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
mapOfOPSAsumHistograms.insert(std::pair<std::string,TH1D*>(nameHist,OPSAhistoSUM)); mapOfOPSAsumHistograms.insert(std::pair<std::string,TH1D*>(nameHist,OPSAhistoSUM));
@ -1026,13 +1036,16 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
// Delete the histograms from memory if they were created // Delete the histograms from memory if they were created
if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) { if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) {
if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting OPSAhisto" <<"\n";
delete OPSAhisto; delete OPSAhisto;
if (bool_pulseShapeExists) { if (bool_pulseShapeExists) {
if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting OPSAshape and CFD" <<"\n";
delete OPSAshape; delete OPSAshape;
delete OPSA_CFD; delete OPSA_CFD;
} }
} }
if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating new signalInfo" <<"\n";
signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_nPhotPrim,OPSA_timeFirst,OPSA_timeSecond,OPSA_timeThird, signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_nPhotPrim,OPSA_timeFirst,OPSA_timeSecond,OPSA_timeThird,
OPSA_timeA,OPSA_timeB,OPSA_timeC,OPSA_timeD,OPSA_timeMean,OPSA_timeLast, OPSA_timeA,OPSA_timeB,OPSA_timeC,OPSA_timeD,OPSA_timeMean,OPSA_timeLast,
OPSA_CFD_time,OPSA_CFD_ampl,timeCFDarray,OPSA_timeC1); OPSA_CFD_time,OPSA_CFD_ampl,timeCFDarray,OPSA_timeC1);
@ -1070,6 +1083,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
} }
// Now delete all mySignalInfo* from OPSA_signal_Map // Now delete all mySignalInfo* from OPSA_signal_Map
for (OPSA_signal_MapType::const_iterator itOPSA = OPSA_signal_Map.begin(); itOPSA != OPSA_signal_Map.end(); itOPSA++) { for (OPSA_signal_MapType::const_iterator itOPSA = OPSA_signal_Map.begin(); itOPSA != OPSA_signal_Map.end(); itOPSA++) {
if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting itOPSA->second" <<"\n";
delete (itOPSA->second); delete (itOPSA->second);
} }
OPSA_signal_Map.clear(); OPSA_signal_Map.clear();
@ -1162,16 +1176,20 @@ G4int musrScintSD::FindAPDcellID(G4Step* aStep) {
void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time, G4int nTruePhe) { void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time, G4int nTruePhe) {
if ( (musrRootOutput::store_phot_time) && (nTruePhe>0) ) myRootOutput->SetPhotDetTime(time); if ( (musrRootOutput::store_phot_time) && (nTruePhe>0) ) myRootOutput->SetPhotDetTime(time);
if (!APDcellsEffectRequested) { if (!APDcellsEffectRequested) {
if (verboseLevel>1) G4cout<<"VERBOSE 2 : storing photon hit (non APD)" <<"\n";
APDidClass tmpAPDid(0,nTruePhe); APDidClass tmpAPDid(0,nTruePhe);
// optHitDetectorMap->insert(std::pair<G4double,G4int>(time,0)); // optHitDetectorMap->insert(std::pair<G4double,G4int>(time,0));
optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid)); optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid));
} }
else { else {
G4bool APDcellAlreadyFired = false; G4bool APDcellAlreadyFired = false;
if (verboseLevel>1) G4cout<<"VERBOSE 2 : searching for the APD element" <<"\n";
for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) { for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
if ((it2->second).GetAPDcellID()==APDcellID) { // this cell already fired before ==> check times if ((it2->second).GetAPDcellID()==APDcellID) { // this cell already fired before ==> check times
if (verboseLevel>1) G4cout<<"VERBOSE 2 : this cell already fired before" <<"\n";
APDcellAlreadyFired = true; APDcellAlreadyFired = true;
if (time<(it2->first)) { // the new cell fired before the already saved cell ==> replace it if (time<(it2->first)) { // the new cell fired before the already saved cell ==> replace it
if (verboseLevel>1) G4cout<<"VERBOSE 2 : and the new photon is earlier, replacing" <<"\n";
G4int tmpAPDcellNphot = (it2->second).GetAPDcellNphot() + nTruePhe; G4int tmpAPDcellNphot = (it2->second).GetAPDcellNphot() + nTruePhe;
APDidClass tmpAPDid(APDcellID,tmpAPDcellNphot); APDidClass tmpAPDid(APDcellID,tmpAPDcellNphot);
optHitDetectorMap->erase(it2); optHitDetectorMap->erase(it2);
@ -1183,6 +1201,7 @@ void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int AP
} }
if (!APDcellAlreadyFired) { if (!APDcellAlreadyFired) {
// optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID)); // optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID));
if (verboseLevel>1) G4cout<<"VERBOSE 2 : this cell hasn't fired, storing" <<"\n";
APDidClass tmpAPDid(APDcellID,nTruePhe); APDidClass tmpAPDid(APDcellID,nTruePhe);
optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid)); optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid));
} }

11
src/musrTabulatedElementField.cc
View File

@ -617,10 +617,12 @@ void musrTabulatedElementField::addFieldValue2D(const G4double point[4],
local = global2local.TransformPoint(global); local = global2local.TransformPoint(global);
double x, z, z_sign; double x, z, z_sign;
if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)|| // if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)||
(strcmp(fieldTableType,"2D_OperaXY"))||(strcmp(fieldTableType,"2DEf")==0)) { // (strcmp(fieldTableType,"2D_OperaXY"))||(strcmp(fieldTableType,"2DEf")==0)) {
// Field is defined in just positive range of z; i.e. it is expected to be "symmetric" // Field is defined in just positive range of z; i.e. it is expected to be "symmetric"
// and the field for negative z is calculated from the positive z half. // and the field for negative z is calculated from the positive z half.
if (minimumz >= 0) {
// JSL: Use field values for z<0 if they're in the table! Reflect otherwise.
x = sqrt(local.x()*local.x()+local.y()*local.y()); x = sqrt(local.x()*local.x()+local.y()*local.y());
z = fabs(local.z()); z = fabs(local.z());
z_sign = (local.z()>0) ? 1.:-1.; z_sign = (local.z()>0) ? 1.:-1.;
@ -631,8 +633,9 @@ void musrTabulatedElementField::addFieldValue2D(const G4double point[4],
z = local.z(); z = local.z();
z_sign = 1; z_sign = 1;
} }
// Check that the point is within the defined region // Check that the point is within the defined region (JSL: 2-sided check)
if ( x<maximumx && z<maximumz ) { if ( x>=minimumx && x<maximumx &&
z>=minimumz && z<maximumz ) {
// if (evNr>evNrKriz) std::cout<<"bol som tu"<<std::endl; // if (evNr>evNrKriz) std::cout<<"bol som tu"<<std::endl;
// Position of given point within region, normalized to the range // Position of given point within region, normalized to the range