start stm study report

stm_study_201611/Makefile

@@ -0,0 +1,15 @@
+DOC=stm_study
+INPUT=$(DOC).tex
+TARGET=$(DOC).pdf
+REF=ref.bib
+TEX=pdflatex -shell-escape
+BIB=bibtex
+
+default: $(TARGET)
+
+$(TARGET): $(INPUT) Makefile 
+	$(TEX) $< && $(BIB) $(DOC) && $(TEX) $< && $(TEX) $<
+
+clean:
+	rm -f $(DOC).{pdf,out,aux,log}
+	rm -f *.{pdf,out,aux,log,bbl,blg}

stm_study_201611/listings/G4EmCaptureCascade.cc

@@ -0,0 +1,181 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4EmCaptureCascade.cc 101422 2016-11-17 10:41:23Z gcosmo $
+//
+//-----------------------------------------------------------------------------
+//
+// GEANT4 Class file 
+//
+// File name:  G4EmCaptureCascade
+//
+// Author:        V.Ivanchenko (Vladimir.Ivantchenko@cern.ch)
+// 
+// Creation date: 22 April 2012 on base of G4MuMinusCaptureCascade
+//
+//
+//-----------------------------------------------------------------------------
+//
+// Modifications: 
+//
+//-----------------------------------------------------------------------------
+
+#include "G4EmCaptureCascade.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "Randomize.hh" 
+#include "G4MuonMinus.hh"
+#include "G4Electron.hh"
+#include "G4Gamma.hh"
+#include "G4NucleiProperties.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4EmCaptureCascade::G4EmCaptureCascade()
+  : G4HadronicInteraction("emCaptureCascade")
+{ 
+  theElectron = G4Electron::Electron();
+  theGamma = G4Gamma::Gamma();
+  fMuMass = G4MuonMinus::MuonMinus()->GetPDGMass();
+  fTime = 0.0;
+
+  // Calculate the Energy of K Mesoatom Level for this Element using
+  // the Energy of Hydrogen Atom taken into account finite size of the
+  // nucleus 
+  static const G4int nlevels = 28;
+  static const G4int listK[nlevels] = {
+      1, 2,  4,  6,  8, 11, 14, 17, 18, 21, 24,
+     26, 29, 32, 38, 40, 41, 44, 49, 53, 55,
+     60, 65, 70, 75, 81, 85, 92};
+  static const G4double listKEnergy[nlevels] = {
+     0.00275, 0.011, 0.043, 0.098, 0.173, 0.326,
+     0.524, 0.765, 0.853, 1.146, 1.472,
+     1.708, 2.081, 2.475, 3.323, 3.627, 
+     3.779, 4.237, 5.016, 5.647, 5.966,
+     6.793, 7.602, 8.421, 9.249, 10.222,
+    10.923,11.984};
+
+  fKLevelEnergy[0] = 0.0;
+  fKLevelEnergy[1] = listKEnergy[0];
+  G4int idx = 1;
+  for(G4int i=1; i<nlevels; ++i) {
+    G4int z1 = listK[idx];
+    G4int z2 = listK[i];
+    if(z1+1 < z2) {
+      G4double dz = G4double(z2 - z1);
+      G4double y1 = listKEnergy[idx]/G4double(z1*z1);
+      G4double y2 = listKEnergy[i]/G4double(z2*z2);
+      for(G4int z=z1+1; z<z2; ++z) {
+        fKLevelEnergy[z] = (y1 + (y2 - y1)*(z - z1)/dz)*z*z;
+      }
+    }
+    fKLevelEnergy[z2] = listKEnergy[i];
+    idx = i;  
+  }
+  for(G4int i = 0; i<14; ++i) { fLevelEnergy[i] = 0.0; }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4EmCaptureCascade::~G4EmCaptureCascade()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4HadFinalState* 
+G4EmCaptureCascade::ApplyYourself(const G4HadProjectile& projectile, 
+				  G4Nucleus& targetNucleus)
+{
+  result.Clear();
+  result.SetStatusChange(isAlive);
+  fTime = projectile.GetGlobalTime();
+
+  G4int Z = targetNucleus.GetZ_asInt(); 
+  G4int A = targetNucleus.GetA_asInt(); 
+  G4double massA = G4NucleiProperties::GetNuclearMass(A, Z);
+  G4double mass = fMuMass * massA / (fMuMass + massA) ;
+  G4double e = 13.6 * eV * (Z * Z) * mass/ electron_mass_c2;
+
+  // precise corrections of energy only for K-shell
+  fLevelEnergy[0] = fKLevelEnergy[std::min(Z, 92)];
+  for(G4int i=1; i<14; ++i) {
+    fLevelEnergy[i] = e/(G4double)((i+1)*(i+1));
+  }
+
+  G4int nElec  = Z;
+  G4int nAuger = 1;
+  G4int nLevel = 13;
+  G4double pGamma = (Z*Z*Z*Z);
+
+  // Capture on 14-th level
+  G4double edep = fLevelEnergy[13];
+  AddNewParticle(theElectron,edep);
+  G4double deltaE;
+
+  // Emit new photon or electron
+  // Simplified model for probabilities
+  // N.C.Mukhopadhyay Phy. Rep. 30 (1977) 1.
+  do {
+
+    // case of Auger electrons
+    if((nAuger < nElec) && ((pGamma + 10000.0) * G4UniformRand() < 10000.0) ) {
+      ++nAuger;
+      deltaE =  fLevelEnergy[nLevel-1] -  fLevelEnergy[nLevel];
+      --nLevel;
+      AddNewParticle(theElectron, deltaE);
+
+    } else {
+
+      // Case of photon cascade, probabilities from
+      // C.S.Wu and L.Wilets, Ann. Rev. Nuclear Sci. 19 (1969) 527.
+
+      G4double var = (10.0 + G4double(nLevel - 1) ) * G4UniformRand();
+      G4int iLevel = nLevel - 1 ;
+      if(var > 10.0) iLevel -= G4int(var-10.0) + 1;
+      if( iLevel < 0 ) iLevel = 0;
+      deltaE =  fLevelEnergy[iLevel] -  fLevelEnergy[nLevel];
+      nLevel = iLevel;
+      AddNewParticle(theGamma, deltaE);
+    }
+    edep += deltaE;
+
+    // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
+  } while( nLevel > 0 );
+
+  result.SetLocalEnergyDeposit(edep);
+  return &result;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4EmCaptureCascade::ModelDescription(std::ostream& outFile) const
+{
+  outFile << "Simulation of electromagnetic cascade from capture level"
+	  << " to K-shell of the mesonic atom\n."
+	  << "Probabilities of gamma and Auger transitions from\n"
+	  << "  N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1.\n";
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

stm_study_201611/listings/code_dir_tree.sh

@@ -0,0 +1,10 @@
+@mu2egpvm02:STM_study_201611 namtran $ pwd
+/mu2e/app/users/namtran/STM_study_201611
+@mu2egpvm02:STM_study_201611 namtran $ tree -L 1
+.
+|-- analysis
+|-- script
+|-- step00
+`-- step01
+
+4 directories, 0 files

stm_study_201611/stm_study.tex

@@ -0,0 +1,229 @@
+\documentclass[11pt]{article}
+\usepackage{mhchem}
+\usepackage{booktabs}
+\usepackage{multirow}
+\usepackage{textcomp}
+\usepackage{epsfig}
+\usepackage{hyperref}
+\usepackage[noabbrev, capitalize]{cleveref} % hyperref must be loaded first
+\usepackage[
+  detect-weight=true,
+  per=slash,
+  detect-family=true,
+  separate-uncertainty=true]{siunitx}
+% \usepackage{listings}
+\usepackage{xcolor}
+\usepackage{upquote}
+% \usepackage{minted}
+\usepackage{tcolorbox}
+\tcbuselibrary{minted,skins}
+\usemintedstyle{monokai}
+
+\definecolor{greybg}{rgb}{0.25,0.25,0.25}
+\newtcblisting{bashcode}{
+  listing engine=minted,
+  colback=bashcodebg,
+  colframe=black!70,
+  listing only,
+  minted style=colorful,
+  minted language=bash,
+  minted options={linenos=true,texcl=true},
+  left=1mm,
+}
+
+% \DeclareSIUnit\eVperc{\eV\per\clight}
+% \DeclareSIUnit\clight{\text{\ensuremath{c}}}
+
+\begin{document}
+\title{Hit rate estimation for STM detectors}
+\author{Nam H. Tran \\ Boston University}
+\date{\today}
+\maketitle
+
+\begin{abstract}
+  This report presents estimated hit rate on STM detectors with the updated
+  geometry in Offline version TODO.
+\end{abstract}
+\section{Overview}
+\label{sec:overview}
+In order to measure the eponymous rate of the Mu2e experiment, the number of
+stopped muons needs to be known to about \SI{10}{\percent}. The most promising
+schemes have involved photon detectors far downstream of the Muon Beam Stop
+measuring emissions from the Stopping Target (ST) at different times with
+respect to each pulse of muons stopping in the target. There are three
+categories we explore here for this stopping-target monitor (STM):
+\begin{itemize}
+  \item Prompt: X-rays emitted when the muon comes to rest in the ST and is
+    captured into atomic orbit (atomic capture, sometimes referred to as “muon
+    stop” in this text). These X-rays are emitted essentially O(ps) with the
+    atomic capture time.
+  \item Semiprompt: $\gamma$-rays emitted upon nuclear capture (sometimes
+    referred to simply as ``capture'' in this text) of the muon after stop.
+    These exhibit timing characteristics of the muonic lifetime ($\tau
+    = \SI{864}{\ns}$).
+  \item Delayed: $\gamma$-rays from activated daughters resulting from muon
+    nuclear capture in the ST.
+\end{itemize}
+The baseline design of the STM has been described by Miller et
+al~\cite{Miller2016}, the estimated hit rate on STM detector was about
+\SI{1.1}{\kHz} in the time window \SIrange{200}{1695}{\ns} from arrival time of
+a proton bunch. But this design would cause a high hit rate on the Cosmic
+Ray Veto (CRV) system, therefore changes have been made in the Offline version
+v6\textunderscore 0\textunderscore 2 by the CRV team to reduce the CRV hit
+rate. The goal of this study is to re-estimate hit rates on STM detectors in
+this new geometry.
+
+The actual changes relevant to the STM are:
+\begin{itemize}
+  \item \SI{10}{\mm} polyethylene liner added to CRV shielding,
+  \item field-of-view (FOV) collimator absorber thickness reduced to
+    \SI{10}{\mm} (from \SI{20}{\mm}).
+\end{itemize}
+
+\section{Simulation details}
+\label{sec:simulation_details}
+The study was done using Mu2e Offline version v6\textunderscore
+0\textunderscore 2 (released on Oct 16,
+2016), hashtag \texttt{3d1e9154d7}. The simulation starts from the entrance of
+TS5 (see \cref{fig:stm_geo_all}), taking \texttt{cd3-beam-g4s2-mubeam.0728a}
+dataset as input. The dataset contains 5098 files, each corresponds to
+\num{1e6} proton-on-target (POT). The dataset were reused 16 times with
+different random seeds, where \SI{97}{\percent} of runs succeeded, equivalent
+to \num{8e11} POTs.
+
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=1.0\textwidth]{figs/stm_geo_all}
+  \caption{Simulation geometry showing the DS region on the left, sweeper magnet,
+    FOV collimator, spot-size collimator, and the STM detectors on the right.
+    Particles saved in the input files are shoot from the TS5 (orange circle),
+    and transported to the STM region.}
+  \label{fig:stm_geo_all}
+\end{figure}
+
+There were 8 virtual detectors (VD) in STM region enabled in this study, their
+identification numbers (\texttt{vdid}), locations, and abbreviation names
+(appear in the simulation output) are listed in \cref{tab:vds_list}.
+Information recorded by the VDs includes: particle type (\texttt{pdgid}),
+global and local coordinates, time, kinetic energy, and parent particle type.
+Only particles considered important to the STM, namely electrons,
+positrons, negative and positive muons, neutrons and photons, were written to
+the output file.
+
+\begin{table}[htbp]
+\centering
+\caption{List of virtual detectors read out in this study}
+\label{tab:vds_list}
+\begin{tabular}{@{}ccll@{}}
+\toprule
+  &VDID & Location               & Abbreviation              \\
+\midrule
+1 & 81  & Exit of neutron shield of the DS      & DSNeutronShieldExit       \\
+2 & 86  & Upstream of the STM system            & STM\textunderscore UpStr                \\
+3 & 87  & Downstream of the sweeper magnet      & STM\textunderscore MagDnStr             \\
+4 & 101 & Upstream of the spot-size collimator  & STM\textunderscore SpotSizeCollUpStr    \\
+5 & 88  & Downstream of the spot-size collimator& STM\textunderscore CollDnStr            \\
+6 & 89  & Upstream of the STM detector 1        & STM\textunderscore Det1UpStr            \\
+7 & 90  & Upstream of the STM detector 2        & STM\textunderscore Det2UpStr            \\
+8 & 100 & Downstream of the FOV collimator      & STM\textunderscore FieldOfViewCollDnStr \\
+\bottomrule
+\end{tabular}
+\end{table}
+
+\section{Simulation and analysis code}
+\label{sec:simulation_and_analysis_code}
+
+The simulation and analysis code are located at:
+\url{/mu2e/app/users/namtran/STM_study_201611}.
+
+% \lstinputlisting[language=bash,frame=single]{listings/code_dir_tree.sh}
+% \inputminted{bash}{listings/code_dir_tree.sh}
+\tcbinputlisting{%
+  listing file=listings/code_dir_tree.sh, colback=red!5!white,
+  colframe=red!25, left=6mm,
+  minted options={style=tcblatex, numbers=left, numberstyle=\tiny\color{red!75!black}}
+}
+
+\texttt{step00} contains configuration files for this simulation and a script to
+submit all 5098 jobs (correspond to number of input files) to the FermiGrid.
+It took about 14 hours to complete a job in average.
+
+The \texttt{analysis} folder contains a script 
+(\texttt{run\textunderscore statistics.sh}) which checks if a job has finished
+successfully, and makes a list of such runs.
+There is a simple analysis code (\texttt{main.cc}) to read the VD records and
+make plots.
+
+\section{Results}
+\label{sec:results}
+\subsection{STM detector spectra}
+\label{sub:stm_detector_spectra}
+
+Energy spectrum of particles hitting STM detectors are presented in
+\cref{fig:stm_det_ke}. There were not many hits, and only the annihilation
+peak stands out. Most of the particles are photons as shown in
+\cref{fig:stm_det_ptype}.
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=0.7\textwidth]{figs/ke_det1UpStr}
+  \includegraphics[width=0.7\textwidth]{figs/ke_det2UpStr}
+  \caption{Kinetic energy of particles hitting STM detectors 1 (top), and
+    2 (bottom).}
+  \label{fig:stm_det_ke}
+\end{figure}
+\begin{figure}[htbp]
+  \centering
+  \includegraphics[width=0.7\textwidth]{figs/ke_pdg_det1UpStr}
+  \includegraphics[width=0.7\textwidth]{figs/ke_pdg_det2UpStr}
+  \caption{Kinetic energy and type of particles hitting STM detectors 1 (top),
+    and 2 (bottom).}
+  \label{fig:stm_det_ptype}
+\end{figure}
+
+\subsection{STM detector hit rate estimation}
+\label{sub:stm_detector_hit_rate_estimation}
+The average number of hits on a STM detector per POT is:
+\begin{equation}
+  \frac{888 + 888}{2 \times 8 \times 10^{11}} = 8.7 \times 10^{-9}. 
+\end{equation}
+There are 3.1 POTs per proton bunch, so the number of hits per bunch is:
+\begin{equation}
+  8.7 \times 10^{-9} \times 3.1 \times 10^7 = 0.27
+\end{equation}
+The instantaneous hit rate, assuming an interval of \SI{1695}{\ns} between
+bunches, is:
+\begin{equation}
+  \frac{0.27}{1695\times 10^{-9}} = \SI{159e3}{\Hz}
+\end{equation}
+
+\appendix
+\section{How to run the simulation and analyze data}
+\label{sec:how_to_run_the_simulation_and_analyze_data}
+
+\section{Muonic X-rays in Geant4}
+\label{sec:muonic_x_rays_in_geant4}
+The muonic energy levels and transition probabilities were calculated using
+a simple model described by Mukhopadhyay~\cite{Mukhopadhyay.1977}.
+
+\begin{itemize}
+  \item Energies of K-shell electrons were precisely corrected based on
+    that of hydrogen atom, taking finite size of the nucleus into account
+    % \lstinputlisting[
+      % language=c++, firstline=64, lastline=93,firstnumber=64,
+      % breaklines=true, breakatwhitespace=true,
+      % frame=single]{listings/G4EmCaptureCascade.cc}
+
+      \inputminted[
+        bgcolor=greybg,
+        breaklines=true,
+        stepnumber=5,
+        linenos=true,
+        firstline=64,
+        fontsize=\footnotesize,
+        lastline=93]{c++}{listings/G4EmCaptureCascade.cc}
+  \item Energies of 
+\end{itemize}
+
+\bibliographystyle{h-physrev}
+\bibliography{stm_study}
+\end{document}