diff --git a/CMakeMacroParseArguments.cmake b/CMakeMacroParseArguments.cmake new file mode 100644 index 0000000..310d722 --- /dev/null +++ b/CMakeMacroParseArguments.cmake @@ -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) diff --git a/CMakeParseArguments.cmake b/CMakeParseArguments.cmake new file mode 100644 index 0000000..cc63604 --- /dev/null +++ b/CMakeParseArguments.cmake @@ -0,0 +1,138 @@ +# CMAKE_PARSE_ARGUMENTS( 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 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 argument contains all keywords for this macro +# which are followed by one value, like e.g. DESTINATION keyword of the +# install() command. +# +# The 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 , and +# a variable composed of the given +# followed by "_" and the name of the respective keyword. +# These variables will then hold the respective value from the argument list. +# For the keywords this will be TRUE or FALSE. +# +# All remaining arguments are collected in a variable +# _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 +# +# 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) diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index 9a52c76..ff30f12 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index e2b1735..8fd8a3f 100644 --- a/doc/musrSim.tex +++ b/doc/musrSim.tex @@ -23,10 +23,11 @@ %\title{GEANT Simulation Program for the %$\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}\\ -$^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)}\\ %Institute of Physics AS\,CR\\ %%Academy of Sciences of the Czech Republic diff --git a/doc/musrSimAna.pdf b/doc/musrSimAna.pdf index 947ee3b..534e434 100644 Binary files a/doc/musrSimAna.pdf and b/doc/musrSimAna.pdf differ diff --git a/doc/musrSimAna.tex b/doc/musrSimAna.tex index 57efd1e..dc8c053 100644 --- a/doc/musrSimAna.tex +++ b/doc/musrSimAna.tex @@ -1,815 +1,861 @@ -\documentclass[twoside]{dis04} -\usepackage{epsfig} -\def\runauthor{PSI} -\def\shorttitle{musrSimAna} -\begin{document} - -\newcommand{\musr}{\ensuremath{\mu}SR} -\newcommand{\musrSim}{\emph{musrSim}} -\newcommand{\musrSimAna}{\emph{musrSimAna}} -\title{Manual of musrSimAna} -\author{Kamil Sedl\'ak$^1$} - -\address{{$^1$ Laboratory for Muon Spin Spectroscopy, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland}} - -\maketitle - -\abstracts{ -``musrSimAna'' is an analysis program, which helps the user to analyse and interpret the output of the -``musrSim'' simulation.} - -%======================================================================================================== -\section{Introduction} -\label{introduction} -The purpose of the \musrSimAna\ program is to analyse the Root tree, the output -of the \musrSim\ program. -In general, the event-by-event information stored in the Root tree can be easily used only for -a very quick and rough tests -- e.g.\ to see, where the muons stop and decay irrespective -of whether they were triggered in the M-counter or not, or to have an idea what energy -spectrum was deposited in a given counter. Typically, however, one is more interested -to study \musr\ signal amplitude, time-independent background, or where the muons decay -only for the ``stopped'' muons (e.g.\ muons triggered by M-counter and not detected -by a veto counter), or only for triggered muons which actually stopped in a cryostat -rather than in a sample. Such a study requires some kind of analysis program, which -loops over all events stored in the Root tree (output of \musrSim), and implements the -logic of the required \musr\ experiment. - -One way to do such analysis is to start from a Root command {\tt MakeClass(``MyAnalysis'')}, -which creates some sort of skeleton c++ class for the given Root tree. This allows the user -to make his/her own and very specific analysis program tailored to his/her needs. On the other -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 -program for \musr\ instruments. -At the moment, however, \musrSimAna\ is tuned to the needs of continuous muon beam facilities. -Some modifications might be necessary for the pulsed muon beams. - -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 -of the steering file was taken over from the TDC setup files used by real \musr\ instruments -at PSI~\cite{acquisitionProkscha,TDCsetup}. This setup file defines the ``logic'' of a given experiment, -namely the coincidences and anti-coincidences of different counters and veto detectors, as well as some -timing parameters. -The setup file for the case of simulation had to be extended by the definitions of -histograms, cuts, fitting, etc. The description of the setup file will be -presented in chapter~\ref{howToAnalyse}-\ref{description}, the chapter~\ref{sec:GPD} -illustrates the whole concept of \musrSimAna\ on the example of an existing \musr\ instrument. - - - -\section{The main components of the \musrSimAna\ setup file} -\label{howToAnalyse} -The parameters for the analysis are stored in the setup file ``{\tt RUNxxx.v1190}'', whose name -typically consists of the run number ({\tt RUN}), -some further identifier ({\tt xxx}), and ends with ``{\tt .v1190}''. -One of the following commands can be used to execute \musrSimAna: \\[1em] -{\tt -\$> musrSimAna RUN RUNxxx \\ -\$> musrSimAna RUN RUNxxx nographic > data/RUN.RUNxxx.out -}\\[1em] -where the first {\tt RUN} stands for the run number (more precisely, it specifies -the \musrSim\ output, which is expected to be stored as ``{\tt data/musr\_RUN.root}'') -and {\tt RUNxxx} specifies the \musrSimAna\ setup file without the ending, which is -expected to be stored in the current directory. The output file of \musrSimAna, -containing the result histograms, will be stored as ``{\tt data/his\_RUN\_RUNxxx.v1190.root}''. - -The syntax of the setup file in \musrSimAna\ is based on the experimental one, -defined in~\cite{TDCsetup}. -At the beginning of this file, some timing parameters are defined:\\[1em] -\begin{tabular}{ll} -{\tt RESOLUTION=100 } & \emph{\ldots one TDC bin corresponds to 100\,ps} \\[0.7em] -{\tt MCOINCIDENCEW=50 } & \emph{\ldots time interval (in TDC bins) to find coincidences with M-counter} \\ -{\tt PCOINCIDENCEW=50 } & \emph{\ldots time interval (in TDC bins) to find coincidences with P-counter} \\ -{\tt VCOINCIDENCEW=100 } & \emph{\ldots time interval (in TDC bins) to find anti-coincidences with M-counter} \\[0.7em] -{\tt MUONRATEFACTOR=0.089 } & \emph{\ldots will be described in section~\ref{sec:eventMixing}} \\[0.7em] -{\tt DATAWINDOWMIN=-2.0 } & \emph{\ldots data interval (in $\mu$s) in which positrons are detected} \\ -{\tt DATAWINDOWMAX=10.0 } & \\[0.7em] -{\tt PILEUPWINDOWMIN=-10.0 } & \emph{\ldots the pileup interval (in $\mu$s) for muons} \\ -{\tt PILEUPWINDOWMAX=10.0 } & \\[1em] -\end{tabular} -% -A good event has exactly one hit in the M-counter within the pile-up window, and exactly -one hit in the positron counters within the data window. Both windows are defined relative to $t_0$, -which is the time of the currently analysed M-counter hit. - -The coincidence and anti-coincidence logic between the different counters -(the detector logic) is also specified in the setup file. -An example may look like this:\\ -\begin{verbatim} -102; "M up"; M; 0.4; 2005; -401; - 1; "B1"; P; 0.4; 2005; 21 -401; B1; 1485; 1515; 50995; - 2; "B2"; P; 0.4; 2005; 22 -401; B2; 1485; 1515; 50995; - 11; "F1"; P; 0.4; 2005; -401 -21 -22; F11; 1485; 1515; 50995; - 12; "F2"; P; 0.4; 2005; -401 -21 -22; F12; 1485; 1515; 50995; - 13; "F3"; P; 0.4; 2005; -401 -21 -22; F13; 1485; 1515; 50995; - 21; "Coinc B1" K; 0.4; 2005; - 22; "Coinc B2" K; 0.4; 2005; -401; "Active Veto" V; 0.1; 2005; -\end{verbatim} -\mbox{} \\ -% -In each row, the first number stands for the detector ID (defined in the steering file -of the \musrSim), the second variable is the name of the counter, -the third variable identifies the type of the counter (M = muon, P = positron, -K = coincidence and V = veto counters, respectively), the forth variable is the threshold in MeV applied -to the energy deposited in the counter, the fifth variable is a time delay (in units -of bin-width) applied to the given counter, the sixth set of parameters defines -which other detectors have to be in coincidence (positive ID) or anti-coincidence -(negative ID) with the given counter, and in the case of positron counters, the -last four parameters define a special histogram for the given counter. (Note -that this histogram is used for a compatibility with the experimental steering file, however -there is a different and more powerful way of how to define histograms, which is -defined later). The only difference with respect to the experimental steering -file is the presence of threshold definition in the fourth column. -In the example above, the M-counter is the volume which has ID 102 in the simulation, -and a valid hit has to have energy in the counter above 0.4\,MeV, and no signal -above 0.1\,MeV in the ``Active Veto'' (ID=401) within a time interval defined by -{\tt VCOINCIDENCEW}. The detailed description of *.v1190 steering files can be found -in chapter~\ref{description} and in~\cite{TDCsetup}. - - - -The main output of the simulation is generated in the form of histograms. -The histograms are defined in the steering file. For example, the following -line defines 1-dimensional histogram of the $z$-position of where the muons stop and decay:\\[1em] -{\tt -musrTH1D hMuDecayPosZ "Where the muons stop;z[mm];N" 100 -5. 5. muDecayPosZ -}\\[1em] -This histogram has 100 bins, spans from -5\,mm to 5\,mm, and the variable filled -in the histogram is {\tt muDecayPosZ}. -In fact, this line does not create one histogram, but an array of histograms -- one -histogram for each ``condition''. An example of five conditions reads:\\[1em] -{\tt -condition 1 oncePerEvent \\ -condition 2 muonDecayedInSample\_gen \\ -condition 4 muonTriggered\_det \\ -condition 6 goodEvent\_det \\ -condition 9 pileupEvent -}\\[1em] -where integer numbers (1,2,4,6 and 9) denote the condition number, and the string at the end -describes the condition -- e.g.\ ``{\tt muonDecayedInSample\_gen}'' specifies that muon -has stopped and decayed in the sample, and ``{\tt muonTriggered\_det}'' specifies that -there had to be a good muon candidate triggered in the M-counter. -The ending ``{\tt \_gen}'' indicates that the variable used in the condition -was ``generated'', i.e.\ it is known in the simulation, but can not be directly -measured experimentally. -On the other hand the ending ``{\tt \_det}'' specifies that the given condition -is based on measurable variables, as e.g.\ energy deposits in a counter. - -There can be up to -30 conditions requested, and for each of them a separate histogram will be filled. -In the output file of \musrSimAna, the histograms corresponding to the previous -set of conditions would be saved as -{\tt hMuDecayPosZ\_1, hMuDecayPosZ\_2, hMuDecayPosZ\_4, hMuDecayPosZ\_6, hMuDecayPosZ\_9}. -The {\tt hMuDecayPosZ\_1} shows where the muons stop irrespective whether they were -detected or not. The {\tt hMuDecayPosZ\_6} contains ``good'' muons, i.e.\ muons that would -be saved in the final histograms in the experiment. The {\tt hMuDecayPosZ\_9} shows -the muons, which contribute to the time-independent background. Note that there are always -two muons, $\mu_1$ and $\mu_2$, contributing to the time-independent background, and {\tt hMuDecayPosZ\_9} -shows the one ($\mu_1$) that started the M-counter. More to this topic is presented -in chapter~\ref{sec:eventMixing}. Here we just mention that in order -to see the $z$-coordinate of the second muon ($\mu_2$), -which was not triggered by the M-counter but whose decay positron hit the positron counter, -the histogram {\tt hpileup\_muDecayPosZ\_9} must be used:\\[1em] -{\tt -musrTH1D hpileup\_muDecayPosZ "Pileup muons;z[mm];N" 100 -5.0 5. pileup\_muDecayPosZ -} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Event mixing -- an unavoidable complication} -\label{sec:eventMixing} -The output of \musrSim\ is stored in a Root tree named ``{\tt t1}''. One event corresponds -to one simulated muon, and it is saved as one row of the tree. The only information that -relates one event (muon) to any other one is the variable {\tt timeToNextEvent}, -which is a randomly generated time difference between two subsequent events. - -The \musrSimAna\ allows one to study the ``time-independent background'', which is -due to mixing of two different events. A simple example of the event mixing affecting -the \musr\ measurement is the following: the first muon, $\mu_1$, hits the M-counter, -and subsequently stops and decays in the sample, however the decay positron escapes -undetected -- most likely because of the limited angular acceptance of positron counters. -The second muon, $\mu_2$, arrives around the same time, misses the M-counter, -and stops and decays in a collimator wall or elsewhere. -Its decay positron hits a positron counter. Thus there are good-looking muon and positron -hits, which however arise from two uncorrelated muons. -This fake (background) event is experimentally treated as a good event, -and contributes to the time-independent background. - -Events can be mixed in \musrSimAna, allowing us to study sources of -the time-independent background. -%\emph{MusrSimAna} loops over the events and checks whether the conditions defining -%the detector logic in the steering file have been fulfilled. However, to take event -%mixing properly into account, the conditions have to be veryfied on several subsequent -%events. Therefore, at the beginning of the analysis of -Before analysing a given event, arrays of hits are filled -for all counters (M, positron, veto, coincidence counters), which store the hits occurring in the -future up to the time equal to $ 3\, \cdot $ pileup window or $ 3\, \cdot $ data window, -whatever is larger\footnote{The data and pile-up windows defined by parameters -{\tt DATAWINDOWMIN, DATAWINDOWMAX, PILEUPWINDOWMIN} and {\tt PILEUPWINDOWMAX} -are applied later on in the analysis.}. -There is one such array for every counter. After this initial filling, -there might be several hits in every array, originating from one or more events. - - -The fraction of the time-independent background to good events depends on the incoming muon rate -measured by the trigger, possibly in the anti-coincidence with a backward veto detector, if -used in the experiment. -Typically, the experimentalists set the incoming muon rate (rate of \emph{stopped muons}) -to $\sim 30\,000\,\mu/s$. The same should be done in the simulation. -However, the rate of stopped muons is known only after the analysis is done, -because, for example, many simulated muons stop and decay in collimators or elsewhere -in the beam-pipe without producing any signal in the M-counter. -Therefore the simulation is normally started with an initial rate of generated muons -of $30\,000\,\mu/s$, which in practise can correspond to much lower rate of \emph{stopped muons}. -The rate of stopped muons is calculated at the end of the \musrSimAna\ execution, and -it is printed out in the \musrSimAna\ output. The user can use this information and rescale the -initial muon rate by changing parameter {\tt MUONRATEFACTOR} in the {\tt *.v1190} setup file. -This is done without the necessity to re-run the CPU consuming -simulation -- only the \musrSimAna\ has to be repeated. The complete simulation and analysis -chain therefore usually consists of three steps: -% -\begin{enumerate} - \item \musrSim\ (takes typically 10 hours). - \item \musrSimAna\ with {\tt MUONRATEFACTOR=1.0} (takes typically 10 minutes). - \item \musrSimAna\ with {\tt MUONRATEFACTOR} determined in step 2. -\end{enumerate} -% -The {\tt MUONRATEFACTOR} specifies a multiplicative factor applied to -the variable {\tt timeToNextEvent} -(the randomly generated time difference between two subsequent events). - -{\bf IMPORTANT NOTE:} In order to get the pile-up effects correctly analysed by \musrSimAna, -it is probably necessary to run the \musrSim\ with no event reweighting (i.e.\ -the command ``{\tt /musr/command logicalVolumeToBeReweighted \ldots}'' must {\bf not} be used). -All events should/have to be (?) saved in the Root tree -(i.e.\ the command ``{\tt /musr/command storeOnlyEventsWithHits false}'' must be used). - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Detailed list of steering file parameters} -\label{description} -\begin{description} -% \item{\bf INSTRUMENT=\emph{string}} \\ -% ignored -% \item{\bf DESCRIPTION=\emph{string}} \\ -% ignored -% \item{\bf TYPE=\emph{string}} \\ -% ignored - \item{\bf RESOLUTION=\emph{value}} \\ - width of the TDC bin in picoseconds. - \item{\bf MDELAY=\emph{value}} \\ - currently not used (probably not needed in the case of simulation). - \item{\bf PDELAY=\emph{value}} \\ - currently not used (probably not needed in the case of simulation). - \item{\bf MCOINCIDENCEW=\emph{value}} \\ - time window for the coincidences of the coincidence detectors (``K'') - with the M-counter. The \emph{value} is given in TDC bins (see {\tt RESOLUTION} above). - \item{\bf PCOINCIDENCEW=\emph{value}} \\ - time window for the coincidences of the coincidence detectors (``K'') - with positron counters. The \emph{value} is given in TDC bins. - \item{\bf VCOINCIDENCEW=\emph{value}} \\ - time window for the coincidences of the anti-coincidence detectors (``V'') - with any other counter. The \emph{value} is given in TDC bins. - \item{\bf MUONRATEFACTOR=\emph{value}} \\ - a multiplicative factor which is used to rescale time differences between subsequent muons. - Setting \emph{value} larger than 1 artificially prolongs the time difference between - 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 - incoming muon rate to a given required value, typically to 30\,000\,$\mu/$s.\\ - See also variable ``INFINITELYLOWMUONRATE''. - \item{\bf INFINITELYLOWMUONRATE} \\ - 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 - no pileup can be observed. The variable ``MUONRATEFACTOR'' becomes irrelevant when - INFINITELYLOWMUONRATE is specified. - \item{\bf DATAWINDOWMIN=\emph{value}} \\ - Beginning of the data interval for the positron counters in $\mu$s. - \item{\bf DATAWINDOWMAX=\emph{value}} \\ - End of the data interval for the positron counters in $\mu$s. - \item{\bf PILEUPWINDOWMIN=\emph{value}} \\ - Beginning of the pileup interval for the M-counter in $\mu$s. - \item{\bf PILEUPWINDOWMAX=\emph{value}} \\ - End of the pileup interval for the M-counter in $\mu$s. - \item{\bf PROMPTPEAKMIN=\emph{value}} \\ - Beginning of the prompt-peak interval in $\mu$s. This variable is used only for the condition - ``{\tt promptPeak}'', ``{\tt promptPeakF}'', etc.\ , and normally does not need to be specified. It becomes useful if - the user wants to investigate, where the prompt-peak originates from (where do the muons, - which give rise to the prompt peak, stop). The default value is -0.01\,$\mu$s. - \item{\bf PROMPTPEAKMAX=\emph{value}} \\ - End of the prompt-peak interval in $\mu$s, the default value is 0.01\,$\mu$s. (See comments - for {\tt PROMPTPEAKMIN}.) - \item{\bf MUSRMODE=\emph{string}} \\ - Defines the mode of \musr\ experiment -- presently only ``D'', corresponding to - the time differential mode is implemented. - \item{\bf REWINDTIMEBINS=\emph{value}} \\ - 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. - \item{\bf DEBUGEVENT \emph{eventID} \emph{debugLevel}}\\ - Prints out debug information about the event with the ID \emph{``eventID''}. - The higher the \emph{debugLevel}, the more details are printed. - (Both \emph{eventID} and \emph{debugLevel} are integers). - \item{\bf CLONECHANNEL \emph{ID1} \emph{ID2}}\\ - Clones the hits detected in counter \emph{ID1} into a new counter \emph{ID2}. - A different (even a lower) threshold can be applied to the cloned counter. - This way the same counter can be used as two independent counters -- e.g.\ once as a veto - detector for the M-counter, and simultaneously as a coincidence detector for - a P-counter. In both cases the energy threshold and time windows are defined independently. - \item{\bf WRITE\_OUT\_DUMP\_FILE \emph{fileNameString} \emph{clockChannel} \emph{randomJitter}}\\ - If present, this command will create two output files, the so-called ``dump files'':\\ - {\tt data/TDC\_V1190\_200\_dump\_\emph{fileNameString}.bin} -- file that can be used - as an input to the PSI analysis front-end of a real experiment.\\ - {\tt data/TDC\_V1190\_200\_dump\_\emph{fileNameString}.txt} -- file that contains the same - information (hits) as the previous file, however in a human-readable form. The first number in the file - stands for the channel number, the second number stands for the time bin in the TDC bin units.\\ - \emph{clockChannel} ... the channel of the clock signal (typically 15).\\ - \emph{randomJitter} ... this value is in TDC bins, typically 0 or 8000 (?). If \emph{randomJitter} is smaller then - 1, then the hits in the dump files will be sorted according to time. If it is larger than 0, then - subsequent hits can be unordered in time, but the time difference never exceeds the value of \emph{randomJitter}. - This is just a technical thing serving to test the analysis software -- it should not - have any effect on the analysis results. - \item{\bf musrTH1D \emph{histoName} \emph{histoTitle} \emph{nBins} \emph{min} \emph{max} \emph{variable} - [{\tt rotreference} $\nu_{\rm RRF}$ $\phi_{\rm RRF}$] $|$ [correctexpdecay]} \\ - Defines a histogram (or more precisely an array of histograms, where the number of histograms - in the array is given by the number of conditions, see section~\ref{howToAnalyse}). - The name of the histogram is defined by \emph{histoName} + the number of the condition. - The string variable \emph{histoTitle} specifies the title of the histogram, - \emph{nBins}, \emph{min} and \emph{max} stand for the number of bins and minimum and maximum - of the $x$-axis of the histogram. \\ - The optional keyword ``{\tt rotreference}'' signals that the given histogram will be filled in - rotating reference frame (RRF) with the frequency of $\nu_{\rm RRF}$ and a phase shift of $\phi_{\rm RRF}$. - \\ - The optional keyword ``{\tt correctexpdecay}'' signals that the given histogram will be corrected - for the muon exponential decay (i.e. multiplied by a factor $\exp(t/2.19703)$. It is meaningful - only for time variables. - \\ - The \emph{variable} stands for the variable that will be - filled into the histogram. The \emph{variable} can be any variable from the output Root tree - of musrSim (see ``Manual of musrSim'') (except for the array variables like - {\tt det\_*[], save*[], odet\_*[]}). In addition, it can be also one of the following: - \begin{description} - \item[muDecayPosR] \ldots $\sqrt{ {\rm muDecayPosX}^2 + {\rm muDecayPosY}^2}$. - \item[wght] \ldots weight of the event. - \item[det\_m0edep] \ldots energy deposited in the M-counter that gives the muon signal. - \item[det\_posEdep] \ldots energy deposited in the P-counter that gives the positron signal. - \item[muIniPosR] \ldots $\sqrt{ {\rm muIniPosX}^2 + {\rm muIniPosY}^2}$. - \item[muIniMomTrans] \ldots $\sqrt{ {\rm muIniMomX}^2 + {\rm muIniMomY}^2}$. - \item[muTargetPol\_Theta] \ldots theta angle of the muon spin when muon enters target (-180,180 deg). - \item[muTargetPol\_Theta360]\ldots theta angle of the muon spin when muon enters target (0,360 deg). - \item[muTargetPol\_Phi] \ldots phi angle of the muon spin when muon enters target (-180,180 deg). - \item[muTargetPol\_Phi360] \ldots phi angle of the muon spin when muon enters target (0,360 deg). - \item[pos\_Momentum] \ldots magnitude of the momentum of the decay positron (``generated'', not measurable variable). - \item[pos\_Trans\_Momentum] \ldots transverse momentum of the decay positron. - \item[pos\_Radius] \ldots positron radius calculated from the decay positron momentum and magnetic - field at the point of decay. - \item[pos\_Theta] \ldots polar angle of the decay positron. - \item[pos\_Phi] \ldots azimuth angle of the decay positron. - \item[det\_time10] \ldots time difference between the positron and muon counters - (measured by the respective counters). - \item[gen\_time10] \ldots the time difference between the muon decay and the first muon hit - in the M-counter (i.e.\ {\tt muDecayTime - muM0Time}). - \item[det\_time10\_MINUS\_gen\_time10] \ldots {\tt det\_time10 - gen\_time10} in picoseconds. - \item[det\_time20] \ldots very similar to {\tt det\_time10}, however taking into - account ``counter phase shift'' defined by {\tt counterPhaseShifts} - variable. This gives the user a possibility to sum up backward and - forward histograms into one histogram (this of course make sense - only in the simulation, where there is ``no physics'' happening - in the sample, just the muon spin rotation). - \item[pileup\_eventID] \ldots eventID of the $\mu_2$, where $\mu_2$ stands for the muon, - which did not give signal in the M-counter, but whose - decay positron gave signal in the positron counter around the time - when a different muon ($\mu_1$) triggered the M-counter. - \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\_muDecayPosR] \ldots radius of where the $\mu_2$ stopped and decayed. - \end{description} - 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). - \item{\bf musrTH2D \emph{histoName} \emph{histoTitle} \emph{nBins} \emph{min} \emph{max} \emph{nBins2} \emph{min2} \emph{max2} \emph{variable}} \\ - Similar to \emph{musrTH1D}, but for a 2-D histogram. - \item{\bf humanDecayHistograms \emph{hist\_decay\_detID} \emph{hist\_decay\_detID\_pileup} \emph{id$_1$} \emph{name$_1$} \ldots \emph{id$_n$} \emph{name$_n$} } \\ - This is a special kind of histogram, which converts two histograms - (\emph{hist\_decay\_detID} \emph{hist\_decay\_detID\_pileup}) - into a human-friendly histograms, where detector ID on the $x$-axis is converted into a string label. - The \emph{id$_i$} is the detector id, and the \emph{name$_i$} is the corresponding label name. - If \emph{name$_i$ = name$_j$}, the corresponding bins of the original histograms will - be summed up together into the same bin. - The \emph{hist\_decay\_detID} and \emph{hist\_decay\_detID\_pileup} have to be defined before - (by the command {\tt musrTH1D}). - \item{\bf condition \emph{conditionID} \emph{conditionName}} \\ - Definition of a condition, which is then used when filling histograms. The \emph{conditionID} - specifies the number of the condition, and must be between 0 and 30 (0 and 30 are also possible). - The \emph{conditionName} is one of the following: - \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[oncePerEvent] \ldots true once for every event (the first hit in M-counter, if any, is considered). - \item[muonDecayedInSample\_gen] \ldots true if muon stopped and decayed in the sample. - \item[muonTriggered\_gen] \ldots true if muon passed through the M-counter (irrespective of the deposited energy) - -- not a measurable variable. - \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. - \item[positronHit\_det] \ldots true if a good positron candidate was found in the positron counter. - Double hits within the data window are excluded. - \item[goodEvent\_det] \ldots true if {\tt muonTriggered\_det} and {\tt positronHit\_det}. - \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 - on the positron is implied, i.e.\ the positron may or may not be detected. - Not a measurable variable. - \item[goodEvent\_det\_AND\_goodEvent\_gen] - \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[goodEvent\_det\_AND\_muonDecayedInSample\_gen] - \item[goodEvent\_F\_det] \ldots {\tt goodEvent\_det}, where the positron was detected in the forward detectors - defined by the command {\tt counterGrouping}. - \item[goodEvent\_B\_det] \ldots like {\tt goodEvent\_F\_det} but for backward positron counters. - \item[goodEvent\_U\_det] \ldots like {\tt goodEvent\_F\_det} but for upper 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\_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\_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\_D\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_D\_det} and {\tt pileupEvent}. - \item[goodEvent\_L\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_L\_det} and {\tt pileupEvent}. - \item[goodEvent\_R\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_R\_det} and {\tt pileupEvent}. - \item[promptPeak] \ldots {\tt goodEvent\_det}, and {\tt PROMPTPEAKMIN < det\_time10 < PROMPTPEAKMAX}. - \item[promptPeakF] \ldots like {\tt goodEvent\_F\_det} and {\tt promptPeak}. - \item[promptPeakB] \ldots like {\tt goodEvent\_B\_det} and {\tt promptPeak}. - \item[promptPeakU] \ldots like {\tt goodEvent\_U\_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[promptPeakR] \ldots like {\tt goodEvent\_R\_det} and {\tt promptPeak}. - \end{description} - Additional conditions may be implemented on request. - \item{\bf draw \emph{histogramName} \emph{conditionID} } \\ - Plot histogram (for a given condition) at the end of the analysis. Note that all histograms - are saved into the output file irrespective whether they are plotted or not. - \item{\bf counterPhaseShifts \emph{ID$_1$} \emph{$\phi_1$} \emph{ID$_2$} \emph{$\phi_2$} \ldots \emph{ID$_n$} \emph{$\phi_n$} }\\ - Defines relative phase shifts of signal between different positron counters, which is used - for calculation variable {\tt det\_time20}. \emph{ID$_i$} is the ID number of the positron counter, - \emph{$\phi_i$} is its phase shift. This gives the user a possibility to sum up backward and - forward histograms into one histogram. - \item{\bf counterGrouping \emph{group} \emph{ID$_1$ ID$_2$ \ldots ID$_n$} } \\ - This defines a group of detectors, where \emph{group} stands for ``B'' (backward), - ``F'' (forward), ``U'' (up), ``D'' (down), ``L'' (left) and ``R'' (right) detectors. - This grouping is used in the definition of some conditions. - \item{\bf sampleID \emph{ID$_1$ ID$_2$ \ldots ID$_n$} } \\ - Defines which volume (or volumes, if there are more) is the sample. Typically, the - sample is just one volume, but sometimes there is a smaller volume inside of the sample, - to which a field is applied, so the small volume also has to be considered as the sample. - This information is needed for the condition ``muonDecayedInSample\_gen''. - \item{\bf setSpecialAnticoincidenceTimeWindow \emph{detectorID} \emph{timeMin} \emph{timeMax} \emph{unit}} \\ - This command sets a special anti-coincidence time window for a detector \emph{detectorID}. - Normally, the anti-coincidence time window is defined by {\tt VCOINCIDENCEW}, and is the same for all anti-coincidence - detectors. However, sometimes it might be interesting to set the anti-coincidence time window - differently for a specific detector (e.g.\ one might test an anti-coincidence of a veto detector with - the M-counter for the whole pile-up time window of $\sim$\,10\,$\mu s$. - Unlike in the case of {\tt VCOINCIDENCEW}, here the \emph{units} are not TDC bins, but - rather time in ``nanosecond'' or ``microsecond''. - \item{\bf fit \emph{histogramName} \emph{function} \emph{option} \emph{min} \emph{max} \emph{p$_1$} \ldots \emph{p$_n$}} \\ - Fits the histogram by a given function, where \emph{min}, \emph{max} define the range of the fit - on the $x$-axis of the histogram, \emph{option} is a string defining fit options (see Root manual for details), - and \emph{p$_1$} \ldots \emph{p$_n$} are (typically, with some exceptions) - the initial values of the function parameters. The following functions are currently predefined: - \begin{description} - \item[pol0] $=p_0$ \ldots a constant (1 parameter) - typically used to fit background. - \item[simpleExpoPLUSconst] $=p_0 \exp(-x/2.19703)+p_1$ - \item[rotFrameTime20] $= p_2 \cos(p_0 x+p_1)$ - \item[funct1] $=p_3 \exp((p_4 - x)/2.19703) \cdot (1+p_2 \cos(p_0 x+p_1))$ - \item[funct2] $=p_3 \exp((p_4 - x)/2.19703) \cdot (1+p_2 \cos(p_0 x+p_1)) + p_5$ - % \item[funct3] the same as {\tt funct2} - \item[funct4] $=p_3 \exp((- x)/2.19703) \cdot (1+p_2 \cos(p_0 x+p_1)) + p_4$ - \item[TFieldCos] $=p_3 (1+p_2 \cos(p_0 x + p_1))$ \hspace{1cm} (this function is useful when the histogram is filled with {\tt ``correctexpdecay''} keyword.) - \item[TFieldCosPLUSbg] $=p_3 (1+p_2 \cos(p_0 x + p_1)) + p_4 \exp(x/2.19703)$ - \hspace{1cm} (this function is useful when the histogram is filled with {\tt ``correctexpdecay''} keyword.) - \item[gaus] ... Gauss distribution - \end{description} - - -\end{description} -%======================================================================================================== -\section{A real-life example: GPD instrument} -\label{sec:GPD} -The simulation of the General Purpose Decay-Channel Spectrometer (GPD) -instrument~\cite{GPD} installed at PSI has been exemplified -in the \musrSim\ manual~\cite{musrSim}. -Here we analyse the output of this simulation using \musrSimAna. -The run number of this simulation is 201, therefore the steering -file names are ``{\tt 201.mac}'' for \musrSim, and ``{\tt 201.v1190}'' for -\musrSimAna, respectively, and the output file name of \musrSim\ is saved as -``{\tt data/musr\_201.root}''. -The detector system is very simple with only six counters -- M-counter, -two backward positron counters and three forward positron counters. -The reader is strongly recommended to see the illustration of the GPD -geometry in the \musrSim\ manual~\cite{musrSim}. - -8\,000\,000 of events were simulated -(i.e.\ 8\,000\,000 of muons were generated 100\,cm in front of -the GPD sample). -In only 949\,759 events (11.9\% out of the 8 million) there was a signal detected -in one or more counters. The remaining muons stopped somewhere (most -often in collimator, as we will see later), decayed, and the decay positron -(and any other particles created thereof) missed the counters. -This is illustrated in more details in Fig.~\ref{det_n}, -% -\begin{figure}[tbp]\centering -\epsfig{file=pict/det_n_infititely_low_muon_rate.eps,width=0.8\linewidth,clip=} -\caption{Number of hits in all counters per event, assuming infinitely low incoming muon -rate. The same detector may be hit more than once (e.g.\ if both the muon and its decay -positron pass through the M-counter).} -\label{det_n} -\end{figure} -% -where number of detector hit per event, assuming infinitely low incoming muon -rate, is shown. This plot was created in Root by executing:\\[1em] -{\tt -root [0] TFile* ff=new TFile("data/musr\_201.root") \\ -root [1] t1->Print() \\ -root [2] t1->Print("det\_n","det\_n>0") \\ -}\\[1em] -% -It has to be pointed out, that the ratio of muons passing through the opening -in collimators to the number of all generated muons strongly depends on the -beam properties -- beam profile, beam convergence, etc. Typically, if we have -too broad muon beam, we simulate -many ``useless'' events. However, the other extreme (simulating too narrow -beam) can lead to underestimating the time-independent background. - -It took approximately 12 hours of the CPU (on PC bought in 2010, where 1 out -of 4 processor cores was running) to simulate these 8\,000\,000 events. -Assuming the 30\,000\,$\mu/$s trigger rate, this corresponds to 26 seconds -of real experimental running. - -\subsection{Where the muons stop and decay} -\label{sect_muons} -The positions, or more precisely the components of the GPD instrument, where the muons -stop and decay, are shown in Fig.~\ref{humanDecayHistograms_1}: -% -\begin{figure}[tbp]\centering -\epsfig{file=pict/Plot201_1.eps,width=0.9\linewidth,% -%%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt,clip=} -clip=} -\caption{This plot indicates, where the muons stopped and decayed. -The dashed histogram shows all generated muons. The full-line histograms show -where stopped the muons, for which either the muon itself or its secondary -particle ($e^+, \gamma$) triggered the M-counter: black histogram stands for -all such muons, corresponding to infinitely low incoming muon rate, while -the red histogram stands for the incoming muon rate of 30\,000\,$\mu/$s. -8\,000\,000 of events were simulated.} -\label{humanDecayHistograms_1} -\end{figure} -% -%Notes to Fig.~\ref{humanDecayHistograms_1}: -\begin{itemize} - \item Figure~\ref{humanDecayHistograms_1} was generated by Root macro file ``Plot201.C''. - \item The labels on the $x$-axis are defined in the file {\tt 201.v1190} by the - command \\ - {\tt humanDecayHistograms \ldots} - \item The dashed-line histogram in Fig.~\ref{humanDecayHistograms_1} - shows where the muons stopped and decayed if no preselection - criteria are applied on the muons, i.e.\ if all generated muons are considered. - This is histogram ``{\tt humanDecayHistograms\_1}''. - \item The full-line histograms show - where stopped the muons, for which either the muon itself or its secondary - particle ($e^+, \gamma$) triggered the M-counter: black histogram stands for - all such muons, corresponding to infinitely low incoming muon rate, while - the red histogram represents the case for the 30\,000\,$\mu/$s incoming muon rate. - An energy deposit of at least 0.4\,MeV in the M-counter is required to fire the trigger. - The number of triggered events decreases with the incoming muon rate, - because some of the events are rejected due to the 10\,$\mu$s pileup gate. - - The histogram name is in both cases ``{\tt humanDecayHistograms\_4}'', - where the black histogram was calculated using the setup file ``{\tt 201a.v1190}'' - with the keyword {\tt INFINITELYLOWMUONRATE}, while the red histogram - was calculated using the setup file ``{\tt 201.v1190}'' - with {\tt MUONRATEFACTOR=0.0965819}. - \item The $\pm 10\,\mu$s pile-up gate at the incoming muon rate of 30\,000$\,\mu/$s - rejects approx.\ 45\% of the triggered events. This number can be calulated - in Root as the ratio of the ``{\tt Integral}'' of the red and black histograms - in Fig.~\ref{humanDecayHistograms_1}:\\[1em] - {\tt \small - root [0] TFile* file1 = new TFile("data/his\_201\_201a.v1190.root") \\ - root [1] humanDecayHistograms\_4->Integral() \\ - root [0] TFile* file2 = new TFile("data/his\_201\_201.v1190.root") \\ - root [1] humanDecayHistograms\_4->Integral() \\ - } - \item The muon sample fraction (ratio of muons stopped in the sample over all muons that fired - the trigger) for the triggered events (full-line histograms) - is 65\%, and it is practically the same - for both infinitely low and 30\,000\,$\mu/$s incoming rate. - This number can be obtained in Root by dividing the first column of histogram - {\tt humanDecayHistograms\_4} by the sum of all entries in this histogram:\\[1em] - {\tt \small - root [0] TFile* file = new TFile("data/his\_201\_201.v1190.root") \\ - root [1] (humanDecayHistograms\_4->GetBinContent(1))/(humanDecayHistograms\_4->Integral()) \\ - } - \item The largest fraction of generated muons (dashed-line histogram) stopped in collimators. - Only a small fraction of them caused a hit in the M-counter (full-line histograms). - \item Despite the high initial muon momentum of $100 \pm 3\,$GeV/c, muons are - significantly scattered in the last 50\,cm region of air. This can be - clearly seen if the magnetic field is off and a point-like muon beam - is used (which can be done by modifying the {\tt 201.mac} file) - -- only 77\% of the muons stop in the sample cell or in the sample, while the - remaining 23\% of the mouns are scattered so much in the air, that they - end up in collimators or elsewhere (not shown here). - \item ``World'' in the histogram label means that the muon decayed in the beampipe vacuum - or somewhere else in the air (on the fly). - \item ``Escaped'' means that the muon left the simulated instrument (more precisely the - ``world'' volume) prior its decay. -\end{itemize} - - -Figure~\ref{humanDecayHistograms_9} shows the ``pile-up events''. -% -\begin{figure}[htbp]\centering -\epsfig{file=pict/Plot201_2_new.eps,width=0.9\linewidth,% -%%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt,clip=} -clip=} -\caption{Pile-up events, i.e.\ the events in which one muon ($\mu_1$) fired the -trigger, while the hit in a positron counter is due to a decay positron from -a different muon ($\mu_2$). Pile-up events look like a good events, and contribute -to the time-independent background.} -\label{humanDecayHistograms_9} -\end{figure} -% -These are events, in which one muon ($\mu_1$) is triggered by the m-counter, -while a positron from a different muon ($\mu_2$) was detected by -a positron counter\footnote{In fact, the trigger may also be triggered by -the decay positron of $\mu_1$ and/or a positron counter may detect -directly $\mu_2$, not its decay positron. Such cases are rare, but they -are implicitly included in the simulation.}. -In addition to this requirement, the decay positron of $\mu_1$ must -escape undetected (e.g.\ it must miss positron counters) and $\mu_2$ must not trigger the m-counter --- otherwise the event would be rejected. -Pile-up events are the source of the time independent background. -Usually $\mu_1$ is a good-looking muon that stops in the sample or in the sample cell -(red histogram in Fig.~\ref{humanDecayHistograms_9}), while $\mu_2$ stops and decays at different places, -mainly in the collimators (green histogram in Fig.~\ref{humanDecayHistograms_9}). - -A nice visualisation of where the background-contributing muons $\mu_2$ stop and decay -is presented in Fig.~\ref{Pileup_muon_decay_map} (histogram ``{\tt hMuDecayMappileup\_9}''). -% -\begin{figure}[htbp]\centering -\epsfig{file=pict/Pileup_muon_decay_map.eps,width=0.7\linewidth,% -%%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt,clip=} -clip=} -\caption{Positions of where the $\mu_2$ stop and decay.} -\label{Pileup_muon_decay_map} -\end{figure} -% -In this two dimensional histogram, different components of the GPD instrument, -like the lead collimator, the copper collimator and the sample cell, can be recognised. -The lead collimator is located at the $z$-position between -115\,mm and -85\,mm. -Due to the high initial muon momentum of $\sim 100\,$MeV/c, -the maximum of muons in Fig.~\ref{Pileup_muon_decay_map} stop quite deep in the -lead collimator, at around $z=-103$\,mm. This might be a little bit surprising to the -\musr\ scientists who are used to work with the surface muons with momentum of 28\,MeV/c. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{The $\mu$SR signal} -% -Figure~\ref{hdet_time10_10} -% -\begin{figure}[htbp]\centering -\epsfig{file=pict/hdet_time10_10.eps,width=0.495\linewidth,% -bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} -\epsfig{file=pict/hdet_time10_11.eps,width=0.495\linewidth,% -bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} -%%clip=} -\caption{MuSR signal for the run 201 (TF$=300\,$gauss). The tree forward positron counters -are summed up in the left histogram, and the two backward counters in the right histogram.} -\label{hdet_time10_10} -\end{figure} -% -shows the $\mu$SR spectra for the same run, -i.e.\ for the transverse field of 300\,gauss, integrated over the three forward positron -counters (left histogram called {\tt hdet\_time10\_10}) -and over the two backward positron counters (right histogram called {\tt hdet\_time10\_11}). -Zero on the time axis corresponds to $t_0$, i.e.\ time of the m-counter hit. -One can see a prompt peak at $t_0$, time independent background at negative times -and an oscillating signal at positive times. -The following function has been fitted to the oscillating part of the signal: -% -\begin{equation} -f=p_3 \cdot e^{-t/2.19703} \cdot (1+p_2 \cdot \cos(t \cdot p_0+p_1))+p_4 -\label{eq_simple} -\end{equation} -The fits were restricted to the time interval of $(t_0+0.05 \mu\rm{s},t_0+9.95\mu\rm{s})$, -and the parameter $p_0$ was fixed (e.g. not fitted). -The fitted amplitude of asymmetry are $p_2 = 0.307 \pm 0.009$ and -$p_2 = 0.290 \pm 0.009$ for the forward and backward counters respectively. - -Parts of the spectra from Fig.~\ref{hdet_time10_10} are shown -in detail in Fig.~\ref{hdet_time10_10_detail}. -% -\begin{figure}[htbp]\centering -\epsfig{file=pict/hdet_time10_10_detail.eps,width=0.495\linewidth,% -bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} -\epsfig{file=pict/hdet_time10_11_pileup.eps,width=0.495\linewidth,% -bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} -%%clip=} -\caption{MuSR signal for the run 201 (TF$=300\,$gauss) -- details of -Fig.~\ref{hdet_time10_10}. The left plot shows the signal in the forward counters around $t_0$, -the right plot shows the (time-independent background) signal at negative times in the -backward counters.} -\label{hdet_time10_10_detail} -\end{figure} -% -The left plot in Fig.~\ref{hdet_time10_10_detail} shows the signal -in the forward counters around $t_0$, the right plot shows the -(time-independent background) signal at negative times in the backward counters. - -An important characteristic of a \musr\ instrument is the time-independent -background. It is usually expressed as -% -\begin{equation} -{\rm Bgr} = p_{-} / p_3 ~~~, -\label{eq_background} -\end{equation} -% -where $p_{-}$ is the fit to the time-independent background, i.e.\ signal at negative times, -and $p_3$ is the parameter from eq.(\ref{eq_simple}), which specifies what the -size of the signal would be at $t_0$ in the absence of oscillations. -In the case of backward counters ${\rm Bgr}_{\rm backw} = 14.47/262 = 5.5\,\%$, -in the case of forward counters ${\rm Bgr}_{\rm forw} = 6.88/267.9 = 2.6\,\%$. - -Note that the histogram on right hand side of Fig.~\ref{hdet_time10_10_detail} -is labelled ``{\tt hdet\_time10\_Bgr\_11}'', not ``{\tt hdet\_time10\_11}''. -In fact, the two histograms are identical, as one can see in the setup file -{\tt 201.v1190}. The only difference is in the fitting -- the same data stored in -both histograms are fitted by different functions in different time ranges. - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{The $\mu$SR signal from individual counters} -% -Figure~\ref{F11} shows the observed signal in the -forward counter\ No.~11 (FW11). -% -\begin{figure}[htbp]\centering -\epsfig{file=pict/F11_rebinned.eps,width=0.495\linewidth,% -bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} -\epsfig{file=pict/F11_B11_prompt_peak_thicker.eps,width=0.495\linewidth,% -bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} -\caption{\musr\ signal in the forward positron counter\ No.~11 (run 201, TF$=300\,$gauss). -The left plot shows the (rebinned) signal in the counter, -the right plot shows the detail of the \emph{prompt peak}, i.e.\ the region -around $t_0$ in the same counter (black line), -compared with the prompt peak in the backward positron counter\ No.~1 (magenta line).} -\label{F11} -\end{figure} -% -Originally, the histogram F11 was defined with the bin width of 100\,ps. -The number of bins was 50995, covering the time interval of approx.\ (-0.2\,$\mu$s, 4.9\,$\mu$s). -In the left hand side plot, however, the histogram was rebinned -(200 bins were summed up into 1 bin). -The right hand side plot shows the detail of the \emph{prompt peak}, i.e.\ the region -around $t_0$, of one forward and one backward positron counters, prior to the rebinning. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Conclusion of the GPD analysis example} -% -The purpose of the example analysis of the GPD simulation was to illustrate -the potential of \musrSim\ and \musrSimAna\ programs to investigate features -like time-independent background, sample muon fraction, prompt peak, \ldots -This information can be used in design and optimisation of \musr\ instruments. -%======================================================================================================== -\section{GPS instrument} -% -It is foreseen that GPS instrument could be arranged in two geometries -after the upgrade (depending from which side the calorimeter would be -inserted). -\begin{itemize} - \item {\tt 50130hb.v1190} -- Calorimeter inserted from one side. - \item {\tt 50130hl.v1190} -- Calorimeter inserted from the other side. - \item {\tt 50130hb1.v1190 -- 50130hb6.v1190} -- All positron counters - analysed individually. -\end{itemize} -See the document about the GPS simulations saved in the directory: \\ -/afs/psi.ch/project/lmu/Facility/musr\_simulations/documentation/GPS/ \\ -for more details. -%======================================================================================================== -\section{Other Examples} -Many different ``*.v1190'' files are stored in the file: -``run\_musrSimAna\_many\_files.tar.gz''. They could serve as additional examples. -Note that the syntax of the ``fit'' command was changed -at some point, and therefore the ``fit'' command might cause problems -(the {\tt ``option''} has to be added in the old ``*.v1190'' files). -%======================================================================================================== -\begin{thebibliography}{0} - -\bibitem{acquisitionProkscha} T.~Prokscha {\it et al.} ``A novel VME based \musr\ data acquisition system at PSI'', -Physica {\bf B~404}, (2009) 1007-1009. - -\bibitem{TDCsetup} ``TDC Manual -- Setting up the required logic'', -http://lmu.web.psi.ch/facilities/electronics/TDC/set\_logic.html - -\bibitem{GPD} -http://lmu.web.psi.ch/facilities/gpd/gpd.html - -\bibitem{musrSim} -K.Sedlak {\it et al.}, ``Manual of musrSim''. - - -\end{thebibliography} - -\end{document} +\documentclass[twoside]{dis04} +\usepackage{epsfig} +\def\runauthor{PSI} +\def\shorttitle{musrSimAna} +\begin{document} + +\newcommand{\musr}{\ensuremath{\mu}SR} +\newcommand{\musrSim}{\emph{musrSim}} +\newcommand{\musrSimAna}{\emph{musrSimAna}} +\title{Manual of musrSimAna} +\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{{$^2$ ISIS Facility, Rutherford Appleton Laboratory, Chilton, Oxon OX11 0QX, U.K.}} + +\maketitle + +\abstracts{ +``musrSimAna'' is an analysis program, which helps the user to analyse and interpret the output of the +``musrSim'' simulation.} + +%======================================================================================================== +\section{Introduction} +\label{introduction} +The purpose of the \musrSimAna\ program is to analyse the Root tree, the output +of the \musrSim\ program. +In general, the event-by-event information stored in the Root tree can be easily used only for +a very quick and rough tests -- e.g.\ to see, where the muons stop and decay irrespective +of whether they were triggered in the M-counter or not, or to have an idea what energy +spectrum was deposited in a given counter. Typically, however, one is more interested +to study \musr\ signal amplitude, time-independent background, or where the muons decay +only for the ``stopped'' muons (e.g.\ muons triggered by M-counter and not detected +by a veto counter), or only for triggered muons which actually stopped in a cryostat +rather than in a sample. Such a study requires some kind of analysis program, which +loops over all events stored in the Root tree (output of \musrSim), and implements the +logic of the required \musr\ experiment. + +One way to do such analysis is to start from a Root command {\tt MakeClass(``MyAnalysis'')}, +which creates some sort of skeleton c++ class for the given Root tree. This allows the user +to make his/her own and very specific analysis program tailored to his/her needs. On the other +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 +program for \musr\ instruments. +\musrSimAna\ was originally written for continuous muon beam facilities. +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 +required by the analysis in an input text file (steering file). The initial structure +of the steering file was taken over from the TDC setup files used by real \musr\ instruments +at PSI~\cite{acquisitionProkscha,TDCsetup}. This setup file defines the ``logic'' of a given experiment, +namely the coincidences and anti-coincidences of different counters and veto detectors, as well as some +timing parameters. +The setup file for the case of simulation had to be extended by the definitions of +histograms, cuts, fitting, etc. The description of the setup file will be +presented in chapter~\ref{howToAnalyse}-\ref{description}, the chapter~\ref{sec:GPD} +illustrates the whole concept of \musrSimAna\ on the example of an existing \musr\ instrument. + + + +\section{The main components of the \musrSimAna\ setup file} +\label{howToAnalyse} +The parameters for the analysis are stored in the setup file ``{\tt RUNxxx.v1190}'', whose name +typically consists of the run number ({\tt RUN}), +some further identifier ({\tt xxx}), and ends with ``{\tt .v1190}''. +One of the following commands can be used to execute \musrSimAna: \\[1em] +{\tt +\$> musrSimAna RUN RUNxxx \\ +\$> musrSimAna RUN RUNxxx nographic > data/RUN.RUNxxx.out +}\\[1em] +where the first {\tt RUN} stands for the run number (more precisely, it specifies +the \musrSim\ output, which is expected to be stored as ``{\tt data/musr\_RUN.root}'') +and {\tt RUNxxx} specifies the \musrSimAna\ setup file without the ending, which is +expected to be stored in the current directory. The output file of \musrSimAna, +containing the result histograms, will be stored as ``{\tt data/his\_RUN\_RUNxxx.v1190.root}''. + +The syntax of the setup file in \musrSimAna\ is based on the experimental one, +defined in~\cite{TDCsetup}. +At the beginning of this file, some timing parameters are defined:\\[1em] +\begin{tabular}{ll} +{\tt RESOLUTION=100 } & \emph{\ldots one TDC bin corresponds to 100\,ps} \\[0.7em] +{\tt MCOINCIDENCEW=50 } & \emph{\ldots time interval (in TDC bins) to find coincidences with M-counter} \\ +{\tt PCOINCIDENCEW=50 } & \emph{\ldots time interval (in TDC bins) to find coincidences with P-counter} \\ +{\tt VCOINCIDENCEW=100 } & \emph{\ldots time interval (in TDC bins) to find anti-coincidences with M-counter} \\[0.7em] +{\tt MUONRATEFACTOR=0.089 } & \emph{\ldots will be described in section~\ref{sec:eventMixing}} \\[0.7em] +{\tt DATAWINDOWMIN=-2.0 } & \emph{\ldots data interval (in $\mu$s) in which positrons are detected} \\ +{\tt DATAWINDOWMAX=10.0 } & \\[0.7em] +{\tt PILEUPWINDOWMIN=-10.0 } & \emph{\ldots the pileup interval (in $\mu$s) for muons} \\ +{\tt PILEUPWINDOWMAX=10.0 } & \\[1em] +\end{tabular} +% +A good event has exactly one hit in the M-counter within the pile-up window, and exactly +one hit in the positron counters within the data window. Both windows are defined relative to $t_0$, +which is the time of the currently analysed M-counter hit. + +The coincidence and anti-coincidence logic between the different counters +(the detector logic) is also specified in the setup file. +An example may look like this:\\ +\begin{verbatim} +102; "M up"; M; 0.4; 2005; -401; + 1; "B1"; P; 0.4; 2005; 21 -401; B1; 1485; 1515; 50995; + 2; "B2"; P; 0.4; 2005; 22 -401; B2; 1485; 1515; 50995; + 11; "F1"; P; 0.4; 2005; -401 -21 -22; F11; 1485; 1515; 50995; + 12; "F2"; P; 0.4; 2005; -401 -21 -22; F12; 1485; 1515; 50995; + 13; "F3"; P; 0.4; 2005; -401 -21 -22; F13; 1485; 1515; 50995; + 21; "Coinc B1" K; 0.4; 2005; + 22; "Coinc B2" K; 0.4; 2005; +401; "Active Veto" V; 0.1; 2005; +\end{verbatim} +\mbox{} \\ +% +In each row, the first number stands for the detector ID (defined in the steering file +of the \musrSim), the second variable is the name of the counter, +the third variable identifies the type of the counter (M = muon, P = positron, +K = coincidence and V = veto counters, respectively), the forth variable is the threshold in MeV applied +to the energy deposited in the counter, the fifth variable is a time delay (in units +of bin-width) applied to the given counter, the sixth set of parameters defines +which other detectors have to be in coincidence (positive ID) or anti-coincidence +(negative ID) with the given counter, and in the case of positron counters, the +last four parameters define a special histogram for the given counter. (Note +that this histogram is used for a compatibility with the experimental steering file, however +there is a different and more powerful way of how to define histograms, which is +defined later). The only difference with respect to the experimental steering +file is the presence of threshold definition in the fourth column. +In the example above, the M-counter is the volume which has ID 102 in the simulation, +and a valid hit has to have energy in the counter above 0.4\,MeV, and no signal +above 0.1\,MeV in the ``Active Veto'' (ID=401) within a time interval defined by +{\tt VCOINCIDENCEW}. The detailed description of *.v1190 steering files can be found +in chapter~\ref{description} and in~\cite{TDCsetup}. + + + +The main output of the simulation is generated in the form of histograms. +The histograms are defined in the steering file. For example, the following +line defines 1-dimensional histogram of the $z$-position of where the muons stop and decay:\\[1em] +{\tt +musrTH1D hMuDecayPosZ "Where the muons stop;z[mm];N" 100 -5. 5. muDecayPosZ +}\\[1em] +This histogram has 100 bins, spans from -5\,mm to 5\,mm, and the variable filled +in the histogram is {\tt muDecayPosZ}. +In fact, this line does not create one histogram, but an array of histograms -- one +histogram for each ``condition''. An example of five conditions reads:\\[1em] +{\tt +condition 1 oncePerEvent \\ +condition 2 muonDecayedInSample\_gen \\ +condition 4 muonTriggered\_det \\ +condition 6 goodEvent\_det \\ +condition 9 pileupEvent +}\\[1em] +where integer numbers (1,2,4,6 and 9) denote the condition number, and the string at the end +describes the condition -- e.g.\ ``{\tt muonDecayedInSample\_gen}'' specifies that muon +has stopped and decayed in the sample, and ``{\tt muonTriggered\_det}'' specifies that +there had to be a good muon candidate triggered in the M-counter. +The ending ``{\tt \_gen}'' indicates that the variable used in the condition +was ``generated'', i.e.\ it is known in the simulation, but can not be directly +measured experimentally. +On the other hand the ending ``{\tt \_det}'' specifies that the given condition +is based on measurable variables, as e.g.\ energy deposits in a counter. + +There can be up to +30 conditions requested, and for each of them a separate histogram will be filled. +In the output file of \musrSimAna, the histograms corresponding to the previous +set of conditions would be saved as +{\tt hMuDecayPosZ\_1, hMuDecayPosZ\_2, hMuDecayPosZ\_4, hMuDecayPosZ\_6, hMuDecayPosZ\_9}. +The {\tt hMuDecayPosZ\_1} shows where the muons stop irrespective whether they were +detected or not. The {\tt hMuDecayPosZ\_6} contains ``good'' muons, i.e.\ muons that would +be saved in the final histograms in the experiment. The {\tt hMuDecayPosZ\_9} shows +the muons, which contribute to the time-independent background. Note that there are always +two muons, $\mu_1$ and $\mu_2$, contributing to the time-independent background, and {\tt hMuDecayPosZ\_9} +shows the one ($\mu_1$) that started the M-counter. More to this topic is presented +in chapter~\ref{sec:eventMixing}. Here we just mention that in order +to see the $z$-coordinate of the second muon ($\mu_2$), +which was not triggered by the M-counter but whose decay positron hit the positron counter, +the histogram {\tt hpileup\_muDecayPosZ\_9} must be used:\\[1em] +{\tt +musrTH1D hpileup\_muDecayPosZ "Pileup muons;z[mm];N" 100 -5.0 5. pileup\_muDecayPosZ +} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Event mixing -- an unavoidable complication} +\label{sec:eventMixing} +The output of \musrSim\ is stored in a Root tree named ``{\tt t1}''. One event corresponds +to one simulated muon, and it is saved as one row of the tree. The only information that +relates one event (muon) to any other one is the variable {\tt timeToNextEvent}, +which is a randomly generated time difference between two subsequent events. + +The \musrSimAna\ allows one to study the ``time-independent background'', which is +due to mixing of two different events. A simple example of the event mixing affecting +the \musr\ measurement is the following: the first muon, $\mu_1$, hits the M-counter, +and subsequently stops and decays in the sample, however the decay positron escapes +undetected -- most likely because of the limited angular acceptance of positron counters. +The second muon, $\mu_2$, arrives around the same time, misses the M-counter, +and stops and decays in a collimator wall or elsewhere. +Its decay positron hits a positron counter. Thus there are good-looking muon and positron +hits, which however arise from two uncorrelated muons. +This fake (background) event is experimentally treated as a good event, +and contributes to the time-independent background. + +Events can be mixed in \musrSimAna, allowing us to study sources of +the time-independent background. +%\emph{MusrSimAna} loops over the events and checks whether the conditions defining +%the detector logic in the steering file have been fulfilled. However, to take event +%mixing properly into account, the conditions have to be veryfied on several subsequent +%events. Therefore, at the beginning of the analysis of +Before analysing a given event, arrays of hits are filled +for all counters (M, positron, veto, coincidence counters), which store the hits occurring in the +future up to the time equal to $ 3\, \cdot $ pileup window or $ 3\, \cdot $ data window, +whatever is larger\footnote{The data and pile-up windows defined by parameters +{\tt DATAWINDOWMIN, DATAWINDOWMAX, PILEUPWINDOWMIN} and {\tt PILEUPWINDOWMAX} +are applied later on in the analysis.}. +There is one such array for every counter. After this initial filling, +there might be several hits in every array, originating from one or more events. + + +The fraction of the time-independent background to good events depends on the incoming muon rate +measured by the trigger, possibly in the anti-coincidence with a backward veto detector, if +used in the experiment. +Typically, the experimentalists set the incoming muon rate (rate of \emph{stopped muons}) +to $\sim 30\,000\,\mu/s$. The same should be done in the simulation. +However, the rate of stopped muons is known only after the analysis is done, +because, for example, many simulated muons stop and decay in collimators or elsewhere +in the beam-pipe without producing any signal in the M-counter. +Therefore the simulation is normally started with an initial rate of generated muons +of $30\,000\,\mu/s$, which in practise can correspond to much lower rate of \emph{stopped muons}. +The rate of stopped muons is calculated at the end of the \musrSimAna\ execution, and +it is printed out in the \musrSimAna\ output. The user can use this information and rescale the +initial muon rate by changing parameter {\tt MUONRATEFACTOR} in the {\tt *.v1190} setup file. +This is done without the necessity to re-run the CPU consuming +simulation -- only the \musrSimAna\ has to be repeated. The complete simulation and analysis +chain therefore usually consists of three steps: +% +\begin{enumerate} + \item \musrSim\ (takes typically 10 hours). + \item \musrSimAna\ with {\tt MUONRATEFACTOR=1.0} (takes typically 10 minutes). + \item \musrSimAna\ with {\tt MUONRATEFACTOR} determined in step 2. +\end{enumerate} +% +The {\tt MUONRATEFACTOR} specifies a multiplicative factor applied to +the variable {\tt timeToNextEvent} +(the randomly generated time difference between two subsequent events). + +{\bf IMPORTANT NOTE:} In order to get the pile-up effects correctly analysed by \musrSimAna, +it is probably necessary to run the \musrSim\ with no event reweighting (i.e.\ +the command ``{\tt /musr/command logicalVolumeToBeReweighted \ldots}'' must {\bf not} be used). +All events should/have to be (?) saved in the Root tree +(i.e.\ the command ``{\tt /musr/command storeOnlyEventsWithHits false}'' must be used). + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Detailed list of steering file parameters} +\label{description} +\begin{description} +% \item{\bf INSTRUMENT=\emph{string}} \\ +% ignored +% \item{\bf DESCRIPTION=\emph{string}} \\ +% ignored +% \item{\bf TYPE=\emph{string}} \\ +% ignored + \item{\bf RESOLUTION=\emph{value}} \\ + width of the TDC bin in picoseconds. + \item{\bf MDELAY=\emph{value}} \\ + currently not used (probably not needed in the case of simulation). + \item{\bf PDELAY=\emph{value}} \\ + currently not used (probably not needed in the case of simulation). + \item{\bf MCOINCIDENCEW=\emph{value}} \\ + time window for the coincidences of the coincidence detectors (``K'') + with the M-counter. The \emph{value} is given in TDC bins (see {\tt RESOLUTION} above). + \item{\bf PCOINCIDENCEW=\emph{value}} \\ + time window for the coincidences of the coincidence detectors (``K'') + with positron counters. The \emph{value} is given in TDC bins. + \item{\bf VCOINCIDENCEW=\emph{value}} \\ + time window for the coincidences of the anti-coincidence detectors (``V'') + with any other counter. The \emph{value} is given in TDC bins. + \item{\bf MUONRATEFACTOR=\emph{value}} \\ + a multiplicative factor which is used to rescale time differences between subsequent muons. + Setting \emph{value} larger than 1 artificially prolongs the time difference between + 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 + 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''. + \item{\bf INFINITELYLOWMUONRATE} \\ + 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 + no pileup can be observed. The variable ``MUONRATEFACTOR'' becomes irrelevant when + 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}} \\ + Beginning of the data interval for the positron counters in $\mu$s. + \item{\bf DATAWINDOWMAX=\emph{value}} \\ + End of the data interval for the positron counters in $\mu$s. + \item{\bf PILEUPWINDOWMIN=\emph{value}} \\ + Beginning of the pileup interval for the M-counter in $\mu$s. + \item{\bf PILEUPWINDOWMAX=\emph{value}} \\ + End of the pileup interval for the M-counter in $\mu$s. + \item{\bf PROMPTPEAKMIN=\emph{value}} \\ + Beginning of the prompt-peak interval in $\mu$s. This variable is used only for the condition + ``{\tt promptPeak}'', ``{\tt promptPeakF}'', etc.\ , and normally does not need to be specified. It becomes useful if + the user wants to investigate, where the prompt-peak originates from (where do the muons, + which give rise to the prompt peak, stop). The default value is -0.01\,$\mu$s. + \item{\bf PROMPTPEAKMAX=\emph{value}} \\ + End of the prompt-peak interval in $\mu$s, the default value is 0.01\,$\mu$s. (See comments + for {\tt PROMPTPEAKMIN}.) + \item{\bf MUSRMODE=\emph{string}} \\ + Defines the mode of \musr\ experiment -- presently only ``D'', corresponding to + the time differential mode, and ``P'', for pulsed mode, are implemented. + \item{\bf REWINDTIMEBINS=\emph{value}} \\ + 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. + \item{\bf DEBUGEVENT \emph{eventID} \emph{debugLevel}}\\ + Prints out debug information about the event with the ID \emph{``eventID''}. + The higher the \emph{debugLevel}, the more details are printed. + (Both \emph{eventID} and \emph{debugLevel} are integers). + \item{\bf CLONECHANNEL \emph{ID1} \emph{ID2}}\\ + Clones the hits detected in counter \emph{ID1} into a new counter \emph{ID2}. + A different (even a lower) threshold can be applied to the cloned counter. + This way the same counter can be used as two independent counters -- e.g.\ once as a veto + detector for the M-counter, and simultaneously as a coincidence detector for + a P-counter. In both cases the energy threshold and time windows are defined independently. + \item{\bf WRITE\_OUT\_DUMP\_FILE \emph{fileNameString} \emph{clockChannel} \emph{randomJitter}}\\ + If present, this command will create two output files, the so-called ``dump files'':\\ + {\tt data/TDC\_V1190\_200\_dump\_\emph{fileNameString}.bin} -- file that can be used + as an input to the PSI analysis front-end of a real experiment.\\ + {\tt data/TDC\_V1190\_200\_dump\_\emph{fileNameString}.txt} -- file that contains the same + information (hits) as the previous file, however in a human-readable form. The first number in the file + stands for the channel number, the second number stands for the time bin in the TDC bin units.\\ + \emph{clockChannel} ... the channel of the clock signal (typically 15).\\ + \emph{randomJitter} ... this value is in TDC bins, typically 0 or 8000 (?). If \emph{randomJitter} is smaller then + 1, then the hits in the dump files will be sorted according to time. If it is larger than 0, then + subsequent hits can be unordered in time, but the time difference never exceeds the value of \emph{randomJitter}. + This is just a technical thing serving to test the analysis software -- it should not + have any effect on the analysis results. + \item{\bf musrTH1D \emph{histoName} \emph{histoTitle} \emph{nBins} \emph{min} \emph{max} \emph{variable} + [{\tt rotreference} $\nu_{\rm RRF}$ $\phi_{\rm RRF}$] $|$ [correctexpdecay]} \\ + Defines a histogram (or more precisely an array of histograms, where the number of histograms + in the array is given by the number of conditions, see section~\ref{howToAnalyse}). + The name of the histogram is defined by \emph{histoName} + the number of the condition. + The string variable \emph{histoTitle} specifies the title of the histogram, + \emph{nBins}, \emph{min} and \emph{max} stand for the number of bins and minimum and maximum + of the $x$-axis of the histogram. \\ + The optional keyword ``{\tt rotreference}'' signals that the given histogram will be filled in + rotating reference frame (RRF) with the frequency of $\nu_{\rm RRF}$ and a phase shift of $\phi_{\rm RRF}$. + \\ + The optional keyword ``{\tt correctexpdecay}'' signals that the given histogram will be corrected + for the muon exponential decay (i.e. multiplied by a factor $\exp(t/2.19703)$. It is meaningful + only for time variables. + \\ + The \emph{variable} stands for the variable that will be + filled into the histogram. The \emph{variable} can be any variable from the output Root tree + of musrSim (see ``Manual of musrSim'') (except for the array variables like + {\tt det\_*[], save*[], odet\_*[]}). In addition, it can be also one of the following: + \begin{description} + \item[muDecayPosR] \ldots $\sqrt{ {\rm muDecayPosX}^2 + {\rm muDecayPosY}^2}$. + \item[wght] \ldots weight of the event. + \item[det\_m0edep] \ldots energy deposited in the M-counter that gives the muon signal. + \item[det\_posEdep] \ldots energy deposited in the P-counter that gives the positron signal. + \item[muIniPosR] \ldots $\sqrt{ {\rm muIniPosX}^2 + {\rm muIniPosY}^2}$. + \item[muIniMomTrans] \ldots $\sqrt{ {\rm muIniMomX}^2 + {\rm muIniMomY}^2}$. + \item[muTargetPol\_Theta] \ldots theta angle of the muon spin when muon enters target (-180,180 deg). + \item[muTargetPol\_Theta360]\ldots theta angle of the muon spin when muon enters target (0,360 deg). + \item[muTargetPol\_Phi] \ldots phi angle of the muon spin when muon enters target (-180,180 deg). + \item[muTargetPol\_Phi360] \ldots phi angle of the muon spin when muon enters target (0,360 deg). + \item[pos\_Momentum] \ldots magnitude of the momentum of the decay positron (``generated'', not measurable variable). + \item[pos\_Trans\_Momentum] \ldots transverse momentum of the decay positron. + \item[pos\_Radius] \ldots positron radius calculated from the decay positron momentum and magnetic + field at the point of decay. + \item[pos\_Theta] \ldots polar angle of the decay positron. + \item[pos\_Phi] \ldots azimuth angle of the decay positron. + \item[det\_time10] \ldots time difference between the positron and muon counters + (measured by the respective counters). + \item[gen\_time10] \ldots the time difference between the muon decay and the first muon hit + in the M-counter (i.e.\ {\tt muDecayTime - muM0Time}). + \item[det\_time10\_MINUS\_gen\_time10] \ldots {\tt det\_time10 - gen\_time10} in picoseconds. + \item[det\_time20] \ldots very similar to {\tt det\_time10}, however taking into + account ``counter phase shift'' defined by {\tt counterPhaseShifts} + variable. This gives the user a possibility to sum up backward and + forward histograms into one histogram (this of course make sense + only in the simulation, where there is ``no physics'' happening + in the sample, just the muon spin rotation). + \item[pileup\_eventID] \ldots eventID of the $\mu_2$, where $\mu_2$ stands for the muon, + which did not give signal in the M-counter, but whose + decay positron gave signal in the positron counter around the time + when a different muon ($\mu_1$) triggered the M-counter. + \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\_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} + 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). + \item{\bf musrTH2D \emph{histoName} \emph{histoTitle} \emph{nBins} \emph{min} \emph{max} \emph{nBins2} \emph{min2} \emph{max2} \emph{variable}} \\ + Similar to \emph{musrTH1D}, but for a 2-D histogram. + \item{\bf humanDecayHistograms \emph{hist\_decay\_detID} \emph{hist\_decay\_detID\_pileup} \emph{id$_1$} \emph{name$_1$} \ldots \emph{id$_n$} \emph{name$_n$} } \\ + This is a special kind of histogram, which converts two histograms + (\emph{hist\_decay\_detID} \emph{hist\_decay\_detID\_pileup}) + into a human-friendly histograms, where detector ID on the $x$-axis is converted into a string label. + The \emph{id$_i$} is the detector id, and the \emph{name$_i$} is the corresponding label name. + If \emph{name$_i$ = name$_j$}, the corresponding bins of the original histograms will + be summed up together into the same bin. + The \emph{hist\_decay\_detID} and \emph{hist\_decay\_detID\_pileup} have to be defined before + (by the command {\tt musrTH1D}). + \item{\bf condition \emph{conditionID} \emph{conditionName}} \\ + Definition of a condition, which is then used when filling histograms. The \emph{conditionID} + specifies the number of the condition, and must be between 0 and 30 (0 and 30 are also possible). + The \emph{conditionName} is one of the following: + \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[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 (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) + -- not a measurable variable. + \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. + \item[positronHit\_det] \ldots true if a good positron candidate was found in the positron counter. + Double hits within the data window are excluded. + \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 + (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. + 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[pileupEventCandidate] \ldots M-counter hit and positron counter hit both come from two different events. + \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\_F\_det] \ldots {\tt goodEvent\_det}, where the positron was detected in the forward detectors + defined by the command {\tt counterGrouping}. + \item[goodEvent\_B\_det] \ldots like {\tt goodEvent\_F\_det} but for backward positron counters. + \item[goodEvent\_U\_det] \ldots like {\tt goodEvent\_F\_det} but for upper 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\_R\_det] \ldots like {\tt goodEvent\_F\_det} but for right positron counters. + \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\_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\_L\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_L\_det} and {\tt pileupEvent}. + \item[goodEvent\_R\_det\_AND\_pileupEvent] \ldots {\tt goodEvent\_R\_det} and {\tt pileupEvent}. + \item[promptPeak] \ldots {\tt goodEvent\_det}, and {\tt PROMPTPEAKMIN < det\_time10 < PROMPTPEAKMAX}. + \item[promptPeakF] \ldots like {\tt goodEvent\_F\_det} and {\tt promptPeak}. + \item[promptPeakB] \ldots like {\tt goodEvent\_B\_det} and {\tt promptPeak}. + \item[promptPeakU] \ldots like {\tt goodEvent\_U\_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[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} + Additional conditions may be implemented on request. + \item{\bf draw \emph{histogramName} \emph{conditionID} } \\ + Plot histogram (for a given condition) at the end of the analysis. Note that all histograms + are saved into the output file irrespective whether they are plotted or not. + \item{\bf counterPhaseShifts \emph{ID$_1$} \emph{$\phi_1$} \emph{ID$_2$} \emph{$\phi_2$} \ldots \emph{ID$_n$} \emph{$\phi_n$} }\\ + Defines relative phase shifts of signal between different positron counters, which is used + for calculation variable {\tt det\_time20}. \emph{ID$_i$} is the ID number of the positron counter, + \emph{$\phi_i$} is its phase shift. This gives the user a possibility to sum up backward and + forward histograms into one histogram. + \item{\bf counterGrouping \emph{group} \emph{ID$_1$ ID$_2$ \ldots ID$_n$} } \\ + This defines a group of detectors, where \emph{group} stands for ``B'' (backward), + ``F'' (forward), ``U'' (up), ``D'' (down), ``L'' (left) and ``R'' (right) detectors. + This grouping is used in the definition of some conditions. + \item{\bf sampleID \emph{ID$_1$ ID$_2$ \ldots ID$_n$} } \\ + Defines which volume (or volumes, if there are more) is the sample. Typically, the + sample is just one volume, but sometimes there is a smaller volume inside of the sample, + to which a field is applied, so the small volume also has to be considered as the sample. + This information is needed for the condition ``muonDecayedInSample\_gen''. + \item{\bf setSpecialAnticoincidenceTimeWindow \emph{detectorID} \emph{timeMin} \emph{timeMax} \emph{unit}} \\ + This command sets a special anti-coincidence time window for a detector \emph{detectorID}. + Normally, the anti-coincidence time window is defined by {\tt VCOINCIDENCEW}, and is the same for all anti-coincidence + detectors. However, sometimes it might be interesting to set the anti-coincidence time window + differently for a specific detector (e.g.\ one might test an anti-coincidence of a veto detector with + the M-counter for the whole pile-up time window of $\sim$\,10\,$\mu s$. + Unlike in the case of {\tt VCOINCIDENCEW}, here the \emph{units} are not TDC bins, but + rather time in ``nanosecond'' or ``microsecond''. + \item{\bf fit \emph{histogramName} \emph{function} \emph{option} \emph{min} \emph{max} \emph{p$_1$} \ldots \emph{p$_n$}} \\ + Fits the histogram by a given function, where \emph{min}, \emph{max} define the range of the fit + on the $x$-axis of the histogram, \emph{option} is a string defining fit options (see Root manual for details), + and \emph{p$_1$} \ldots \emph{p$_n$} are (typically, with some exceptions) + the initial values of the function parameters. The following functions are currently predefined: + \begin{description} + \item[pol0] $=p_0$ \ldots a constant (1 parameter) - typically used to fit background. + \item[simpleExpoPLUSconst] $=p_0 \exp(-x/2.19703)+p_1$ + \item[rotFrameTime20] $= p_2 \cos(p_0 x+p_1)$ + \item[funct1] $=p_3 \exp((p_4 - x)/2.19703) \cdot (1+p_2 \cos(p_0 x+p_1))$ + \item[funct2] $=p_3 \exp((p_4 - x)/2.19703) \cdot (1+p_2 \cos(p_0 x+p_1)) + p_5$ + % \item[funct3] the same as {\tt funct2} + \item[funct4] $=p_3 \exp((- x)/2.19703) \cdot (1+p_2 \cos(p_0 x+p_1)) + p_4$ + \item[TFieldCos] $=p_3 (1+p_2 \cos(p_0 x + p_1))$ \hspace{1cm} (this function is useful when the histogram is filled with {\tt ``correctexpdecay''} keyword.) + \item[TFieldCosPLUSbg] $=p_3 (1+p_2 \cos(p_0 x + p_1)) + p_4 \exp(x/2.19703)$ + \hspace{1cm} (this function is useful when the histogram is filled with {\tt ``correctexpdecay''} keyword.) + \item[gaus] ... Gauss distribution + \end{description} + + +\end{description} +%======================================================================================================== +\section{A real-life example: GPD instrument} +\label{sec:GPD} +The simulation of the General Purpose Decay-Channel Spectrometer (GPD) +instrument~\cite{GPD} installed at PSI has been exemplified +in the \musrSim\ manual~\cite{musrSim}. +Here we analyse the output of this simulation using \musrSimAna. +The run number of this simulation is 201, therefore the steering +file names are ``{\tt 201.mac}'' for \musrSim, and ``{\tt 201.v1190}'' for +\musrSimAna, respectively, and the output file name of \musrSim\ is saved as +``{\tt data/musr\_201.root}''. +The detector system is very simple with only six counters -- M-counter, +two backward positron counters and three forward positron counters. +The reader is strongly recommended to see the illustration of the GPD +geometry in the \musrSim\ manual~\cite{musrSim}. + +8\,000\,000 of events were simulated +(i.e.\ 8\,000\,000 of muons were generated 100\,cm in front of +the GPD sample). +In only 949\,759 events (11.9\% out of the 8 million) there was a signal detected +in one or more counters. The remaining muons stopped somewhere (most +often in collimator, as we will see later), decayed, and the decay positron +(and any other particles created thereof) missed the counters. +This is illustrated in more details in Fig.~\ref{det_n}, +% +\begin{figure}[tbp]\centering +\epsfig{file=pict/det_n_infititely_low_muon_rate.eps,width=0.8\linewidth,clip=} +\caption{Number of hits in all counters per event, assuming infinitely low incoming muon +rate. The same detector may be hit more than once (e.g.\ if both the muon and its decay +positron pass through the M-counter).} +\label{det_n} +\end{figure} +% +where number of detector hit per event, assuming infinitely low incoming muon +rate, is shown. This plot was created in Root by executing:\\[1em] +{\tt +root [0] TFile* ff=new TFile("data/musr\_201.root") \\ +root [1] t1->Print() \\ +root [2] t1->Print("det\_n","det\_n>0") \\ +}\\[1em] +% +It has to be pointed out, that the ratio of muons passing through the opening +in collimators to the number of all generated muons strongly depends on the +beam properties -- beam profile, beam convergence, etc. Typically, if we have +too broad muon beam, we simulate +many ``useless'' events. However, the other extreme (simulating too narrow +beam) can lead to underestimating the time-independent background. + +It took approximately 12 hours of the CPU (on PC bought in 2010, where 1 out +of 4 processor cores was running) to simulate these 8\,000\,000 events. +Assuming the 30\,000\,$\mu/$s trigger rate, this corresponds to 26 seconds +of real experimental running. + +\subsection{Where the muons stop and decay} +\label{sect_muons} +The positions, or more precisely the components of the GPD instrument, where the muons +stop and decay, are shown in Fig.~\ref{humanDecayHistograms_1}: +% +\begin{figure}[tbp]\centering +\epsfig{file=pict/Plot201_1.eps,width=0.9\linewidth,% +%%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt,clip=} +clip=} +\caption{This plot indicates, where the muons stopped and decayed. +The dashed histogram shows all generated muons. The full-line histograms show +where stopped the muons, for which either the muon itself or its secondary +particle ($e^+, \gamma$) triggered the M-counter: black histogram stands for +all such muons, corresponding to infinitely low incoming muon rate, while +the red histogram stands for the incoming muon rate of 30\,000\,$\mu/$s. +8\,000\,000 of events were simulated.} +\label{humanDecayHistograms_1} +\end{figure} +% +%Notes to Fig.~\ref{humanDecayHistograms_1}: +\begin{itemize} + \item Figure~\ref{humanDecayHistograms_1} was generated by Root macro file ``Plot201.C''. + \item The labels on the $x$-axis are defined in the file {\tt 201.v1190} by the + command \\ + {\tt humanDecayHistograms \ldots} + \item The dashed-line histogram in Fig.~\ref{humanDecayHistograms_1} + shows where the muons stopped and decayed if no preselection + criteria are applied on the muons, i.e.\ if all generated muons are considered. + This is histogram ``{\tt humanDecayHistograms\_1}''. + \item The full-line histograms show + where stopped the muons, for which either the muon itself or its secondary + particle ($e^+, \gamma$) triggered the M-counter: black histogram stands for + all such muons, corresponding to infinitely low incoming muon rate, while + the red histogram represents the case for the 30\,000\,$\mu/$s incoming muon rate. + An energy deposit of at least 0.4\,MeV in the M-counter is required to fire the trigger. + The number of triggered events decreases with the incoming muon rate, + because some of the events are rejected due to the 10\,$\mu$s pileup gate. + + The histogram name is in both cases ``{\tt humanDecayHistograms\_4}'', + where the black histogram was calculated using the setup file ``{\tt 201a.v1190}'' + with the keyword {\tt INFINITELYLOWMUONRATE}, while the red histogram + was calculated using the setup file ``{\tt 201.v1190}'' + with {\tt MUONRATEFACTOR=0.0965819}. + \item The $\pm 10\,\mu$s pile-up gate at the incoming muon rate of 30\,000$\,\mu/$s + rejects approx.\ 45\% of the triggered events. This number can be calulated + in Root as the ratio of the ``{\tt Integral}'' of the red and black histograms + in Fig.~\ref{humanDecayHistograms_1}:\\[1em] + {\tt \small + root [0] TFile* file1 = new TFile("data/his\_201\_201a.v1190.root") \\ + root [1] humanDecayHistograms\_4->Integral() \\ + root [0] TFile* file2 = new TFile("data/his\_201\_201.v1190.root") \\ + root [1] humanDecayHistograms\_4->Integral() \\ + } + \item The muon sample fraction (ratio of muons stopped in the sample over all muons that fired + the trigger) for the triggered events (full-line histograms) + is 65\%, and it is practically the same + for both infinitely low and 30\,000\,$\mu/$s incoming rate. + This number can be obtained in Root by dividing the first column of histogram + {\tt humanDecayHistograms\_4} by the sum of all entries in this histogram:\\[1em] + {\tt \small + root [0] TFile* file = new TFile("data/his\_201\_201.v1190.root") \\ + root [1] (humanDecayHistograms\_4->GetBinContent(1))/(humanDecayHistograms\_4->Integral()) \\ + } + \item The largest fraction of generated muons (dashed-line histogram) stopped in collimators. + Only a small fraction of them caused a hit in the M-counter (full-line histograms). + \item Despite the high initial muon momentum of $100 \pm 3\,$GeV/c, muons are + significantly scattered in the last 50\,cm region of air. This can be + clearly seen if the magnetic field is off and a point-like muon beam + is used (which can be done by modifying the {\tt 201.mac} file) + -- only 77\% of the muons stop in the sample cell or in the sample, while the + remaining 23\% of the mouns are scattered so much in the air, that they + end up in collimators or elsewhere (not shown here). + \item ``World'' in the histogram label means that the muon decayed in the beampipe vacuum + or somewhere else in the air (on the fly). + \item ``Escaped'' means that the muon left the simulated instrument (more precisely the + ``world'' volume) prior its decay. +\end{itemize} + + +Figure~\ref{humanDecayHistograms_9} shows the ``pile-up events''. +% +\begin{figure}[htbp]\centering +\epsfig{file=pict/Plot201_2_new.eps,width=0.9\linewidth,% +%%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt,clip=} +clip=} +\caption{Pile-up events, i.e.\ the events in which one muon ($\mu_1$) fired the +trigger, while the hit in a positron counter is due to a decay positron from +a different muon ($\mu_2$). Pile-up events look like a good events, and contribute +to the time-independent background.} +\label{humanDecayHistograms_9} +\end{figure} +% +These are events, in which one muon ($\mu_1$) is triggered by the m-counter, +while a positron from a different muon ($\mu_2$) was detected by +a positron counter\footnote{In fact, the trigger may also be triggered by +the decay positron of $\mu_1$ and/or a positron counter may detect +directly $\mu_2$, not its decay positron. Such cases are rare, but they +are implicitly included in the simulation.}. +In addition to this requirement, the decay positron of $\mu_1$ must +escape undetected (e.g.\ it must miss positron counters) and $\mu_2$ must not trigger the m-counter +-- otherwise the event would be rejected. +Pile-up events are the source of the time independent background. +Usually $\mu_1$ is a good-looking muon that stops in the sample or in the sample cell +(red histogram in Fig.~\ref{humanDecayHistograms_9}), while $\mu_2$ stops and decays at different places, +mainly in the collimators (green histogram in Fig.~\ref{humanDecayHistograms_9}). + +A nice visualisation of where the background-contributing muons $\mu_2$ stop and decay +is presented in Fig.~\ref{Pileup_muon_decay_map} (histogram ``{\tt hMuDecayMappileup\_9}''). +% +\begin{figure}[htbp]\centering +\epsfig{file=pict/Pileup_muon_decay_map.eps,width=0.7\linewidth,% +%%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt,clip=} +clip=} +\caption{Positions of where the $\mu_2$ stop and decay.} +\label{Pileup_muon_decay_map} +\end{figure} +% +In this two dimensional histogram, different components of the GPD instrument, +like the lead collimator, the copper collimator and the sample cell, can be recognised. +The lead collimator is located at the $z$-position between -115\,mm and -85\,mm. +Due to the high initial muon momentum of $\sim 100\,$MeV/c, +the maximum of muons in Fig.~\ref{Pileup_muon_decay_map} stop quite deep in the +lead collimator, at around $z=-103$\,mm. This might be a little bit surprising to the +\musr\ scientists who are used to work with the surface muons with momentum of 28\,MeV/c. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{The $\mu$SR signal} +% +Figure~\ref{hdet_time10_10} +% +\begin{figure}[htbp]\centering +\epsfig{file=pict/hdet_time10_10.eps,width=0.495\linewidth,% +bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} +\epsfig{file=pict/hdet_time10_11.eps,width=0.495\linewidth,% +bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} +%%clip=} +\caption{MuSR signal for the run 201 (TF$=300\,$gauss). The tree forward positron counters +are summed up in the left histogram, and the two backward counters in the right histogram.} +\label{hdet_time10_10} +\end{figure} +% +shows the $\mu$SR spectra for the same run, +i.e.\ for the transverse field of 300\,gauss, integrated over the three forward positron +counters (left histogram called {\tt hdet\_time10\_10}) +and over the two backward positron counters (right histogram called {\tt hdet\_time10\_11}). +Zero on the time axis corresponds to $t_0$, i.e.\ time of the m-counter hit. +One can see a prompt peak at $t_0$, time independent background at negative times +and an oscillating signal at positive times. +The following function has been fitted to the oscillating part of the signal: +% +\begin{equation} +f=p_3 \cdot e^{-t/2.19703} \cdot (1+p_2 \cdot \cos(t \cdot p_0+p_1))+p_4 +\label{eq_simple} +\end{equation} +The fits were restricted to the time interval of $(t_0+0.05 \mu\rm{s},t_0+9.95\mu\rm{s})$, +and the parameter $p_0$ was fixed (e.g. not fitted). +The fitted amplitude of asymmetry are $p_2 = 0.307 \pm 0.009$ and +$p_2 = 0.290 \pm 0.009$ for the forward and backward counters respectively. + +Parts of the spectra from Fig.~\ref{hdet_time10_10} are shown +in detail in Fig.~\ref{hdet_time10_10_detail}. +% +\begin{figure}[htbp]\centering +\epsfig{file=pict/hdet_time10_10_detail.eps,width=0.495\linewidth,% +bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} +\epsfig{file=pict/hdet_time10_11_pileup.eps,width=0.495\linewidth,% +bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} +%%clip=} +\caption{MuSR signal for the run 201 (TF$=300\,$gauss) -- details of +Fig.~\ref{hdet_time10_10}. The left plot shows the signal in the forward counters around $t_0$, +the right plot shows the (time-independent background) signal at negative times in the +backward counters.} +\label{hdet_time10_10_detail} +\end{figure} +% +The left plot in Fig.~\ref{hdet_time10_10_detail} shows the signal +in the forward counters around $t_0$, the right plot shows the +(time-independent background) signal at negative times in the backward counters. + +An important characteristic of a \musr\ instrument is the time-independent +background. It is usually expressed as +% +\begin{equation} +{\rm Bgr} = p_{-} / p_3 ~~~, +\label{eq_background} +\end{equation} +% +where $p_{-}$ is the fit to the time-independent background, i.e.\ signal at negative times, +and $p_3$ is the parameter from eq.(\ref{eq_simple}), which specifies what the +size of the signal would be at $t_0$ in the absence of oscillations. +In the case of backward counters ${\rm Bgr}_{\rm backw} = 14.47/262 = 5.5\,\%$, +in the case of forward counters ${\rm Bgr}_{\rm forw} = 6.88/267.9 = 2.6\,\%$. + +Note that the histogram on right hand side of Fig.~\ref{hdet_time10_10_detail} +is labelled ``{\tt hdet\_time10\_Bgr\_11}'', not ``{\tt hdet\_time10\_11}''. +In fact, the two histograms are identical, as one can see in the setup file +{\tt 201.v1190}. The only difference is in the fitting -- the same data stored in +both histograms are fitted by different functions in different time ranges. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{The $\mu$SR signal from individual counters} +% +Figure~\ref{F11} shows the observed signal in the +forward counter\ No.~11 (FW11). +% +\begin{figure}[htbp]\centering +\epsfig{file=pict/F11_rebinned.eps,width=0.495\linewidth,% +bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} +\epsfig{file=pict/F11_B11_prompt_peak_thicker.eps,width=0.495\linewidth,% +bbllx=13pt,bblly=5pt,bburx=520pt,bbury=351pt,clip=} +\caption{\musr\ signal in the forward positron counter\ No.~11 (run 201, TF$=300\,$gauss). +The left plot shows the (rebinned) signal in the counter, +the right plot shows the detail of the \emph{prompt peak}, i.e.\ the region +around $t_0$ in the same counter (black line), +compared with the prompt peak in the backward positron counter\ No.~1 (magenta line).} +\label{F11} +\end{figure} +% +Originally, the histogram F11 was defined with the bin width of 100\,ps. +The number of bins was 50995, covering the time interval of approx.\ (-0.2\,$\mu$s, 4.9\,$\mu$s). +In the left hand side plot, however, the histogram was rebinned +(200 bins were summed up into 1 bin). +The right hand side plot shows the detail of the \emph{prompt peak}, i.e.\ the region +around $t_0$, of one forward and one backward positron counters, prior to the rebinning. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{Conclusion of the GPD analysis example} +% +The purpose of the example analysis of the GPD simulation was to illustrate +the potential of \musrSim\ and \musrSimAna\ programs to investigate features +like time-independent background, sample muon fraction, prompt peak, \ldots +This information can be used in design and optimisation of \musr\ instruments. +%======================================================================================================== +\section{GPS instrument} +% +It is foreseen that GPS instrument could be arranged in two geometries +after the upgrade (depending from which side the calorimeter would be +inserted). +\begin{itemize} + \item {\tt 50130hb.v1190} -- Calorimeter inserted from one side. + \item {\tt 50130hl.v1190} -- Calorimeter inserted from the other side. + \item {\tt 50130hb1.v1190 -- 50130hb6.v1190} -- All positron counters + analysed individually. +\end{itemize} +See the document about the GPS simulations saved in the directory: \\ +/afs/psi.ch/project/lmu/Facility/musr\_simulations/documentation/GPS/ \\ +for more details. +%======================================================================================================== +\section{Other Examples} +Many different ``*.v1190'' files are stored in the file: +``run\_musrSimAna\_many\_files.tar.gz''. They could serve as additional examples. +Note that the syntax of the ``fit'' command was changed +at some point, and therefore the ``fit'' command might cause problems +(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} + +\bibitem{acquisitionProkscha} T.~Prokscha {\it et al.} ``A novel VME based \musr\ data acquisition system at PSI'', +Physica {\bf B~404}, (2009) 1007-1009. + +\bibitem{TDCsetup} ``TDC Manual -- Setting up the required logic'', +http://lmu.web.psi.ch/facilities/electronics/TDC/set\_logic.html + +\bibitem{GPD} +http://lmu.web.psi.ch/facilities/gpd/gpd.html + +\bibitem{musrSim} +K.Sedlak {\it et al.}, ``Manual of musrSim''. + + +\end{thebibliography} + +\end{document} diff --git a/musrSim.cc b/musrSim.cc index 55e4853..c09c16c 100644 --- a/musrSim.cc +++ b/musrSim.cc @@ -10,7 +10,9 @@ #include "G4RunManager.hh" #include "G4UImanager.hh" #include "G4UIterminal.hh" -#include "G4UItcsh.hh" +#ifdef G4UI_USE_TCSH +#include "G4UItcsh.hh" //JSL: should be commented on windows ? +#endif //#include //#include @@ -24,7 +26,7 @@ #include "Randomize.hh" -#include +// #include //JSL #ifdef G4VIS_USE // #include "musrVisManager.hh" diff --git a/musrSimAna/musrAnalysis.cxx b/musrSimAna/musrAnalysis.cxx index fa70903..49958ce 100644 --- a/musrSimAna/musrAnalysis.cxx +++ b/musrSimAna/musrAnalysis.cxx @@ -1,1522 +1,2094 @@ -#define musrAnalysis_cxx -#include "musrAnalysis.hh" -#include "musrTH.hh" -//#include -#include -#include -#include -#include -//#include - - -musrWriteDump* musrAnalysis::myWriteDump = NULL; - -void musrAnalysis::Loop(char* runChar, char* v1190FileName, Int_t nrEvents) -{ -// In a ROOT session, you can do: -// Root > .L musrAnalysis.C -// Root > musrAnalysis t -// Root > t.GetEntry(12); // Fill t data members with entry number 12 -// Root > t.Show(); // Show values of entry 12 -// Root > t.Show(16); // Read and show values of entry 16 -// Root > t.Loop(); // Loop on all entries -// - -// This is the loop skeleton where: -// jentry is the global entry number in the chain -// ientry is the entry number in the current Tree -// Note that the argument to GetEntry must be: -// jentry for TChain::GetEntry -// ientry for TTree::GetEntry and TBranch::GetEntry -// -// To read only selected branches, Insert statements like: -// METHOD1: -// fChain->SetBranchStatus("*",0); // disable all branches -// fChain->SetBranchStatus("branchname",1); // activate branchname -// METHOD2: replace line -// fChain->GetEntry(jentry); //read all branches -//by b_branchname->GetEntry(ientry); //read only this branch - if (fChain == 0) return; - - std::cerr<<"musrSimAnalysis::Loop() : Analysing run "<GetEntriesFast(); - if (nrEvents!=0) nentries = nrEvents; - - Long64_t nbytes = 0, nb = 0; - for (Long64_t jentry=0; jentryGetEntry(jentry); InitialiseEvent(); nbytes += nb; - // if (Cut(ientry) < 0) continue; - if (jentry%100000==0) std::cout << "Analysing event nr. jentry=" << jentry << std::endl; - AnalyseEvent(jentry); - } - SaveHistograms(runChar,v1190FileName); -} - - - -//================================================================ - -void musrAnalysis::ReadInInputParameters(char* charV1190FileName) { - std::cout<<"musrAnalysis::ReadInInputParameters()"< tmpCoincidenceMultimap; - // for (int i=0; i("a", 1)); - // - // } - // } - while (!feof(fSteeringFile)) { - fgets(line,10000,fSteeringFile); - if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')&&(line[0]!='$')&&(line[0]!='!')) { - char tmpString0[200]="Unset"; - char tmpString1[200]="Unset"; - sscanf(&line[0],"%s",tmpString0); - if (strncmp(tmpString0,"INSTRUMENT=",strlen("INSTRUMENT="))==0) { - instrument=&line[strlen("INSTRUMENT=")]; - } - else if (strncmp(tmpString0,"DESCRIPTION=",strlen("DESCRIPTION="))==0) { - description=&line[strlen("DESCRIPTION=")]; - } - else if (strncmp(tmpString0,"TYPE=",strlen("TYPE="))==0) { - tdchwtype=&line[strlen("TYPE=")]; - } - else if (strncmp(tmpString0,"RESOLUTION=",strlen("RESOLUTION="))==0) { - sscanf(&line[strlen("RESOLUTION=")],"%g",&tmpfvar); - tdcresolution = tmpfvar * picosecond; - if (tdcresolution<1*picosecond) { - std::cout<<"\n !!! ERROR !!! Time resolution of TDC must be larger than 1 ps (you requested "< too low time resolution might lead to numerical problems in musrCounter ==> S T O P"<(ichannel_orig_tmp,ichannel_new_tmp)); - } - else if (strcmp(tmpString0,"DEBUGEVENT")==0) { - int ieventToDebug_tmp, iLevelToDebug_tmp; - sscanf(&line[0],"%*s %d %d",&ieventToDebug_tmp,&iLevelToDebug_tmp); - bool_debugingRequired=true; - if (ieventToDebug_tmp>-1) { - debugEventMap.insert(std::pair(ieventToDebug_tmp,iLevelToDebug_tmp)); - } - else { - for (int j=0; j<=-ieventToDebug_tmp;j++) { - debugEventMap.insert(std::pair(j,iLevelToDebug_tmp)); - } - } - } - else if (strcmp(tmpString0,"WRITE_OUT_DUMP_FILE")==0) { - int clock_channelID; - int max_timeBin_jitter; - // int clock_interval_tmp; - // sscanf(&line[0],"%*s %s %d %d",tmpString1,&clock_channelID_tmp,&clock_interval_tmp); - sscanf(&line[0],"%*s %s %d %d",tmpString1,&clock_channelID,&max_timeBin_jitter); - // clock_channelID = clock_channelID_tmp; - // clock_interval = clock_interval_tmp; - musrCounter::bool_WriteDataToDumpFile = true; - // std::cout<<"tmpString1="< STOP !!!"< S T O P !!!"< ignored, histogram not defined! S T O P !!!"< ignored, histogram not defined! S T O P !!!"<second; - if (strcmp(tmpString0,"musrTH1D")==0) { - // std::cout<<"pointerToVariable1="<SetRotatingReferenceFrame(rot_ref_frequency, rot_ref_phase); - hInfo->Fill(21, rot_ref_frequency); // value may be overwritten - just the last rot. ref. frame will be saved - hInfo->Fill(22, rot_ref_phase); // value may be overwritten - just the last rot. ref. frame will be saved - } - else if (strcmp(furtherOption,"correctexpdecay")==0) { - myTH->SetExpDecayCorrection(); - } - } - else { - Double_t* pointerToVariable2 = iter2->second; - musrTH* myTH = new musrTH("2D",histoName,histoTitle,nrBinsX,minX,maxX,nrBinsY,minY,maxY,pointerToVariable1,pointerToVariable2); - listOfAllHistograms2D.push_back(myTH); - } - } - else if (strncmp(tmpString0,"humanDecayHistograms",strlen("humanDecayHistograms"))==0) { - char motherHistoName[200]; - char motherHistoPileupName[200]; - sscanf(&line[0],"%*s %s %s",motherHistoName,motherHistoPileupName); - std::cout<<" Defining humanDecayHistograms (motherHistoName="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { - if ((*it)->IsThisHistoNamed(motherHistoName)) { - motherOfHumanDecayHistograms = *it; - } - if ((*it)->IsThisHistoNamed(motherHistoPileupName)) { - motherOfHumanDecayPileupHistograms = *it; - } - } - if (motherOfHumanDecayHistograms==NULL) { - std::cout<<"\n\nmusrAnalysis::ReadInInputParameters Mother histogram (motherOfHumanDecayHistograms)" - <<"of the \"Human decay histogram\" not found!"< S T O P "< S T O P "<second)==(std::string) NAME1) { - nameDefined = true; - iBinHuman = it->first; - } - } - if (!nameDefined) { - nBinHuman++; - iBinHuman=nBinHuman; - humanDecayMap.insert(std::pair(iBinHuman,NAME1)); - } - humanDecayMultimap.insert(std::pair(iBinHuman,N1)); - char* pch2 = strstr(pch,NAME1)+strlen(NAME1); - pch=pch2; - nscan = sscanf(pch,"%d %s",&N1,NAME1); - // indexHuman++; - } while (nscan==2); - - Int_t nrBinsHumanHist=humanDecayMap.size(); - // Int_t nrBinsHumanHist=0; - // std::string oldString="blaBlaoldStringUndefiNED"; - // for (humanDecayMultimapType::const_iterator it = humanDecayMultimap.begin(); it!=humanDecayMultimap.end(); ++it) { - // if (oldString!=(it->first)) nrBinsHumanHist++; - // } - std::cout<<"\nnrBinsHumanHist="<SetBinLabel1D(k+1,sHumanDecayArray[k]); - // } - // oldString="blaBlaoldStringUndefiNED"; - for (humanDecayMapType::const_iterator it = humanDecayMap.begin(); it!=humanDecayMap.end(); ++it) { - humanDecayHistograms -> SetBinLabel1D(it->first,it->second); - humanDecayPileupHistograms -> SetBinLabel1D(it->first,it->second); - } - } - else if (strncmp(tmpString0,"condition",strlen("condition"))==0) { - // Selection of the predefined conditions to be applied - int iConditionTMP; - char conditionNameTMP[200]; - sscanf(&line[0],"%*s %d %s",&iConditionTMP,conditionNameTMP); - if (iConditionTMP<0) { - std::cout<<" !!! ERROR: Condition number must be >= 0 (not "<=nrConditions) { - std::cout<<" !!! ERROR: Condition number must be < "< Add it in the musrAnalysis.cxx S T O P !!!"<= -1 (not "<=nrConditions) { - std::cout<<" !!! ERROR: draw: condition number must be < "<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { - if ((*it)->IsThisHistoNamed(histoNameTMP)) { - (*it)->SetDrawListOfHistograms(iConditionTMP); - } - } - for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { - if ((*it)->IsThisHistoNamed(histoNameTMP)) { - (*it)->SetDrawListOfHistograms(iConditionTMP); - } - } - if (humanDecayHistograms->IsThisHistoNamed(histoNameTMP)) { - humanDecayHistograms->SetDrawListOfHistograms(iConditionTMP); - } - if (humanDecayPileupHistograms->IsThisHistoNamed(histoNameTMP)) { - humanDecayPileupHistograms->SetDrawListOfHistograms(iConditionTMP); - } - } - else if (strcmp(tmpString0,"counterPhaseShifts")==0) { - int nscan; int N1; float PHASE_SHIFT; char NAME1[100]; char NAME2[100]; - char *pch = line + strlen("counterPhaseShifts"); - do { - nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT); - phaseShiftMap.insert(std::pair(N1,PHASE_SHIFT/180*pi)); - std::cout<<"N1="< S T O P F O R C E D"<::const_iterator it=SampleDetIDList.begin(); it!=SampleDetIDList.end(); ++it) { - std::cout<<" "<< *it; - } - std::cout< S T O P "< S T O P "<second; - counter->SetAntiCoincidenceTimeWindowMin(lmin); - counter->SetAntiCoincidenceTimeWindowMax(lmax); - } - else if (strcmp(tmpString0,"artificiallyChangeMuDecayTime")==0) { - float min, max, mmmin, mmmax; - sscanf(&line[0],"%*s %g %g %g %g %g %g %g",&min,&max,&mmmin,&mmmax); - bool_muDecayTimeTransformation = true; - muDecayTime_t_min = min; - muDecayTime_t_max = max; - muDecayTime_Transformation_min = mmmin; - muDecayTime_Transformation_max = mmmax; - if ((muDecayTime_t_max <= muDecayTime_t_min) || (muDecayTime_Transformation_max <= muDecayTime_Transformation_min)) { - std::cout<<" ERROR! musrAnalysis: error when setting the \"artificiallyChangeMuDecayTime\" parameters! Min > Max !!!"< S T O P "< SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - funct -> SetParameter(3,p3); - funct -> FixParameter(4,xMin); - } - else if (strcmp(functionName,"funct2")==0) { - funct = new TF1("funct2","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - funct -> SetParameter(3,p3); - funct -> FixParameter(4,xMin); - funct -> SetParameter(5,p4); - } - else if (strcmp(functionName,"funct3")==0) { - funct = new TF1("funct3","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - funct -> SetParameter(3,p3); - funct -> FixParameter(4,xMin); - funct -> FixParameter(5,p4); - } - else if (strcmp(functionName,"funct4")==0) { - funct = new TF1("funct4","[3]*exp(-x/2.19703)*(1+[2]*cos(x*[0]+[1]))+[4]"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - funct -> SetParameter(3,p3); - funct -> SetParameter(4,p4); - } - else if (strcmp(functionName,"pol0")==0) { - funct = new TF1("pol0","pol0"); - } - else if (strcmp(functionName,"simpleExpoPLUSconst")==0) { - funct = new TF1("simpleExpoPLUSconst","[0]*exp(-x/2.19703)+[1]"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - } - else if (strcmp(functionName,"TFieldCosPLUSbg")==0) { - funct = new TF1("TFieldCosPLUSbg","[3]*(1+[2]*cos(x*[0]+[1]))+[4]*exp(x/2.19703)"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - funct -> SetParameter(3,p3); - funct -> SetParameter(4,p4); - } - else if (strcmp(functionName,"TFieldCos")==0) { - funct = new TF1("TFieldCos","[3]*(1+[2]*cos(x*[0]+[1]))"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - funct -> SetParameter(3,p3); - } - else if (strcmp(functionName,"rotFrameTime20")==0) { - // funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) "); - funct = new TF1("rotFrameTime20","[2]*cos(x*[0]+[1]) "); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - } - else if (strcmp(functionName,"gaus")==0) { - funct = new TF1("gaus","gaus"); - funct -> SetParameter(0,p0); - funct -> SetParameter(1,p1); - funct -> SetParameter(2,p2); - } - else { - std::cout<<"musrAnalysis::ReadInInputParameters: function \""< S T O P"<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { - if ((*it)->IsThisHistoNamed(histoName)) { - tmpHistograms = *it; - } - } - if (tmpHistograms==NULL) { - std::cout<<"\n\nmusrAnalysis::ReadInInputParameters: Fit required, but the corresponding histogram " - <<"not found!"< S T O P "<AssignFunction(funct, functOption, xMin, xMax); - } - - else { - //------------------------------------------------------------------ - // Second part of the configuration file - // - // Replace spaces between quotation marks by '_' - // and remove quatation marks and semicolon by ' '. - // Remember position of the beginning and end of the defintion of the coincidence counters - int nrOfSemicolons=0, beginningOfCoincidenceArray=0, endOfCoincidenceArray =0; - for (Int_t i=0; i<(int)strlen(line); i++) { - if (line[i]=='"') { - line[i] = ' '; - Int_t ii=i; - do { - ii++; - if (isspace(line[ii])) line[ii] = '_'; - } while (line[ii]!='"'); - line[ii]= ' '; - } - else if (line[i] == ';') { - line[i] = ' '; - nrOfSemicolons++; - if (nrOfSemicolons==5) beginningOfCoincidenceArray=i; - if (nrOfSemicolons==6) endOfCoincidenceArray=i; - } - } - - // Read in the values from the configuration file - int tmpChannel=-1, tmpTimeShift=0; - char tmpName[200], tmpType[200]; - float tmpThreshold =0.; - sscanf(&line[0],"%d %s %s %g %d",&tmpChannel,tmpName,tmpType,&tmpThreshold,&tmpTimeShift); - musrCounter* counter = new musrCounter(tmpChannel,tmpName,tmpType[0],tmpThreshold,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]=='K') {kCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} - else if (tmpType[0]=='V') {vCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} - - // Read in the values for the coincidence detectors - // (Tricky while the number of detectors in coincidence is not known) - if ((endOfCoincidenceArray>0)&&(endOfCoincidenceArray!=beginningOfCoincidenceArray)) { - char tmpCoinc[100]; strcpy(tmpCoinc,"987654321"); - int jj=beginningOfCoincidenceArray; - do { - do {jj++;} while (isspace(line[jj])); // skip space characters; - if (jj>endOfCoincidenceArray) {break;} - sscanf(&line[jj],"%s",tmpCoinc); - jj+=strlen(tmpCoinc); - int iCoinc=atoi(tmpCoinc); - if (iCoinc==987654321) continue; - tmpCoincidenceMultimap.insert(std::pair(tmpChannel,iCoinc)); - // tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc; - } while (jjSetTDChistogram(tmpHistoName,tmpT0,tmpT1,tmpT2,tmpHistoaddNr,tmpHistoNameAdd); - } - } - } - // std::cout << line; - } - } - std::cout << "INSTRUMENT=" << instrument << std::endl; - std::cout << "DESCRIPTION=" << description << std::endl; - std::cout << "TYPE=" << tdchwtype << std::endl; - std::cout << "RESOLUTION="<second)->GetCounterNr(); - for (std::multimap::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) { - if ((*itCoinc).first == iChanNr) { - int iChn2 = itCoinc->second; - counterMapType::const_iterator itTwo = allCounterMap.find(abs(iChn2)); - if (itTwo==allCounterMap.end()) { - std::cout<<" Pointer to coincidence counter ("< S T O P"<second)->SetCoincidenceCounter(itTwo->second,iChn2); - } - } - } - } - - for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) { - char DetectorType = (it->second)->GetCounterType(); - if (DetectorType=='M') { - (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin); - // (it->second)->SetCoincidenceTimeWindow_M(pileupWindowBinMin,pileupWindowBinMax); - (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin); - (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin); - } - else if (DetectorType=='P') { - (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin); - // (it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax); - (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin); - (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin); - } - // else if (DetectorType=='V') { - // (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin); - // (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin); - // } - } - - // for (int j=0; j::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) { - // std::cout << " houby [" << (*itCoinc).first << ", " << (*itCoinc).second << "]" << std::endl; - - // - // if (tmpCoincidenceMultimap[iChanNr][j]!=0) { -// int j = (*itCoinc).first; -// counterMapType::const_iterator itTwo = allCounterMap.find(j); -// if (itTwo==allCounterMap.end()) { -// std::cout<<" Pointer to coincidence counter ("< S T O P"<second; -// (it->second)->SetCoincidenceCounter(itTwo->second,iChn); -// } -// // } -// } -// } - - - - // - // - // char tmpString1[100]="Unset",tmpString2[100]="Unset",tmpString3[100]="Unset"; - // sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2); - -} - -//================================================================ - -void musrAnalysis::CreateHistograms() { - std::cout<<"musrAnalysis::CreateHistograms()"<pileupWindowMax)? 3*dataWindowMax: 3*pileupWindowMax; - std::cout<<"musrAnalysis::CreateHistograms(): safeTimeWindow = "<GetEntry(1); Double_t fieldValue = fieldNomVal[0]; if (nFieldNomVal>1) fieldValue += fieldNomVal[1]; - // omega = 851.372*fabs(fieldValue); - // omega = 851.372*fieldValue; // value set originally in this program - // omega = 850.62*fieldValue; // value used in Geant ? - omega = 851.610*fieldValue; // value from the fits of the data - - hInfo->Fill(1, fieldValue); - hInfo->Fill(6, runID); - hInfo->Fill(7, hGeantParameters->GetBinContent(7)); - std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<DrawTH1D("hist",0); - // (*it)->DrawTH1D("hist",1); - } - for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { - (*it)->DrawTH2DdrawList("hist"); - } - for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) { - (it->second)->DrawTDChistogram(); - } - if (humanDecayHistograms) { - // humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,indexHuman,nHumanDecayArray); - humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,humanDecayMap,humanDecayMultimap); - humanDecayHistograms->DrawTH1DdrawList("text"); - humanDecayPileupHistograms->FillHumanDecayArray(motherOfHumanDecayPileupHistograms,humanDecayMap,humanDecayMultimap); - humanDecayPileupHistograms->DrawTH1DdrawList("text"); - } - - Long64_t numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates(); - Long64_t numberOfMuonCandidatesAfterVK = mCounter->GetNumberOfMuonCandidatesAfterVK(); - Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval = mCounter->GetNumberOfMuonCandidatesAfterVKandDoubleHitRemoval(); - Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000; - hInfo->Fill(5,(Double_t) numberOfMuonCandidates); // number of "triggers" - hInfo->Fill(4,(Double_t) numberOfMuonCandidatesAfterVK); // number of "triggers" - hInfo->Fill(3,durationOfExperiment); - hInfo->Fill(8,(Double_t) numberOfMuonCandidatesAfterVKandDoubleHitRemoval); - //============================== - // Write out the histograms - std::cout<<"musrAnalysis::SaveHistograms()"<GetList()); - while ( TObject *obj = next() ) { - if (obj->InheritsFrom(TH1::Class()) ) {obj->Write();} - } - fHistograms->Close(); - delete fHistograms; - - // if (musrCounter::bool_WriteDataToDumpFile==true) { - // musrCounter::dumpFile.close(); - // } - - if (musrCounter::bool_WriteDataToDumpFile) { - if (myWriteDump!=NULL) { - delete myWriteDump; - } - } - - //============================== - - std::cout<<"==================================================================="<0) {std::cout<<"DEBUGEVENT "<lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)GetEntry(iiiEntry); InitialiseEvent(); - if (bool_debugingRequired) { - if (debugEventMap[eventID]>2) PrintHitsInAllCounters(); - // if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<myPrintThisCounter(eventID); - // } - if (bool_debugingRequired) { - if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<1) {std::cout<<"DEBUGEVENT "<1) {std::cout<<"____________________________________________________________________"<::max()/2); - if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere; - // RewindAllTimeInfo(rewindTimeBins); - RewindAllTimeInfo(); - } -} - -//================================================================ - -//void musrAnalysis::RemoveOldHitsFromCounters(Double_t timeLimit) { -void musrAnalysis::RemoveOldHitsFromCounters(Long64_t timeBinLimit) { - // Loop over all Counters and remove hits that happen at or before the "timeBinLimit". - // std::cout<<"musrAnalysis::RemoveOldHitsFromCounters: timeBinLimit="<RemoveHitsInCounter(timeBinLimit); - } -} - -//================================================================ - -//void musrAnalysis::RewindAllTimeInfo(Double_t timeToRewind) { -//void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) { -void musrAnalysis::RewindAllTimeInfo() { - Long64_t timeBinsToRewind = rewindTimeBins; - if (musrCounter::bool_WriteDataToDumpFile) { - // Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution); - // mCounter->CheckClockInfo(currentTimeBin); - mCounter->WriteRewindIntoDumpFile(); - // musrCounter::previousClock -= timeBinsToRewind; - } - currentTime -= timeBinsToRewind*tdcresolution; - nextUnfilledEventTime -= timeBinsToRewind*tdcresolution; - numberOfRewinds++; - // Loop over all timing information and rewind it by "timeToRewind"; - // std::cout<<"musrAnalysis::RewindAllTimeInfo()"<RewindHitsInCounter(timeBinsToRewind); - (*it).second->RewindHitsInCounter(); - } -} - -//================================================================ - -void musrAnalysis::PrintHitsInAllCounters() { - // std::cout<<"___________________\n"; - for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) { - (*it).second->myPrintThisCounter(eventID,0); - } -} -//================================================================ - -void musrAnalysis::FillHistograms(Int_t iiiEntry) { - // std::cout<<"musrAnalysis::FillHistograms() event="< dataBinMax !!! - pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,positronBinMax,timeBin1,timeBin2,timeBinDoubleHit,posEntry,idetP,idetP_ID,idetP_edep,pos_detID_doubleHit_INT); - pos_detID = Double_t(idetP_ID); - pos_detID_doubleHit = Double_t(pos_detID_doubleHit_INT); - if (debugEventMap[eventID]>2) { - if (pCounterHitExistsForThisEventID) {std::cout<<"FillHistograms: GOOD positron candidate found: timeBin1="< -1000; - muonTriggered_gen_AND_muonDecayedInSample_gen = muonTriggered_gen && muonDecayedInSample_gen; - muonTriggered_det = mCounterHitExistsForThisEventID; - positronHit_det = pCounterHitExistsForThisEventID; - goodEvent_det = muonTriggered_det && positronHit_det; - goodEvent_gen = (muDecayTime>-999)&&(muM0Time>-999); - 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() - - goodEvent_F_det = goodEvent_det && ( (find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID)) != F_posCounterList.end() ); - goodEvent_B_det = goodEvent_det && ( (find(B_posCounterList.begin(), B_posCounterList.end(), idetP_ID)) != B_posCounterList.end() ); - goodEvent_U_det = goodEvent_det && ( (find(U_posCounterList.begin(), U_posCounterList.end(), idetP_ID)) != U_posCounterList.end() ); - goodEvent_D_det = goodEvent_det && ( (find(D_posCounterList.begin(), D_posCounterList.end(), idetP_ID)) != D_posCounterList.end() ); - goodEvent_L_det = goodEvent_det && ( (find(L_posCounterList.begin(), L_posCounterList.end(), idetP_ID)) != L_posCounterList.end() ); - goodEvent_R_det = goodEvent_det && ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() ); - // std::cout<<"goodEvent_F_det="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { - (*it)->FillTH1D(wght,&condition[0]); - } - for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { - (*it)->FillTH2D(wght,&condition[0]); - } - 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!!!"<GetBranch("weight") ) ? weight : 1.; -} - -//================================================================ -Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) { - // std::cout<<"musrAnalysis::PreprocessEvent()"<GetEntry(iEn); InitialiseEvent(); - - // Clone some channels into different one, if requested by user - // (This is usefull when e.g. user splits a signal from a veto - // and uses it in two different ways - e.g. once for vetoing - // muons, and second (with a different threshold) for validating - // a positron candidate. This is initiated by the - // keyword "CLONECHANNEL" in the *.v1190 file - - - if (bool_clonedChannelsMultimap_NotEmpty) { - // std::cout<<"det_n="< ret = clonedChannelsMultimap.equal_range(chNumTMP); - for (clonedChannelsMultimapType::const_iterator ittt=ret.first; ittt!=ret.second; ++ittt) { - // std::cout << " ittt->second=" << ittt->second; - int chNumNewTMP = ittt->second; - det_ID[det_n] = chNumNewTMP; - det_edep[det_n] = det_edep[i]; - det_edep_el[det_n] = det_edep_el[i]; - det_edep_pos[det_n] = det_edep_pos[i]; - det_edep_gam[det_n] = det_edep_gam[i]; - det_edep_mup[det_n] = det_edep_mup[i]; - det_nsteps[det_n] = det_nsteps[i]; - det_length[det_n] = det_length[i]; - det_time_start[det_n] = det_time_start[i]; - det_time_end[det_n] = det_time_end[i]; - det_x[det_n] = det_x[i]; - det_y[det_n] = det_y[i]; - det_z[det_n] = det_z[i]; - det_kine[det_n] = det_kine[i]; - det_VrtxKine[det_n] = det_VrtxKine[i]; - det_VrtxX[det_n] = det_VrtxX[i]; - det_VrtxY[det_n] = det_VrtxY[i]; - det_VrtxZ[det_n] = det_VrtxZ[i]; - det_VrtxVolID[det_n] = det_VrtxVolID[i]; - det_VrtxProcID[det_n] = det_VrtxProcID[i]; - det_VrtxTrackID[det_n] = det_VrtxTrackID[i]; - det_VrtxParticleID[det_n]= det_VrtxParticleID[i]; - det_VvvKine[det_n] = det_VvvKine[i]; - det_VvvX[det_n] = det_VvvX[i]; - det_VvvY[det_n] = det_VvvY[i]; - det_VvvZ[det_n] = det_VvvZ[i]; - det_VvvVolID[det_n] = det_VvvVolID[i]; - det_VvvProcID[det_n] = det_VvvProcID[i]; - det_VvvTrackID[det_n] = det_VvvTrackID[i]; - det_VvvParticleID[det_n] = det_VvvParticleID[i]; - det_n++; - - } - } - } - } - - // std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<::iterator it; - it = allCounterMap.find(det_ID[i]); - if (it==allCounterMap.end()) { - // std::cout<<"Active detector with det_ID="<FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i],eventID); - } - } - - // std::cout<<"lastPreprocessedEntry+1="<GetNextGoodMuon(evID,timeBinMin,timeBin0,kEntry,idet,idetID,idetEdep); -// if (!mCounterHitCanditateExists) return false; -// // Check for other muons within the pileup window: -// if ( mCounter->CheckForPileupMuons(timeBin0,kEntry) ) return false; // This muon candidate is killed due to a double hit rejection. -// return true; -//} -// -//================================================================ -Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t positronBinMax, Long64_t& tBin1, Long64_t& tBin2, Long64_t& tBinDoubleHit, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Int_t& idetP_ID_doubleHit) { - - if (bool_debugingRequired) { - if (debugEventMap[eventID]>4) {std::cout<<"PositronCounterHit: pCounterMap.size()="<second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); - positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMin,positronBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); - if (positronQuality==3) { // double hit was found in the same counter - if (debugEventMap[eventID]>3) {std::cout<<"PositronCounterHit: positron candidate killed - double hit in the same counter"<3) {std::cout<<"PositronCounterHit: positron candidate killed - double hit in a different counter"<first; - } - std::cout<second); - } - std::cout<-1000.00001) ) return -1000.; - if ( (x<-999.99999) && (x>-1000.00001) ) return -1000.; - if ( (x==0) && (y==0) ) return -1000.; - return atan2(y,x); -} -//================================================================ -Double_t musrAnalysis::deltaAngle(Double_t alpha, Double_t beta) { - // Calculates the difference between angle alpha and beta. - // The angles alpha and beta are in degrees. - // The difference will be in the range of (-180,180> degrees. - if ((alpha<-998.)||(beta<-998.)) {return -1001.;} // one of the angles was undefined; - Double_t delta = alpha - beta; - if (delta<=-180) {delta+=360;} - else {if (delta>180) delta-=360;} - return delta; -} -//================================================================ +#define musrAnalysis_cxx +#include "musrAnalysis.hh" +#include "musrTH.hh" +//#include +#include +#include +#include +#include +//#include + + +musrWriteDump* musrAnalysis::myWriteDump = NULL; + +void musrAnalysis::Loop(char* runChar, char* v1190FileName, Int_t nrEvents) +{ +// In a ROOT session, you can do: +// Root > .L musrAnalysis.C +// Root > musrAnalysis t +// Root > t.GetEntry(12); // Fill t data members with entry number 12 +// Root > t.Show(); // Show values of entry 12 +// Root > t.Show(16); // Read and show values of entry 16 +// Root > t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + std::cerr<<"musrSimAnalysis::Loop() : Analysing run "<GetEntriesFast(); + if (nrEvents!=0) nentries = nrEvents; + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); InitialiseEvent(); nbytes += nb; + // if (Cut(ientry) < 0) continue; + if (jentry%10000==0) std::cout << "Analysing event nr. jentry=" << jentry << " and found " << numberOfGoodMuons << " muons so far" << std::endl; + AnalyseEvent(jentry); + } + SaveHistograms(runChar,v1190FileName); +} + + + +//================================================================ + +void musrAnalysis::ReadInInputParameters(char* charV1190FileName) { + std::cout<<"musrAnalysis::ReadInInputParameters()"< tmpCoincidenceMultimap; + // for (int i=0; i("a", 1)); + // + // } + // } + while (!feof(fSteeringFile)) { + fgets(line,10000,fSteeringFile); + if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')&&(line[0]!='$')&&(line[0]!='!')) { + char tmpString0[200]="Unset"; + char tmpString1[200]="Unset"; + sscanf(&line[0],"%s",tmpString0); + if (strncmp(tmpString0,"INSTRUMENT=",strlen("INSTRUMENT="))==0) { + instrument=&line[strlen("INSTRUMENT=")]; + } + else if (strncmp(tmpString0,"DESCRIPTION=",strlen("DESCRIPTION="))==0) { + description=&line[strlen("DESCRIPTION=")]; + } + else if (strncmp(tmpString0,"TYPE=",strlen("TYPE="))==0) { + tdchwtype=&line[strlen("TYPE=")]; + } + else if (strncmp(tmpString0,"RESOLUTION=",strlen("RESOLUTION="))==0) { + sscanf(&line[strlen("RESOLUTION=")],"%g",&tmpfvar); + tdcresolution = tmpfvar * picosecond; + if (tdcresolution<1*picosecond) { + std::cout<<"\n !!! ERROR !!! Time resolution of TDC must be larger than 1 ps (you requested "< too low time resolution might lead to numerical problems in musrCounter ==> S T O P"<(ichannel_orig_tmp,ichannel_new_tmp)); + } + else if (strcmp(tmpString0,"DEBUGEVENT")==0) { + int ieventToDebug_tmp, iLevelToDebug_tmp; + sscanf(&line[0],"%*s %d %d",&ieventToDebug_tmp,&iLevelToDebug_tmp); + bool_debugingRequired=true; + if (ieventToDebug_tmp>-1) { + debugEventMap.insert(std::pair(ieventToDebug_tmp,iLevelToDebug_tmp)); + } + else { + for (int j=0; j<=-ieventToDebug_tmp;j++) { + debugEventMap.insert(std::pair(j,iLevelToDebug_tmp)); + } + } + } + else if (strcmp(tmpString0,"WRITE_OUT_DUMP_FILE")==0) { + int clock_channelID; + int max_timeBin_jitter; + // int clock_interval_tmp; + // sscanf(&line[0],"%*s %s %d %d",tmpString1,&clock_channelID_tmp,&clock_interval_tmp); + sscanf(&line[0],"%*s %s %d %d",tmpString1,&clock_channelID,&max_timeBin_jitter); + // clock_channelID = clock_channelID_tmp; + // clock_interval = clock_interval_tmp; + musrCounter::bool_WriteDataToDumpFile = true; + // std::cout<<"tmpString1="< STOP !!!"< S T O P !!!"< ignored, histogram not defined! S T O P !!!"< ignored, histogram not defined! S T O P !!!"<second; + if (strcmp(tmpString0,"musrTH1D")==0) { + // std::cout<<"pointerToVariable1="<SetRotatingReferenceFrame(rot_ref_frequency, rot_ref_phase); + hInfo->Fill(21, rot_ref_frequency); // value may be overwritten - just the last rot. ref. frame will be saved + hInfo->Fill(22, rot_ref_phase); // value may be overwritten - just the last rot. ref. frame will be saved + } + else if (strcmp(furtherOption,"correctexpdecay")==0) { + myTH->SetExpDecayCorrection(); + } + } + else { + Double_t* pointerToVariable2 = iter2->second; + musrTH* myTH = new musrTH("2D",histoName,histoTitle,nrBinsX,minX,maxX,nrBinsY,minY,maxY,pointerToVariable1,pointerToVariable2); + listOfAllHistograms2D.push_back(myTH); + } + } + else if (strncmp(tmpString0,"humanDecayHistograms",strlen("humanDecayHistograms"))==0) { + char motherHistoName[200]; + char motherHistoPileupName[200]; + sscanf(&line[0],"%*s %s %s",motherHistoName,motherHistoPileupName); + std::cout<<" Defining humanDecayHistograms (motherHistoName="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + if ((*it)->IsThisHistoNamed(motherHistoName)) { + motherOfHumanDecayHistograms = *it; + } + if ((*it)->IsThisHistoNamed(motherHistoPileupName)) { + motherOfHumanDecayPileupHistograms = *it; + } + } + if (motherOfHumanDecayHistograms==NULL) { + std::cout<<"\n\nmusrAnalysis::ReadInInputParameters Mother histogram (motherOfHumanDecayHistograms)" + <<"of the \"Human decay histogram\" not found!"< S T O P "< S T O P "<second)==(std::string) NAME1) { + nameDefined = true; + iBinHuman = it->first; + } + } + if (!nameDefined) { + nBinHuman++; + iBinHuman=nBinHuman; + humanDecayMap.insert(std::pair(iBinHuman,NAME1)); + } + humanDecayMultimap.insert(std::pair(iBinHuman,N1)); + char* pch2 = strstr(pch,NAME1)+strlen(NAME1); + pch=pch2; + nscan = sscanf(pch,"%d %s",&N1,NAME1); + // indexHuman++; + } while (nscan==2); + + Int_t nrBinsHumanHist=humanDecayMap.size(); + // Int_t nrBinsHumanHist=0; + // std::string oldString="blaBlaoldStringUndefiNED"; + // for (humanDecayMultimapType::const_iterator it = humanDecayMultimap.begin(); it!=humanDecayMultimap.end(); ++it) { + // if (oldString!=(it->first)) nrBinsHumanHist++; + // } + std::cout<<"\nnrBinsHumanHist="<SetBinLabel1D(k+1,sHumanDecayArray[k]); + // } + // oldString="blaBlaoldStringUndefiNED"; + for (humanDecayMapType::const_iterator it = humanDecayMap.begin(); it!=humanDecayMap.end(); ++it) { + humanDecayHistograms -> SetBinLabel1D(it->first,it->second); + humanDecayPileupHistograms -> SetBinLabel1D(it->first,it->second); + } + } + else if (strncmp(tmpString0,"condition",strlen("condition"))==0) { + // Selection of the predefined conditions to be applied + int iConditionTMP; + char conditionNameTMP[200]; + sscanf(&line[0],"%*s %d %s",&iConditionTMP,conditionNameTMP); + if (iConditionTMP<0) { + std::cout<<" !!! ERROR: Condition number must be >= 0 (not "<=nrConditions) { + std::cout<<" !!! ERROR: Condition number must be < "< Add it in the musrAnalysis.cxx S T O P !!!"<= -1 (not "<=nrConditions) { + std::cout<<" !!! ERROR: draw: condition number must be < "<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + if ((*it)->IsThisHistoNamed(histoNameTMP)) { + (*it)->SetDrawListOfHistograms(iConditionTMP); + } + } + for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { + if ((*it)->IsThisHistoNamed(histoNameTMP)) { + (*it)->SetDrawListOfHistograms(iConditionTMP); + } + } + if (humanDecayHistograms->IsThisHistoNamed(histoNameTMP)) { + humanDecayHistograms->SetDrawListOfHistograms(iConditionTMP); + } + if (humanDecayPileupHistograms->IsThisHistoNamed(histoNameTMP)) { + humanDecayPileupHistograms->SetDrawListOfHistograms(iConditionTMP); + } + } + else if (strcmp(tmpString0,"counterPhaseShifts")==0) { + int nscan; int N1; float PHASE_SHIFT; char NAME1[100]; char NAME2[100]; + char *pch = line + strlen("counterPhaseShifts"); + do { + nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT); + phaseShiftMap.insert(std::pair(N1,PHASE_SHIFT/180*pi)); + std::cout<<"N1="< S T O P F O R C E D"<::const_iterator it=SampleDetIDList.begin(); it!=SampleDetIDList.end(); ++it) { + std::cout<<" "<< *it; + } + std::cout< S T O P "< S T O P "<second; + counter->SetAntiCoincidenceTimeWindowMin(lmin); + counter->SetAntiCoincidenceTimeWindowMax(lmax); + } + else if (strcmp(tmpString0,"artificiallyChangeMuDecayTime")==0) { + float min, max, mmmin, mmmax; + sscanf(&line[0],"%*s %g %g %g %g %g %g %g",&min,&max,&mmmin,&mmmax); + bool_muDecayTimeTransformation = true; + muDecayTime_t_min = min; + muDecayTime_t_max = max; + muDecayTime_Transformation_min = mmmin; + muDecayTime_Transformation_max = mmmax; + if ((muDecayTime_t_max <= muDecayTime_t_min) || (muDecayTime_Transformation_max <= muDecayTime_Transformation_min)) { + std::cout<<" ERROR! musrAnalysis: error when setting the \"artificiallyChangeMuDecayTime\" parameters! Min > Max !!!"< S T O P "< SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> FixParameter(4,xMin); + } + else if (strcmp(functionName,"funct2")==0) { + funct = new TF1("funct2","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> FixParameter(4,xMin); + funct -> SetParameter(5,p4); + } + else if (strcmp(functionName,"funct3")==0) { + funct = new TF1("funct3","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> FixParameter(4,xMin); + funct -> FixParameter(5,p4); + } + else if (strcmp(functionName,"funct4")==0) { + funct = new TF1("funct4","[3]*exp(-x/2.19703)*(1+[2]*cos(x*[0]+[1]))+[4]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> SetParameter(4,p4); + } + else if (strcmp(functionName,"pol0")==0) { + funct = new TF1("pol0","pol0"); + } + else if (strcmp(functionName,"simpleExpoPLUSconst")==0) { + funct = new TF1("simpleExpoPLUSconst","[0]*exp(-x/2.19703)+[1]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + } + else if (strcmp(functionName,"TFieldCosPLUSbg")==0) { + funct = new TF1("TFieldCosPLUSbg","[3]*(1+[2]*cos(x*[0]+[1]))+[4]*exp(x/2.19703)"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> SetParameter(4,p4); + } + else if (strcmp(functionName,"TFieldCos")==0) { + funct = new TF1("TFieldCos","[3]*(1+[2]*cos(x*[0]+[1]))"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + } + else if (strcmp(functionName,"rotFrameTime20")==0) { + // funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) "); + funct = new TF1("rotFrameTime20","[2]*cos(x*[0]+[1]) "); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + } + else if (strcmp(functionName,"gaus")==0) { + funct = new TF1("gaus","gaus"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + } + else { + std::cout<<"musrAnalysis::ReadInInputParameters: function \""< S T O P"<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + if ((*it)->IsThisHistoNamed(histoName)) { + tmpHistograms = *it; + } + } + if (tmpHistograms==NULL) { + std::cout<<"\n\nmusrAnalysis::ReadInInputParameters: Fit required, but the corresponding histogram " + <<"not found!"< S T O P "<AssignFunction(funct, functOption, xMin, xMax); + } + + else { + //------------------------------------------------------------------ + // Second part of the configuration file + // + // Replace spaces between quotation marks by '_' + // and remove quatation marks and semicolon by ' '. + // Remember position of the beginning and end of the defintion of the coincidence counters + int nrOfSemicolons=0, beginningOfCoincidenceArray=0, endOfCoincidenceArray =0; + for (Int_t i=0; i<(int)strlen(line); i++) { + if (line[i]=='"') { + line[i] = ' '; + Int_t ii=i; + do { + ii++; + if (isspace(line[ii])) line[ii] = '_'; + } while (line[ii]!='"'); + line[ii]= ' '; + } + else if (line[i] == ';') { + line[i] = ' '; + nrOfSemicolons++; + if (nrOfSemicolons==5) beginningOfCoincidenceArray=i; + if (nrOfSemicolons==6) endOfCoincidenceArray=i; + } + } + + // Read in the values from the configuration file + int tmpChannel=-1, tmpTimeShift=0; + char tmpName[200], tmpType[200]; + 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); + 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;} + 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]=='V') {vCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} + + // Read in the values for the coincidence detectors + // (Tricky while the number of detectors in coincidence is not known) + if ((endOfCoincidenceArray>0)&&(endOfCoincidenceArray!=beginningOfCoincidenceArray)) { + char tmpCoinc[100]; strcpy(tmpCoinc,"987654321"); + int jj=beginningOfCoincidenceArray; + do { + do {jj++;} while (isspace(line[jj])); // skip space characters; + if (jj>endOfCoincidenceArray) {break;} + sscanf(&line[jj],"%s",tmpCoinc); + jj+=strlen(tmpCoinc); + int iCoinc=atoi(tmpCoinc); + if (iCoinc==987654321) continue; + tmpCoincidenceMultimap.insert(std::pair(tmpChannel,iCoinc)); + // tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc; + } while (jjSetTDChistogram(tmpHistoName,tmpT0,tmpT1,tmpT2,tmpHistoaddNr,tmpHistoNameAdd); + } + } + } + // std::cout << line; + } + } + std::cout << "INSTRUMENT=" << instrument << std::endl; + std::cout << "DESCRIPTION=" << description << std::endl; + std::cout << "TYPE=" << tdchwtype << std::endl; + std::cout << "RESOLUTION="<second)->GetCounterNr(); + for (std::multimap::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) { + if ((*itCoinc).first == iChanNr) { + int iChn2 = itCoinc->second; + counterMapType::const_iterator itTwo = allCounterMap.find(abs(iChn2)); + if (itTwo==allCounterMap.end()) { + std::cout<<" Pointer to coincidence counter ("< S T O P"<second)->SetCoincidenceCounter(itTwo->second,iChn2); + } + } + } + } + + for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) { + char DetectorType = (it->second)->GetCounterType(); + if (DetectorType=='M') { + (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin); + // (it->second)->SetCoincidenceTimeWindow_M(pileupWindowBinMin,pileupWindowBinMax); + (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin); + (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin); + } + else if (DetectorType=='P') { + (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin); + // (it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax); + (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin); + (it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin); + } + // else if (DetectorType=='V') { + // (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin); + // (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin); + // } + } + + // for (int j=0; j::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) { + // std::cout << " houby [" << (*itCoinc).first << ", " << (*itCoinc).second << "]" << std::endl; + + // + // if (tmpCoincidenceMultimap[iChanNr][j]!=0) { +// int j = (*itCoinc).first; +// counterMapType::const_iterator itTwo = allCounterMap.find(j); +// if (itTwo==allCounterMap.end()) { +// std::cout<<" Pointer to coincidence counter ("< S T O P"<second; +// (it->second)->SetCoincidenceCounter(itTwo->second,iChn); +// } +// // } +// } +// } + + + + // + // + // char tmpString1[100]="Unset",tmpString2[100]="Unset",tmpString3[100]="Unset"; + // sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2); + +} + +//================================================================ + +void musrAnalysis::CreateHistograms() { + std::cout<<"musrAnalysis::CreateHistograms()"<pileupWindowMax)? 3*dataWindowMax: 3*pileupWindowMax; + std::cout<<"musrAnalysis::CreateHistograms(): safeTimeWindow = "<GetEntry(1); InitialiseEvent(); Double_t fieldValue = fieldNomVal[0]; if (nFieldNomVal>1) fieldValue += fieldNomVal[1]; + // omega = 851.372*fabs(fieldValue); + // omega = 851.372*fieldValue; // value set originally in this program + // omega = 850.62*fieldValue; // value used in Geant ? + omega = 851.610*fieldValue; // value from the fits of the data + + hInfo->Fill(1, fieldValue); + hInfo->Fill(6, runID); + hInfo->Fill(7, hGeantParameters->GetBinContent(7)); + std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<DrawTH1D("hist",0); + // (*it)->DrawTH1D("hist",1); + } + for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { + (*it)->DrawTH2DdrawList("hist"); + } + if(musrMode == 'D') { + for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) { + (it->second)->DrawTDChistogram(); + } + } + if (humanDecayHistograms) { + // humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,indexHuman,nHumanDecayArray); + humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,humanDecayMap,humanDecayMultimap); + humanDecayHistograms->DrawTH1DdrawList("text"); + humanDecayPileupHistograms->FillHumanDecayArray(motherOfHumanDecayPileupHistograms,humanDecayMap,humanDecayMultimap); + humanDecayPileupHistograms->DrawTH1DdrawList("text"); + } + + // need to re-define these for Pulsed acquisition (and anything else) + Long64_t numberOfMuonCandidates; + Long64_t numberOfMuonCandidatesAfterVK; + 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(4,(Double_t) numberOfMuonCandidatesAfterVK); // number of "triggers" + hInfo->Fill(3,durationOfExperiment); + hInfo->Fill(8,(Double_t) numberOfMuonCandidatesAfterVKandDoubleHitRemoval); + } + //============================== + // Write out the histograms + std::cout<<"musrAnalysis::SaveHistograms()"<GetList()); + while ( TObject *obj = next() ) { + if (obj->InheritsFrom(TH1::Class()) ) {obj->Write();} + } + fHistograms->Close(); + delete fHistograms; + + // if (musrCounter::bool_WriteDataToDumpFile==true) { + // musrCounter::dumpFile.close(); + // } + + if (musrCounter::bool_WriteDataToDumpFile) { + if (myWriteDump!=NULL) { + delete myWriteDump; + } + } + + // condition counts + std::cout << "Number of times each condition was recorded" << std::endl; + for(int i=0; i0) {std::cout<<"DEBUGEVENT "<lastPreprocessedEntry)||((nextUnfilledEventTime=nentries) && (nextUnfilledEventTimelastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)GetEntry(iiiEntry); InitialiseEvent(); + if (bool_debugingRequired) { + if (debugEventMap[eventID]>2) PrintHitsInAllCounters(); + // if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<myPrintThisCounter(eventID); + // } + if (bool_debugingRequired) { + if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<1) {std::cout<<"DEBUGEVENT "<1) {std::cout<<"____________________________________________________________________"<0)) { + nextUnfilledEventTime=0.0; + // std::cout << "resetting NextUnfilledEventTime"<::max()/2); + if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere; + // RewindAllTimeInfo(rewindTimeBins); + RewindAllTimeInfo(); + } + } +} + +//================================================================ + +//void musrAnalysis::RemoveOldHitsFromCounters(Double_t timeLimit) { +void musrAnalysis::RemoveOldHitsFromCounters(Long64_t timeBinLimit) { + // Loop over all Counters and remove hits that happen at or before the "timeBinLimit". + // std::cout<<"musrAnalysis::RemoveOldHitsFromCounters: timeBinLimit="<RemoveHitsInCounter(timeBinLimit); + } +} + +//================================================================ + +//void musrAnalysis::RewindAllTimeInfo(Double_t timeToRewind) { +//void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) { +void musrAnalysis::RewindAllTimeInfo() { + Long64_t timeBinsToRewind = rewindTimeBins; + if (musrCounter::bool_WriteDataToDumpFile) { + // Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution); + // mCounter->CheckClockInfo(currentTimeBin); + mCounter->WriteRewindIntoDumpFile(); + // musrCounter::previousClock -= timeBinsToRewind; + } + currentTime -= timeBinsToRewind*tdcresolution; + nextUnfilledEventTime -= timeBinsToRewind*tdcresolution; + numberOfRewinds++; + // Loop over all timing information and rewind it by "timeToRewind"; + // std::cout<<"musrAnalysis::RewindAllTimeInfo()"<RewindHitsInCounter(timeBinsToRewind); + (*it).second->RewindHitsInCounter(); + } +} + +//================================================================ + +void musrAnalysis::PrintHitsInAllCounters() { + // std::cout<<"___________________\n"; + for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) { + (*it).second->myPrintThisCounter(eventID,0); + } +} +//================================================================ + +void musrAnalysis::FillHistograms(Int_t iiiEntry) { + // std::cout<<"musrAnalysis::FillHistograms() event="< dataBinMax !!! + pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,positronBinMax,timeBin1,timeBin2,timeBinDoubleHit,posEntry,idetP,idetP_ID,idetP_edep,pos_detID_doubleHit_INT); + pos_detID = Double_t(idetP_ID); + pos_detID_doubleHit = Double_t(pos_detID_doubleHit_INT); + if (debugEventMap[eventID]>2) { + if (pCounterHitExistsForThisEventID) {std::cout<<"FillHistograms: GOOD positron candidate found: timeBin1="< -1000; + muonTriggered_gen_AND_muonDecayedInSample_gen = muonTriggered_gen && muonDecayedInSample_gen; + muonTriggered_det = mCounterHitExistsForThisEventID; + positronHit_det = pCounterHitExistsForThisEventID; + goodEvent_det = muonTriggered_det && positronHit_det; + goodEvent_gen = (muDecayTime>-999)&&(muM0Time>-999); + 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() + + goodEvent_F_det = goodEvent_det && ( (find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID)) != F_posCounterList.end() ); + goodEvent_B_det = goodEvent_det && ( (find(B_posCounterList.begin(), B_posCounterList.end(), idetP_ID)) != B_posCounterList.end() ); + goodEvent_U_det = goodEvent_det && ( (find(U_posCounterList.begin(), U_posCounterList.end(), idetP_ID)) != U_posCounterList.end() ); + goodEvent_D_det = goodEvent_det && ( (find(D_posCounterList.begin(), D_posCounterList.end(), idetP_ID)) != D_posCounterList.end() ); + goodEvent_L_det = goodEvent_det && ( (find(L_posCounterList.begin(), L_posCounterList.end(), idetP_ID)) != L_posCounterList.end() ); + goodEvent_R_det = goodEvent_det && ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() ); + // std::cout<<"goodEvent_F_det="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + (*it)->FillTH1D(wght,&condition[0]); + } + for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { + (*it)->FillTH2D(wght,&condition[0]); + } + 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!!!"< unrecordedEventMap; +// unrecordedEventMap unrecordedEvents; +// unrecordedEventMap::const_iterator it2; + + typedef std::multimap 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(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(i,true)); + //} + + // std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<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 ="< 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="< -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="<promptPeakWindowMin) && (det_time10 -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="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + (*it)->FillTH1D(wght,&condition[0]); + } + for(std::list::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!!!"<GetBranch("weight") ) ? weight : 1.; +} + +//================================================================ +Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) { + // std::cout<<"musrAnalysis::PreprocessEvent()"<GetEntry(iEn); InitialiseEvent(); + + // Clone some channels into different one, if requested by user + // (This is usefull when e.g. user splits a signal from a veto + // and uses it in two different ways - e.g. once for vetoing + // muons, and second (with a different threshold) for validating + // a positron candidate. This is initiated by the + // keyword "CLONECHANNEL" in the *.v1190 file + + + if (bool_clonedChannelsMultimap_NotEmpty) { + // std::cout<<"det_n="< ret = clonedChannelsMultimap.equal_range(chNumTMP); + for (clonedChannelsMultimapType::const_iterator ittt=ret.first; ittt!=ret.second; ++ittt) { + // std::cout << " ittt->second=" << ittt->second; + int chNumNewTMP = ittt->second; + det_ID[det_n] = chNumNewTMP; + det_edep[det_n] = det_edep[i]; + det_edep_el[det_n] = det_edep_el[i]; + det_edep_pos[det_n] = det_edep_pos[i]; + det_edep_gam[det_n] = det_edep_gam[i]; + det_edep_mup[det_n] = det_edep_mup[i]; + det_nsteps[det_n] = det_nsteps[i]; + det_length[det_n] = det_length[i]; + det_time_start[det_n] = det_time_start[i]; + det_time_end[det_n] = det_time_end[i]; + det_x[det_n] = det_x[i]; + det_y[det_n] = det_y[i]; + det_z[det_n] = det_z[i]; + det_kine[det_n] = det_kine[i]; + det_VrtxKine[det_n] = det_VrtxKine[i]; + det_VrtxX[det_n] = det_VrtxX[i]; + det_VrtxY[det_n] = det_VrtxY[i]; + det_VrtxZ[det_n] = det_VrtxZ[i]; + det_VrtxVolID[det_n] = det_VrtxVolID[i]; + det_VrtxProcID[det_n] = det_VrtxProcID[i]; + det_VrtxTrackID[det_n] = det_VrtxTrackID[i]; + det_VrtxParticleID[det_n]= det_VrtxParticleID[i]; + det_VvvKine[det_n] = det_VvvKine[i]; + det_VvvX[det_n] = det_VvvX[i]; + det_VvvY[det_n] = det_VvvY[i]; + det_VvvZ[det_n] = det_VvvZ[i]; + det_VvvVolID[det_n] = det_VvvVolID[i]; + det_VvvProcID[det_n] = det_VvvProcID[i]; + det_VvvTrackID[det_n] = det_VvvTrackID[i]; + det_VvvParticleID[det_n] = det_VvvParticleID[i]; + det_n++; + + } + } + } + } + + // std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<(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]; + std::map::iterator it; + it = allCounterMap.find(det_ID[i]); + if (it==allCounterMap.end()) { + std::cout<<"Active detector with det_ID="<FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i],eventID, multiCtr); + // std::cout << "pushed event into counter "<GetNextGoodMuon(evID,timeBinMin,timeBin0,kEntry,idet,idetID,idetEdep); +// if (!mCounterHitCanditateExists) return false; +// // Check for other muons within the pileup window: +// if ( mCounter->CheckForPileupMuons(timeBin0,kEntry) ) return false; // This muon candidate is killed due to a double hit rejection. +// return true; +//} +// +//================================================================ +Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t positronBinMax, Long64_t& tBin1, Long64_t& tBin2, Long64_t& tBinDoubleHit, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Int_t& idetP_ID_doubleHit) { + + if (bool_debugingRequired) { + if (debugEventMap[eventID]>4) {std::cout<<"PositronCounterHit: pCounterMap.size()="<second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); + positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMin,positronBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep); + if (positronQuality==3) { // double hit was found in the same counter + if (debugEventMap[eventID]>3) {std::cout<<"PositronCounterHit: positron candidate killed - double hit in the same counter"<3) {std::cout<<"PositronCounterHit: positron candidate killed - double hit in a different counter"<first; + } + std::cout<second); + } + std::cout<-1000.00001) ) return -1000.; + if ( (x<-999.99999) && (x>-1000.00001) ) return -1000.; + if ( (x==0) && (y==0) ) return -1000.; + return atan2(y,x); +} +//================================================================ +Double_t musrAnalysis::deltaAngle(Double_t alpha, Double_t beta) { + // Calculates the difference between angle alpha and beta. + // The angles alpha and beta are in degrees. + // The difference will be in the range of (-180,180> degrees. + if ((alpha<-998.)||(beta<-998.)) {return -1001.;} // one of the angles was undefined; + Double_t delta = alpha - beta; + if (delta<=-180) {delta+=360;} + else {if (delta>180) delta-=360;} + return delta; +} +//================================================================ diff --git a/musrSimAna/musrAnalysis.hh b/musrSimAna/musrAnalysis.hh index a330f59..38706de 100644 --- a/musrSimAna/musrAnalysis.hh +++ b/musrSimAna/musrAnalysis.hh @@ -23,6 +23,12 @@ //#include "musrTH.hh" #include +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; //#include "musrSimGlobal.hh" @@ -192,6 +198,7 @@ public : Double_t eventID_double; Double_t muDecayDetID_double; Double_t det_n_double; + Double_t det_multifirst_double; Double_t muDecayPosR; Double_t wght; Double_t det_m0edep; @@ -219,6 +226,7 @@ public : Double_t pos_Phi_MINUS_muDecayPol_Phi360; Double_t pos_detID; Double_t pos_detID_doubleHit; + Double_t pos_doubleHit_dPhi; // Double_t det_time0; // Double_t get_time0; // Double_t det_time1; @@ -229,6 +237,7 @@ public : Double_t det_time20; Double_t det_time31; Double_t det_time1_MINUS_muDecayTime; + Double_t det_multi_interval; Double_t detP_x; Double_t detP_y; Double_t detP_z; @@ -254,6 +263,7 @@ public : Bool_t alwaysTrue; Bool_t oncePerEvent; Bool_t muonDecayedInSample_gen; + Bool_t muonDecayedInSampleOnce_gen; Bool_t muonTriggered_gen; Bool_t muonTriggered_gen_AND_muonDecayedInSample_gen; Bool_t muonTriggered_det; @@ -264,6 +274,12 @@ public : Bool_t pileupEventCandidate; Bool_t pileupEvent; 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_B_det; Bool_t goodEvent_U_det; @@ -284,6 +300,25 @@ public : Bool_t promptPeakL; Bool_t promptPeakR; 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); virtual ~musrAnalysis(); @@ -298,6 +333,7 @@ public : virtual void CreateHistograms(); virtual void AnalyseEvent(Long64_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 RemoveOldHitsFromCounters(Long64_t timeBinLimit); // virtual void RewindAllTimeInfo(Double_t timeToRewind); @@ -319,11 +355,13 @@ public : TH1D* hGeantParameters; TH1D* hInfo; + typedef std::multimap doubleHitSorterType; + // typedef std::map debugEventMapType; // debugEventMapType debugEventMap; // Bool_t bool_debugingRequired; - static const Double_t pi=3.14159265358979324; + //JSL static const Double_t pi=3.14159265358979324; static musrWriteDump* myWriteDump; private: @@ -335,9 +373,10 @@ public : Int_t mdelay; Int_t pdelay; 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; - Double_t muonRateFactor; + Double_t muonRateFactor; // also mean rate for pulsed mode + Double_t muonPulseWidthFactor; // pulsed Double_t dataWindowMin; Double_t dataWindowMax; Double_t pileupWindowMin; @@ -346,16 +385,22 @@ public : Double_t promptPeakWindowMax; Int_t overallBinDelay; Bool_t boolInfinitelyLowMuonRate; - - char musrMode; // D = time diferential; I = time integral; + 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; P = pulsed time differential Double_t safeTimeWindow; Double_t currentTime; // Double_t nextEventTime; Double_t nextUnfilledEventTime; - Long64_t numberOfRewinds; + Long64_t numberOfRewinds; // JSL: number of frames for Pulsed Mode Long64_t numberOfGoodMuons; // Int_t currentEventID; Long64_t lastPreprocessedEntry; + Long64_t firstPreprocessedEntry; // JSL: first in this frame musrCounter* mCounter; typedef std::map counterMapType; counterMapType pCounterMap; @@ -367,9 +412,9 @@ public : Bool_t bool_muDecayTimeTransformation; Double_t muDecayTime_Transformation_min, muDecayTime_Transformation_max, muDecayTime_t_min, muDecayTime_t_max; - static const Double_t microsecond=1.; - static const Double_t nanosecond=0.001; - static const Double_t picosecond=0.000001; + // static const Double_t microsecond=1.; + // static const Double_t nanosecond=0.001; + // static const Double_t picosecond=0.000001; // static const Double_t rewindTime=1000000.; // static const Double_t rewindTime=1000.; // static const Long64_t rewindTimeBins=1000000000000000; // Max Long64_t can be +9,223,372,036,854,775,807 @@ -383,6 +428,7 @@ public : public: static const Int_t nrConditions = 31; Bool_t condition[nrConditions]; + Long64_t conditionCounter[nrConditions]; static Long64_t rewindTimeBins; // static Int_t clock_channelID; // static Long64_t clock_interval; @@ -488,7 +534,10 @@ musrAnalysis::musrAnalysis(TTree *tree) variableMap["posIniMomZ"]=&posIniMomZ; // variableMap["nFieldNomVal"]=&nFieldNomVal_double; // 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_multifirst"]=&det_multifirst_double; // variableMap["muDecayPosR"]=&muDecayPosR; variableMap["wght"]=&wght; @@ -517,6 +566,7 @@ musrAnalysis::musrAnalysis(TTree *tree) variableMap["pos_Phi_MINUS_muDecayPol_Phi360"]=&pos_Phi_MINUS_muDecayPol_Phi360; variableMap["pos_detID"]=&pos_detID; variableMap["pos_detID_doubleHit"]=&pos_detID_doubleHit; + variableMap["pos_doubleHit_dPhi"]=&pos_doubleHit_dPhi; // variableMap["det_time0"]=&det_time0; // variableMap["gen_time0"]=&gen_time0; // variableMap["det_time1"]=&det_time1; @@ -525,6 +575,7 @@ musrAnalysis::musrAnalysis(TTree *tree) variableMap["gen_time10"]=&gen_time10; variableMap["det_time10_MINUS_gen_time10"]=&det_time10_MINUS_gen_time10; variableMap["det_time1_MINUS_muDecayTime"]=&det_time1_MINUS_muDecayTime; + variableMap["multiHitInterval"]=&det_multi_interval; variableMap["detP_x"]=&detP_x; variableMap["detP_y"]=&detP_x; variableMap["detP_z"]=&detP_x; diff --git a/musrSimAna/musrCounter.cxx b/musrSimAna/musrCounter.cxx index 8fe3567..59efc31 100644 --- a/musrSimAna/musrCounter.cxx +++ b/musrSimAna/musrCounter.cxx @@ -1,413 +1,523 @@ -//#include -#include "musrCounter.hh" -#include "TCanvas.h" -#include "musrAnalysis.hh" - -typedef std::map debugEventMapType; -debugEventMapType debugEventMap; -Bool_t bool_debugingRequired; - -Bool_t musrCounter::bool_ignoreUnperfectMuons = true; -Bool_t musrCounter::bool_ignoreUnperfectPositrons = true; -Bool_t musrCounter::bool_WriteDataToDumpFile = false; -//Long64_t musrCounter::previousClock = -1; -//Long64_t musrCounter::CLOCK_INTERVAL = 512000; -//ofstream musrCounter::dumpFile; - -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - -musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT) { - // pointerToAnalysis=this; - counterNr = CHANNEL_NR; - strcpy(counterName,CHANNEL_NAME); - counterType = CHANNEL_TYPE; - couterEThreshold = E_THRESH; - counterTimeShift = (Long64_t) TIME_SHIFT; - std::cout<<"musrCounter::musrCounter: Creating counter "<Fill(variable,vaha); -} - -//================================================================ -void musrCounter::DrawTDChistogram() { - char canvasName[501]; - sprintf(canvasName,"c%s",TDC_histoName); - TCanvas* cTmp = new TCanvas(canvasName,canvasName); - histTDC->Draw(); -} - -//================================================================ -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){ - //cDEL std::cout<<"FillHitInCounter: timeBin="<=couterEThreshold) { - - hitInfo* hInfo = new hitInfo(kEntry,eveID,iDet,detectorID,edep,timeBin2-counterTimeShift); - //cDEL std::cout<first; - hitInfo* tempEvnr= it->second; - tempEvnr->RewindTimeBin2(musrAnalysis::rewindTimeBins); - hitMap_TMP.insert( std::pair(tempBinT-musrAnalysis::rewindTimeBins,tempEvnr) ); - } - hitMap.swap(hitMap_TMP); -} - -//================================================================ -Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, char motherCounter){ - // timeBin ... time bin, at which the coincidence is searched - - if (hitMap.empty()) return false; - Long64_t timeBinMin; - Long64_t timeBinMax; - - // If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors). - // Otherwise use timeBinMinimum and timeBinMaximum (e.g.coincidence of a positron counter with other positron counters). - if (counterType == 'V') timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector - else if (motherCounter=='M') timeBinMin = timeBin + coincidenceTimeWindowMin_M; // this is coinc. detector connected to M - else if (motherCounter=='P') timeBinMin = timeBin + coincidenceTimeWindowMin_P; // this is coinc. detector connected to P - - if (counterType == 'V') timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector - else if (motherCounter=='M') timeBinMax = timeBin + coincidenceTimeWindowMax_M; // this is coinc. detector connected to M - else if (motherCounter=='P') timeBinMax = timeBin + coincidenceTimeWindowMax_P; // this is coinc. detector connected to P - - for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { - Long64_t timeBinOfCount_tmp = it->first; - if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) { - // if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) { - return true; - // } - } - else if (timeBinOfCount_tmp > timeBinMax) return false; - } - return false; -} - -//================================================================ -Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep) { - // This function searches for a good muon, i.e. the muon - // 1) belongs to the currently analysed event - // 2) is in coincidence with all required coincidence detectors - // 3) is not in coincidence with veto detectors - // INPUT PARAMETERS: evtID, timeBinMin - // OUTPUT PARAMETERS: timeBinOfNextHit - // - // Loop over the hits in the counter - // std::cout<<" musrCounter::GetNextGoodMuon timeBinMin="< S T O P !!!\n"; exit(1);} - - std::list it_muon_hits_to_be_deleted; - for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { - - Long64_t timeBinOfCount_tmp = it->first; - timeBinOfNextHit = timeBinOfCount_tmp; - Int_t modulo = timeBinOfNextHit % 524288; - if (timeBinOfCount_tmp <= timeBinMin) continue; // This hit was already processed previously ==> skip it - - Int_t eventNumber = (it->second)->eventIDnumber; - if (eventNumber!=evtID) continue; // This trigger hit does not correspond to the currently processed event - // ==> skip it, because it was already proceesed or will be processed in future - numberOfMuonCandidates++; - if (bool_debugingRequired) {if (debugEventMap[evtID]>3) std::cout<<"GetNextGoodMuon: muon candidate found ("<second)->IsInCoincidence(timeBinOfCount_tmp,'M') )) { // no coincidence found ==> skip hit - // if (bool_ignoreUnperfectMuons) hitMap.erase(it); - if (bool_ignoreUnperfectMuons) it_muon_hits_to_be_deleted.push_back(it); - bool_coincidenceConditions = false; - if (bool_debugingRequired) {if (debugEventMap[evtID]>3) std::cout<<"GetNextGoodMuon: muon candidate not in koincidence ("<second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) { // coincidence with veto found ==> skip hit - // if (bool_ignoreUnperfectMuons) hitMap.erase(it); - if (bool_ignoreUnperfectMuons) it_muon_hits_to_be_deleted.push_back(it); - bool_coincidenceConditions = false; - if (bool_debugingRequired) {if (debugEventMap[evtID]>3) std::cout<<"GetNextGoodMuon: muon candidate vetoed ("<3) std::cout<<"GetNextGoodMuon: muon candidate after VK ("<3) std::cout<<"GetNextGoodMuon: muon candidate killed by pileup muon ("<second)->eventEntry; - idet = (it->second)->det_i; - idetID = (it->second)->det_id; - idetEdep = (it->second)->det_edep; - if (bool_debugingRequired) {if (debugEventMap[evtID]>1) std::cout<<"GetNextGoodMuon: GOOD muon candidate found (" - <::iterator itt = it_muon_hits_to_be_deleted.begin(); itt != it_muon_hits_to_be_deleted.end(); ++itt) { - hitMap.erase(*itt); - } - return true; - } - for(std::list::iterator itt = it_muon_hits_to_be_deleted.begin(); itt != it_muon_hits_to_be_deleted.end(); ++itt) { - hitMap.erase(*itt); - } - return false; -} - -//================================================================ -Bool_t musrCounter::CheckForPileupMuons(Long64_t timeBin0) { - // Check for pileup muons. If double hit in M-counter is found, return true. - Long64_t timeBinMinimum = timeBin0; - for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { - Long64_t timeBinOfCount_tmp = it->first; - // std::cout<<"timeBin0="< S T O P !!!\n"; exit(1);} - std::list it_positron_hits_to_be_deleted; - for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { - // Int_t eventNumber = (it->second)->eventIDnumber; - Long64_t timeBinOfCount_tmp = it->first; - Int_t modulo = timeBinOfCount_tmp % 524288; - if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) { - if (bool_debugingRequired) { - if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Hit out of data interval"< skip it - } - - // Hit candidate was found. Now check its coincidences and vetos - Bool_t bool_coincidenceConditions = true; - for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) { - if (bool_debugingRequired) { - if (debugEventMap[evtID]>4) { (itCounter->second)->myPrintThisCounter(evtID); } - } - if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') )) { - if (bool_debugingRequired) { - if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Coincidence required but not found (timeBin="< remove the candidate. - if (bool_ignoreUnperfectPositrons) it_positron_hits_to_be_deleted.push_back(it); - bool_coincidenceConditions = false; - goto CoincidencesChecked; - // goto endOfThisHit; // no coincidence found ==> skip hit - } - } - for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) { - if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') ) { - if (bool_debugingRequired) { - if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Coincidence vith veto detector found (timeBin="< remove the candidate. - if (bool_ignoreUnperfectPositrons) it_positron_hits_to_be_deleted.push_back(it); - bool_coincidenceConditions = false; - goto CoincidencesChecked; - // goto endOfThisHit; // coincidence with veto found ==> skip hit - } - } - - CoincidencesChecked: - if (!bool_coincidenceConditions) continue; // This hit does not fulfill coincidence and veto criteria - if (positronQuality==2) { - for(std::list::iterator itt = it_positron_hits_to_be_deleted.begin(); itt != it_positron_hits_to_be_deleted.end(); ++itt) { - hitMap.erase(*itt); - } - return 3; // An electron was already found before, and now again ==> double hit - } - else positronQuality=2; - - kEntry = (it->second)->eventEntry; - idet = (it->second)->det_i; - idetID = (it->second)->det_id; - idetEdep = (it->second)->det_edep; - timeBinOfNextGoodHit = timeBinOfCount_tmp; - timeBinOfNextGoodHit_phaseShifted = (it->second) -> timeBin2; - if (bool_debugingRequired) { - if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Good positron candidate found in this counter. (timeBin="<second)->eventIDnumber<<")"<::iterator itt = it_positron_hits_to_be_deleted.begin(); itt != it_positron_hits_to_be_deleted.end(); ++itt) { - hitMap.erase(*itt); - } - - return positronQuality; -} - -//================================================================ -void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max) { - // std::cout<<"QQQQQQQQQQQQQQQQQQQQQ koincidenceCounterMap.size()="<second)->GetMaxCoincidenceTimeWindow()); - if (maxCoinc < maxCoinc_AlreadySet) (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc); - - if (motherCounter=='M') (it->second)->SetCoincidenceTimeWindow_M(min,max); - else if (motherCounter=='P') (it->second)->SetCoincidenceTimeWindow_P(min,max); - else { - std::cout<<"musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors ERROR: Strange motherCounter " - < S T O P "<second; - Long64_t maxCoinc_AlreadySet = counter->GetMaxCoincidenceTimeWindow(); - Long64_t min_AlreadySet = counter->GetAntiCoincidenceTimeWindowMin(); - Long64_t max_AlreadySet = counter->GetAntiCoincidenceTimeWindowMax(); - if (maxCoinc < maxCoinc_AlreadySet) counter->SetMaxCoincidenceTimeWindow(maxCoinc); - if (min < min_AlreadySet) counter->SetAntiCoincidenceTimeWindowMin(min); - if (max > max_AlreadySet) counter->SetAntiCoincidenceTimeWindowMax(max); - } -} - -//================================================================ -void musrCounter::myPrintThisCounter(Int_t evtID, Int_t detail) { - Bool_t eventMixing=false; - if ((hitMap.begin()==hitMap.end()) && (detail<=1) ) return; - // if (detail>1) std::cout<<"musrCounter::myPrintThisCounter: counterNr = "<first<<" ("<<(it->first)%524288 <<") "<<(it->second)->eventIDnumber<<","; - if (evtID != (it->second)->eventIDnumber) {eventMixing=true;} - } - if (eventMixing) {std::cout<<" Potential event mixing";} - std::cout< (previousClock+musrAnalysis::clock_interval)) { -// previousClock += musrAnalysis::clock_interval; -// // dumpFile<<"-1\t"<=musrAnalysis::rewindTimeBins) tdcBin -= musrAnalysis::rewindTimeBins; - else if (tdcBin<0) tdcBin += musrAnalysis::rewindTimeBins; - // Int_t tdc = (tdcBinsend_to_dump(detID,tdcBin,false); -} -//================================================================ -void musrCounter::WriteRewindIntoDumpFile() { - DumpInfoToDumpFile(-2,1000,0); -} -//================================================================ +//#include +#include "musrCounter.hh" +#include "TCanvas.h" +#include "musrAnalysis.hh" + +typedef std::map debugEventMapType; +debugEventMapType debugEventMap; +Bool_t bool_debugingRequired; + +Bool_t musrCounter::bool_ignoreUnperfectMuons = true; +Bool_t musrCounter::bool_ignoreUnperfectPositrons = true; +Bool_t musrCounter::bool_WriteDataToDumpFile = false; +//Long64_t musrCounter::previousClock = -1; +//Long64_t musrCounter::CLOCK_INTERVAL = 512000; +//ofstream musrCounter::dumpFile; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +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; + counterNr = CHANNEL_NR; + strcpy(counterName,CHANNEL_NAME); + counterType = CHANNEL_TYPE; + couterEThreshold = E_THRESH; + counterTimeShift = (Long64_t) TIME_SHIFT; + keepBelowThresholdHits = KEEP_BELOW_THRESH; + discriminatorRecoveryTime = DEADTIME; + std::cout<<"musrCounter::musrCounter: Creating counter "<Fill(variable,vaha); +} + +//================================================================ +void musrCounter::DrawTDChistogram() { + char canvasName[501]; + sprintf(canvasName,"c%s",TDC_histoName); + TCanvas* cTmp = new TCanvas(canvasName,canvasName); + histTDC->Draw(); +} + +//================================================================ +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="<=couterEThreshold) || ((counterType=='P') && keepBelowThresholdHits)) { + // 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 "<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<first; + hitInfo* tempEvnr= it->second; + tempEvnr->RewindTimeBin2(musrAnalysis::rewindTimeBins); + hitMap_TMP.insert( std::pair(tempBinT-musrAnalysis::rewindTimeBins,tempEvnr) ); + } + hitMap.swap(hitMap_TMP); +} + +//================================================================ +Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, char motherCounter){ + // timeBin ... time bin, at which the coincidence is searched + + if (hitMap.empty()) return false; + Long64_t timeBinMin; + Long64_t timeBinMax; + + // If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors). + // Otherwise use timeBinMinimum and timeBinMaximum (e.g.coincidence of a positron counter with other positron counters). + if (counterType == 'V') timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector + else if (motherCounter=='M') timeBinMin = timeBin + coincidenceTimeWindowMin_M; // this is coinc. detector connected to M + else if (motherCounter=='P') timeBinMin = timeBin + coincidenceTimeWindowMin_P; // this is coinc. detector connected to P + + if (counterType == 'V') timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector + else if (motherCounter=='M') timeBinMax = timeBin + coincidenceTimeWindowMax_M; // this is coinc. detector connected to M + else if (motherCounter=='P') timeBinMax = timeBin + coincidenceTimeWindowMax_P; // this is coinc. detector connected to P + + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + Long64_t timeBinOfCount_tmp = it->first; + if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) { + // if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) { + return true; + // } + } + else if (timeBinOfCount_tmp > timeBinMax) return false; + } + return false; +} + +//================================================================ +Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep) { + // This function searches for a good muon, i.e. the muon + // 1) belongs to the currently analysed event + // 2) is in coincidence with all required coincidence detectors + // 3) is not in coincidence with veto detectors + // INPUT PARAMETERS: evtID, timeBinMin + // OUTPUT PARAMETERS: timeBinOfNextHit + // + // Loop over the hits in the counter + // std::cout<<" musrCounter::GetNextGoodMuon timeBinMin="< S T O P !!!\n"; exit(1);} + + std::list it_muon_hits_to_be_deleted; + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + + Long64_t timeBinOfCount_tmp = it->first; + timeBinOfNextHit = timeBinOfCount_tmp; + Int_t modulo = timeBinOfNextHit % 524288; + if (timeBinOfCount_tmp <= timeBinMin) continue; // This hit was already processed previously ==> skip it + + Int_t eventNumber = (it->second)->eventIDnumber; + if (eventNumber!=evtID) continue; // This trigger hit does not correspond to the currently processed event + // ==> skip it, because it was already proceesed or will be processed in future + numberOfMuonCandidates++; + if (bool_debugingRequired) {if (debugEventMap[evtID]>3) std::cout<<"GetNextGoodMuon: muon candidate found ("<second)->IsInCoincidence(timeBinOfCount_tmp,'M') )) { // no coincidence found ==> skip hit + // if (bool_ignoreUnperfectMuons) hitMap.erase(it); + if (bool_ignoreUnperfectMuons) it_muon_hits_to_be_deleted.push_back(it); + bool_coincidenceConditions = false; + if (bool_debugingRequired) {if (debugEventMap[evtID]>3) std::cout<<"GetNextGoodMuon: muon candidate not in koincidence ("<second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) { // coincidence with veto found ==> skip hit + // if (bool_ignoreUnperfectMuons) hitMap.erase(it); + if (bool_ignoreUnperfectMuons) it_muon_hits_to_be_deleted.push_back(it); + bool_coincidenceConditions = false; + if (bool_debugingRequired) {if (debugEventMap[evtID]>3) std::cout<<"GetNextGoodMuon: muon candidate vetoed ("<3) std::cout<<"GetNextGoodMuon: muon candidate after VK ("<3) std::cout<<"GetNextGoodMuon: muon candidate killed by pileup muon ("<second)->eventEntry; + idet = (it->second)->det_i; + idetID = (it->second)->det_id; + idetEdep = (it->second)->det_edep; + if (bool_debugingRequired) {if (debugEventMap[evtID]>1) std::cout<<"GetNextGoodMuon: GOOD muon candidate found (" + <::iterator itt = it_muon_hits_to_be_deleted.begin(); itt != it_muon_hits_to_be_deleted.end(); ++itt) { + hitMap.erase(*itt); + } + return true; + } + for(std::list::iterator itt = it_muon_hits_to_be_deleted.begin(); itt != it_muon_hits_to_be_deleted.end(); ++itt) { + hitMap.erase(*itt); + } + return false; +} + +//================================================================ +Bool_t musrCounter::CheckForPileupMuons(Long64_t timeBin0) { + // Check for pileup muons. If double hit in M-counter is found, return true. + Long64_t timeBinMinimum = timeBin0; + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + Long64_t timeBinOfCount_tmp = it->first; + // std::cout<<"timeBin0="< S T O P !!!\n"; exit(1);} + std::list it_positron_hits_to_be_deleted; + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + // Int_t eventNumber = (it->second)->eventIDnumber; + Long64_t timeBinOfCount_tmp = it->first; + Int_t modulo = timeBinOfCount_tmp % 524288; + if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) { + if (bool_debugingRequired) { + if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Hit out of data interval"< skip it + } + + // Hit candidate was found. Now check its coincidences and vetos + Bool_t bool_coincidenceConditions = true; + for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) { + if (bool_debugingRequired) { + if (debugEventMap[evtID]>4) { (itCounter->second)->myPrintThisCounter(evtID); } + } + if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') )) { + if (bool_debugingRequired) { + if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Coincidence required but not found (timeBin="< remove the candidate. + if (bool_ignoreUnperfectPositrons) it_positron_hits_to_be_deleted.push_back(it); + bool_coincidenceConditions = false; + goto CoincidencesChecked; + // goto endOfThisHit; // no coincidence found ==> skip hit + } + } + for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) { + if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') ) { + if (bool_debugingRequired) { + if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Coincidence vith veto detector found (timeBin="< remove the candidate. + if (bool_ignoreUnperfectPositrons) it_positron_hits_to_be_deleted.push_back(it); + bool_coincidenceConditions = false; + goto CoincidencesChecked; + // goto endOfThisHit; // coincidence with veto found ==> skip hit + } + } + + CoincidencesChecked: + if (!bool_coincidenceConditions) continue; // This hit does not fulfill coincidence and veto criteria + if (positronQuality==2) { + for(std::list::iterator itt = it_positron_hits_to_be_deleted.begin(); itt != it_positron_hits_to_be_deleted.end(); ++itt) { + hitMap.erase(*itt); + } + return 3; // An electron was already found before, and now again ==> double hit + } + else positronQuality=2; + + kEntry = (it->second)->eventEntry; + idet = (it->second)->det_i; + idetID = (it->second)->det_id; + idetEdep = (it->second)->det_edep; + timeBinOfNextGoodHit = timeBinOfCount_tmp; + timeBinOfNextGoodHit_phaseShifted = (it->second) -> timeBin2; + if (bool_debugingRequired) { + if (debugEventMap[evtID]>3) {std::cout<<"GetNextGoodPositron: Good positron candidate found in this counter. (timeBin="<second)->eventIDnumber<<")"<::iterator itt = it_positron_hits_to_be_deleted.begin(); itt != it_positron_hits_to_be_deleted.end(); ++itt) { + hitMap.erase(*itt); + } + + 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) { + // std::cout<<"QQQQQQQQQQQQQQQQQQQQQ koincidenceCounterMap.size()="<second)->GetMaxCoincidenceTimeWindow()); + if (maxCoinc < maxCoinc_AlreadySet) (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc); + + if (motherCounter=='M') (it->second)->SetCoincidenceTimeWindow_M(min,max); + else if (motherCounter=='P') (it->second)->SetCoincidenceTimeWindow_P(min,max); + else { + std::cout<<"musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors ERROR: Strange motherCounter " + < S T O P "<second; + Long64_t maxCoinc_AlreadySet = counter->GetMaxCoincidenceTimeWindow(); + Long64_t min_AlreadySet = counter->GetAntiCoincidenceTimeWindowMin(); + Long64_t max_AlreadySet = counter->GetAntiCoincidenceTimeWindowMax(); + if (maxCoinc < maxCoinc_AlreadySet) counter->SetMaxCoincidenceTimeWindow(maxCoinc); + if (min < min_AlreadySet) counter->SetAntiCoincidenceTimeWindowMin(min); + if (max > max_AlreadySet) counter->SetAntiCoincidenceTimeWindowMax(max); + } +} + +//================================================================ +void musrCounter::myPrintThisCounter(Int_t evtID, Int_t detail) { + Bool_t eventMixing=false; + if ((hitMap.begin()==hitMap.end()) && (detail<=1) ) return; + // if (detail>1) std::cout<<"musrCounter::myPrintThisCounter: counterNr = "<first<<" ("<<(it->first)%524288 <<") "<<(it->second)->eventIDnumber<<","; + if (evtID != (it->second)->eventIDnumber) {eventMixing=true;} + } + if (eventMixing) {std::cout<<" Potential event mixing";} + std::cout< (previousClock+musrAnalysis::clock_interval)) { +// previousClock += musrAnalysis::clock_interval; +// // dumpFile<<"-1\t"<=musrAnalysis::rewindTimeBins) tdcBin -= musrAnalysis::rewindTimeBins; + else if (tdcBin<0) tdcBin += musrAnalysis::rewindTimeBins; + // Int_t tdc = (tdcBinsend_to_dump(detID,tdcBin,false); +} +//================================================================ +void musrCounter::WriteRewindIntoDumpFile() { + DumpInfoToDumpFile(-2,1000,0); +} +//================================================================ diff --git a/musrSimAna/musrCounter.hh b/musrSimAna/musrCounter.hh index 0fca77d..b4df8cd 100644 --- a/musrSimAna/musrCounter.hh +++ b/musrSimAna/musrCounter.hh @@ -1,112 +1,123 @@ -#ifndef musrCounter_h -#define musrCounter_h 1 -#include -#include -#include -#include -#include -#include -#include -#include - -class hitInfo { -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() {} - void RewindTimeBin2(Long64_t timeBinsToRewind) {timeBin2-=timeBinsToRewind;} - Int_t eventEntry; - Int_t eventIDnumber; - Int_t det_i; - Int_t det_id; - Double_t det_edep; - Long64_t timeBin2; - - // extern double GlobalKamil; - // typedef std::map debugEventMapType; - // extern debugEventMapType debugEventMap; - // extern Bool_t bool_debugingRequired; -}; - -class musrCounter { - public: - musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT); - ~musrCounter(); - int GetCounterNr() {return counterNr;} - char GetCounterType() {return counterType;} - void SetCoincidenceCounter(musrCounter* c, int icNr) { - int cNr = abs(icNr); - if (icNr>0) { - std::cout<<"SetCoincidenceCounter: Adding counter ="< counterMapType; - counterMapType koincidenceCounterMap; - counterMapType vetoCounterMap; - char TDC_histoName[200]; - int TDC_t0, TDC_t1, TDC_t2, TDC_histoNrAdd; - char TDC_histoNameAdd[200]; - TH1D* histTDC; - Long64_t antiCoincidenceTimeWindowMin, coincidenceTimeWindowMin_M, coincidenceTimeWindowMin_P; - Long64_t antiCoincidenceTimeWindowMax, coincidenceTimeWindowMax_M, coincidenceTimeWindowMax_P; - Long64_t maxCoincidenceTimeWindow; - // typedef std::map hitMap_TYPE; // Long64_t = timeBin, Int_t=eventID - typedef std::map hitMap_TYPE; // Long64_t = timeBin, hitInfo = eventID and det_i - hitMap_TYPE hitMap; - // std::list timeOfHitsList; - - Int_t doubleHitN; - Long64_t numberOfMuonCandidates; - Long64_t numberOfMuonCandidatesAfterVK; - Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval; - - public: - static Bool_t bool_ignoreUnperfectMuons; - static Bool_t bool_ignoreUnperfectPositrons; - static Bool_t bool_WriteDataToDumpFile; - // static ofstream dumpFile; - // static Long64_t previousClock; - // static Long64_t CLOCK_INTERVAL; -}; - -#endif +#ifndef musrCounter_h +#define musrCounter_h 1 +#include +#include +#include +#include +#include +#include +#include +#include + +class hitInfo { +public: + 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() {} + void RewindTimeBin2(Long64_t timeBinsToRewind) {timeBin2-=timeBinsToRewind;} + Int_t eventEntry; + Int_t eventIDnumber; + Int_t det_i; + Int_t det_id; + Double_t det_edep; + Long64_t timeBin1; // repeated here in case this is not stored in a Map with First=timeBin1 + Long64_t timeBin2; + Int_t multiHitCtr; + hitInfo* firstMulti; + Int_t posQual; + + // extern double GlobalKamil; + // typedef std::map debugEventMapType; + // extern debugEventMapType debugEventMap; + // extern Bool_t bool_debugingRequired; +}; + +class musrCounter { + public: + 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(); + int GetCounterNr() {return counterNr;} + char GetCounterType() {return counterType;} + void SetCoincidenceCounter(musrCounter* c, int icNr) { + int cNr = abs(icNr); + if (icNr>0) { + std::cout<<"SetCoincidenceCounter: Adding counter ="< counterMapType; + counterMapType koincidenceCounterMap; + counterMapType vetoCounterMap; + char TDC_histoName[200]; + int TDC_t0, TDC_t1, TDC_t2, TDC_histoNrAdd; + char TDC_histoNameAdd[200]; + TH1D* histTDC; + Long64_t antiCoincidenceTimeWindowMin, coincidenceTimeWindowMin_M, coincidenceTimeWindowMin_P; + Long64_t antiCoincidenceTimeWindowMax, coincidenceTimeWindowMax_M, coincidenceTimeWindowMax_P; + Long64_t maxCoincidenceTimeWindow; + // typedef std::map hitMap_TYPE; // Long64_t = timeBin, Int_t=eventID + typedef std::map hitMap_TYPE; // Long64_t = timeBin, hitInfo = eventID and det_i + hitMap_TYPE hitMap; + // std::list timeOfHitsList; + // for pulsed double hit processing + double discriminatorLastVolts; + Long64_t discriminatorLastTime; + hitMap_TYPE::iterator discriminatorIterator; + Double_t discriminatorRecoveryTime; // in TDC bins! + + Int_t doubleHitN; + Long64_t numberOfMuonCandidates; + Long64_t numberOfMuonCandidatesAfterVK; + Long64_t numberOfMuonCandidatesAfterVKandDoubleHitRemoval; + + public: + static Bool_t bool_ignoreUnperfectMuons; + static Bool_t bool_ignoreUnperfectPositrons; + static Bool_t bool_WriteDataToDumpFile; + // static ofstream dumpFile; + // static Long64_t previousClock; + // static Long64_t CLOCK_INTERVAL; +}; + +#endif diff --git a/musrSimAna/musrTH.cxx b/musrSimAna/musrTH.cxx index 8103c4c..ce0410d 100644 --- a/musrSimAna/musrTH.cxx +++ b/musrSimAna/musrTH.cxx @@ -49,7 +49,7 @@ void musrTH::FillTH1D(Double_t vaha, Bool_t* cond){ for (Int_t i=0; iGetName()<<"\", funct_xMin="<FixParameter(0,omega);} diff --git a/src/meyer.cc b/src/meyer.cc index 2abd560..c3567f4 100644 --- a/src/meyer.cc +++ b/src/meyer.cc @@ -60,6 +60,10 @@ #include #include #include +// James Lord: +// already defines its own constant Pi so no need to look for M_PI anywhere +// #define _USE_MATH_DEFINES 1 +// #include using namespace std; meyer::meyer() @@ -258,7 +262,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2, } // std::cout<< "M4: "<50.) { 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 // ! 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: @@ -305,7 +309,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2, } // ! 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; 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; 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 @@ -348,7 +352,7 @@ void meyer::Get_F_Function_Meyer(double tau, double Ekin, double Z1, double Z2, //! Berechne die Werte fuer theta=0: 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; fValuesFolded[1] = 0.; @@ -463,7 +467,7 @@ void meyer:: Get_F_Function(double tau,double theta, double Ekin, double Z1, dou nBin = 0; - thetaSchlange = Meyer_faktor4 * Ekin * theta *M_PI/180; + thetaSchlange = Meyer_faktor4 * Ekin * theta *Pi/180; 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; } - *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++) { @@ -486,7 +490,7 @@ void meyer:: Get_F_Function(double tau,double theta, double Ekin, double Z1, dou //! Berechne die Werte fuer theta=0: double norm; 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; } diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 5bd84b2..4c67088 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -44,6 +44,7 @@ #include "G4Box.hh" #include "G4Cons.hh" #include "G4Polycone.hh" +#include "G4Torus.hh" #include "G4LogicalVolume.hh" #include "G4PVPlacement.hh" #include "G4Tubs.hh" @@ -329,6 +330,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { solidName+=name; 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 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); @@ -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); solidName+=name; 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(G4ThreeVector(0,0,1),45*deg); // G4RotationMatrix* rot = new G4RotationMatrix(); @@ -673,6 +681,28 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { solid = new G4SubtractionSolid(solidName,solidGPDmHolder,solidGPDcutOut,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); G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm); @@ -1631,6 +1661,41 @@ void musrDetectorConstruction::DefineMaterials() Air->AddElement(N, fractionmass=0.7); 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 // @@ -1661,6 +1726,9 @@ void musrDetectorConstruction::DefineMaterials() 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("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) { diff --git a/src/musrErrorMessage.cc b/src/musrErrorMessage.cc index e220369..60286be 100644 --- a/src/musrErrorMessage.cc +++ b/src/musrErrorMessage.cc @@ -47,6 +47,7 @@ void musrErrorMessage::musrError(SEVERITY severity, G4String message, G4bool sil actualErrorMessage.mesSeverity = severity; actualErrorMessage.nTimes = 1; ErrorMapping[message]=actualErrorMessage; + it = ErrorMapping.find(message); G4cout<<"!!!"<AbortRun(true); } - if (thisEventNr != 0 and thisEventNr%nHowOftenToPrintEvent == 0) { + if ((thisEventNr != 0) && (thisEventNr%nHowOftenToPrintEvent == 0)) { time_t curr=time(0); //char * ctime(const time_t * tp); G4cout << ">>> Event " << evt->GetEventID() <<". Running already for "<GetTrajectoryContainer()))[i]); - trj->DrawTrajectory(1000); + trj->DrawTrajectory(); //removed argument 1000 (JSL) } } } diff --git a/src/musrPhysicsList.cc b/src/musrPhysicsList.cc index aa40a46..2f05521 100644 --- a/src/musrPhysicsList.cc +++ b/src/musrPhysicsList.cc @@ -264,6 +264,7 @@ void musrPhysicsList::ConstructProcess() // For optical photons #include "G4Scintillation.hh" +#include "G4Cerenkov.hh" // For models #include "G4ProcessTable.hh" @@ -336,6 +337,13 @@ void musrPhysicsList::ConstructEM() else if (stringProcessName=="G4OpRayleigh") pManager->AddDiscreteProcess(new G4OpRayleigh); else if (stringProcessName=="G4OpBoundaryProcess") pManager->AddDiscreteProcess(new G4OpBoundaryProcess); 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 { sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess. It can be easily added.", charProcessName); @@ -373,6 +381,14 @@ void musrPhysicsList::ConstructEM() if (musrParameters::finiteRiseTimeInScintillator) myScintillation->SetFiniteRiseTime(true); 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 // in his original "musrPhysicsList.cc about implementing musrMuScatter. // else if (stringProcessName=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3); diff --git a/src/musrRootOutput.cc b/src/musrRootOutput.cc index f25879d..23ea2ba 100644 --- a/src/musrRootOutput.cc +++ b/src/musrRootOutput.cc @@ -833,8 +833,9 @@ void musrRootOutput::SetPhotDetTime(G4double time) { nOptPhotDet++; } else { + nOptPhotDet++; // still want to know how many in total (JSL) char message[200]; 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) } } diff --git a/src/musrScintSD.cc b/src/musrScintSD.cc index 00b621d..77e384a 100644 --- a/src/musrScintSD.cc +++ b/src/musrScintSD.cc @@ -35,6 +35,7 @@ #include "musrErrorMessage.hh" #include "musrSteppingAction.hh" //#include "TCanvas.h" +#include "TH1D.h" #include "TMath.h" #include "TF1.h" #include @@ -131,6 +132,7 @@ musrScintSD::musrScintSD(G4String name) APDcellsTimeVariationRequested=false; APDcrossTalkRequested=false; APDcrossTalkProb=0.; + // verboseLevel = 9; // JSL } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -167,8 +169,10 @@ void musrScintSD::Initialize(G4HCofThisEvent* HCE) { if (musrParameters::boolG4OpticalPhotons) { if (!musrParameters::boolG4OpticalPhotonsUnprocess) { for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) { + if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting optHitDetectorMap" <<"\n"; delete (it->second); } + if (verboseLevel>1) G4cout<<"VERBOSE 2 : clearing optHitMap" <<"\n"; 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(); // newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName()); newHit->SetParticleName(particleName); @@ -286,6 +291,7 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) { // G4cout<<" detID ="<(detID,optHitDetectorMapTmp)); iter = optHitMap.find(detID); @@ -380,7 +386,7 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) { // in which case times are huge and due to the limited rounding // precision become equal --> map ignores the same "keys", // multimap does not. - std::map::iterator it; + std::multimap::iterator it; for (G4int i=0; iGetGlobalTime(); @@ -670,6 +676,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() { 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); + if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new TH1D for OPSAhisto" <<"\n"; OPSAhisto = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); // poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2 // poiss->SetParameter(0,1); @@ -677,9 +684,11 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() { if (bool_pulseShapeExists) { sprintf(nameHist,"OPSAshape_%d_%d_%d",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); 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); + if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating a new TH1D for OPSA_CFD" <<"\n"; 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 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); + if (verboseLevel>1) G4cout<<"VERBOSE 2 : creating new TH1Ds for OPSAhistoSUM,0" <<"\n"; OPSAhistoSUM = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); OPSAhistoSUM0= new TH1D(nameHist0,nameHistTitle0,OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax); mapOfOPSAsumHistograms.insert(std::pair(nameHist,OPSAhistoSUM)); @@ -1026,13 +1036,16 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() { // Delete the histograms from memory if they were created if ((boolStoreThisOPSAhist)||(bool_pulseShapeExists)) { + if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting OPSAhisto" <<"\n"; delete OPSAhisto; if (bool_pulseShapeExists) { + if (verboseLevel>1) G4cout<<"VERBOSE 2 : deleting OPSAshape and CFD" <<"\n"; delete OPSAshape; 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, OPSA_timeA,OPSA_timeB,OPSA_timeC,OPSA_timeD,OPSA_timeMean,OPSA_timeLast, 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 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); } OPSA_signal_Map.clear(); @@ -1162,16 +1176,20 @@ G4int musrScintSD::FindAPDcellID(G4Step* aStep) { void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time, G4int nTruePhe) { if ( (musrRootOutput::store_phot_time) && (nTruePhe>0) ) myRootOutput->SetPhotDetTime(time); if (!APDcellsEffectRequested) { + if (verboseLevel>1) G4cout<<"VERBOSE 2 : storing photon hit (non APD)" <<"\n"; APDidClass tmpAPDid(0,nTruePhe); // optHitDetectorMap->insert(std::pair(time,0)); optHitDetectorMap->insert(std::pair(time,tmpAPDid)); } else { 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++ ) { 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; 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; APDidClass tmpAPDid(APDcellID,tmpAPDcellNphot); optHitDetectorMap->erase(it2); @@ -1183,6 +1201,7 @@ void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int AP } if (!APDcellAlreadyFired) { // optHitDetectorMap->insert(std::pair(time,APDcellID)); + if (verboseLevel>1) G4cout<<"VERBOSE 2 : this cell hasn't fired, storing" <<"\n"; APDidClass tmpAPDid(APDcellID,nTruePhe); optHitDetectorMap->insert(std::pair(time,tmpAPDid)); } diff --git a/src/musrTabulatedElementField.cc b/src/musrTabulatedElementField.cc index 47ec04c..49a1f6e 100644 --- a/src/musrTabulatedElementField.cc +++ b/src/musrTabulatedElementField.cc @@ -617,10 +617,12 @@ void musrTabulatedElementField::addFieldValue2D(const G4double point[4], local = global2local.TransformPoint(global); double x, z, z_sign; - if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)|| - (strcmp(fieldTableType,"2D_OperaXY"))||(strcmp(fieldTableType,"2DEf")==0)) { +// if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==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" // 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()); z = fabs(local.z()); z_sign = (local.z()>0) ? 1.:-1.; @@ -631,8 +633,9 @@ void musrTabulatedElementField::addFieldValue2D(const G4double point[4], z = local.z(); z_sign = 1; } - // Check that the point is within the defined region - if ( x=minimumx && x=minimumz && zevNrKriz) std::cout<<"bol som tu"<