Kamil Sedlak 2009-06-05

Added possibility to use mu- (or other particles) as the primary particles
  in the G4ParticleGun  (this option can be switched on by 
  /gun/primaryparticle ).

doc/musrSim.tex

@@ -45,13 +45,13 @@ Geant4}. The root output variables are also described.
 %========================================================================================================
 \section{Scope of the musrSim program}
 The program ``musrSim'' is a relatively general program that can be used to simulate 
-the response of a $\mu$SR~\cite{Blundel} instruments (detectors) to muons and their decay particles 
+the response of a $\mu$SR~\cite{Blundel:1999} instruments (detectors) to muons and their decay particles 
 (electrons, positrons and gammas).  Even though musrSim is tailored to the needs of 
 the $\mu$SR technique~\cite{shirokaGeant}, it has been used also in the detector development
 studies without any muons involved, e.g.\ to test the response of an APD-based 
 scintillator counters to the irradiation of Sr radioactive source~\cite{AlexeyTestAPD}.
 
-The program is based on the Geant4~\cite{Geant} and Root~\cite{Root} libraries.
+The program is based on the Geant4~\cite{geant} and Root~\cite{root} libraries.
 Geant4 is Monte Carlo toolkit used (not only) in particle physics to simulate
 the passage of particles through the detectors.
 Root is an analysis tool that allows to manipulate and analyse the simulated data, 
@@ -65,7 +65,7 @@ In our view, the main advantages of musrSim are:
 	\item Limited (ideally no) need to modify and/or recompile the source code.
 	\item Implementation of the  $\mu$SR-specific classes (muon spin rotation
 	      in magnetic fields, muonium formation and decay, ...).
-	\item Possibility to read in the output files of the TURTLE~\cite{turle}
+	\item Possibility to read in the output files of the TURTLE~\cite{turtle}
 	      program for the beam-line simulation.
 	\item Simple way how to define (overlapping) electromagnetic fields.
 	\item Output in the Root tree.
@@ -97,9 +97,14 @@ To install and run musrSim, one has to install Geant4 and Root first.
 The present version of musrSim has been tested with Geant version 4.9.1, patch no.~3
 and Root version 5.20.00.
 
-The musrSim installation package can be downloaded from the web page http://lmu.web.psi.ch/simulation/index.html.
+Once Geant4 has been successfully installed and some of the default Geant4 exmples
-Once compiled, the program can be executed by \emph{``musrSim NNN.mac''}, where NNN.mac is a ``macro file''
+has been run, the musrSim installation package can be downloaded from the web page 
-containing the information about the instrument setup.  The string NNN represents the integer run number.
+http://lmu.web.psi.ch/simulation/index.html.
+Usually the ``env.sh'' script has to be run to set-up the environment variables
+appropriatelly before the musrSim or any other Geant4 application can be compiled
+or run.
+The simulation can be executed by \emph{``musrSim RUNNUMBER.mac''}, where RUNNUMBER.mac is a ``macro file''
+containing the information about the instrument setup.  The string ``RUNNUMBER'' represents the integer run number.
 
 In order to simulate a new instrument, the user has to define the following blocks of information
 in the macro file:
@@ -120,8 +125,28 @@ in the macro file:
 \end{itemize}
 
 By default, the output of the simulation is written out in the subdirectory ``data'' with
-the name ``musr\_NNN.root''.
+the name ``musr\_RUNNUMBER.root''.
 
+\section{Conventions}
+The default units of the musrSim in both the macro file (RUNNUMBER.mac) and in the Root tree 
+are summarised in table~\ref{tab:units}.
+\begin{table}[htb]\centering
+\renewcommand{\arraystretch}{1.05}
+  \begin{tabular}{lp{5mm}l}
+    \hline
+    \lower 1mm \hbox{\textbf{Quantity}}  &&  \lower 1mm \hbox{\textbf{Default unit}} \\[5pt]
+    \hline
+     Length      &&  mm  \\
+     Time        &&  ns  \\
+     Energy      &&  MeV \\
+     Momentum    &&  MeV/c\\
+     Magnetic field   && tesla \\
+     Electric field   && keV/mm \\
+     \hline
+  \end{tabular}
+  \caption{The default units in musrSim.}
+  \label{tab:units}
+\end{table}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Detector construction}
@@ -322,7 +347,7 @@ Three special volumes ``Target, M0, M1 and M2''.
 		factor.  The units are tesla for the magnetic field and  kV/mm
 		for the electric field.
 	   \item \emph{[fieldValueFinal]} and \emph{[fieldNrOfSteps]} (floats) 
-		-- an experimental feature, which allows the user to ramp up (down) the field
+		-- these optional parameters allow the user to ramp up (down) the field
 		during a single run.  The \emph{fieldValue} serves as the initial field value,
 		the \emph{[fieldValueFinal]} is the final value and \emph{[fieldNrOfSteps]}
 		specifies number of steps, in which the rump up/down will happen.
@@ -381,6 +406,11 @@ Three special volumes ``Target, M0, M1 and M2''.
 \section{Initial (muon) beam parameters}
 
 \begin{description}
+\item{\bf /gun/primaryparticle \emph{primaryParticleName}}\\
+        (default: /gun/primaryparticle mu+)\\
+	Set the primary particle type, if it is not positive muon.  For example, the negative muons
+	are specified by ``/gun/primaryparticle mu-''.
+
 \item{\bf /gun/vertex \emph{x0} \emph{y0} \emph{z0} \emph{unit}}\\
 	(default: /gun/vertex 0 0 -100 mm) \\
 	Set mean values of the $x$, $y$ and $z$ coordinates of the generated particles (muons).
@@ -844,6 +874,23 @@ in~\cite{Aktas:2004px}.
 
 \begin{thebibliography}{0}
 
+\bibitem{Blundel:1999}  S.J.~Blundel,  Contemp. Phys. 40 (1999) 175.
+
+\bibitem{AlexeyTestAPD} A.~Stoykov {\it et al.}, ``First experience with G-APDs in $\mu$SR instrumentation'',
+	NDIP08, to be published in Nucl. Instrum. Meth. A.
+
+\bibitem{geant} S.~Agostinelli, et al., Nucl. Instr. and Meth. A 506 (2003) 250.\\ %-303. \\
+                J.~Allison, et al., IEEE Trans. Nucl.\ Sci.\  53 (2006) 270. %-278.
+\bibitem{root} R.~Brun, F.~Rademakers ``ROOT - An Object Oriented Data Analysis Framework'',
+%Proceedings AIHENP'96 Workshop, Lausanne, Sep. 1996, 
+Nucl. Inst. and Meth. in Phys. Res. A 389 (1997) 81.%-86. 
+See also http://root.cern.ch/.
+
+\bibitem{turtle}
+K.L.~Brown, Ch.~Iselin, D.C.~Carey, {\it``Decay Turtle''}, CERN 74-2 (1974). \\
+U.~Rohrer, {\it ``Compendium of Decay Turtle Enhancements''},
+http://pc532.psi.ch/turtcomp.htm
+
 \bibitem{shirokaGeant} 
 T.~Shiroka {\it et al.} ``GEANT4 as a simulation framework in muSR'',
 Physica {\bf B~404}, (2009) 966-969

include/musrPrimaryGeneratorAction.hh

@@ -74,7 +74,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     void SetOrReadTheRandomNumberSeeds(G4int eventID);
     void SetTurtleMomentumBite (G4ThreeVector smearingParam)
                                {turtleMomentumBite=true; turtleMomentumP0=smearingParam[0]*MeV; turtleSmearingFactor=smearingParam[1]*0.01;}
-
+    void SetPrimaryParticule(G4String particleName);
     static G4String GetPrimaryName();                
 
   private:
@@ -85,7 +85,8 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     musrPrimaryGeneratorMessenger* gunMessenger;  // messenger of this class
     G4String                      rndmFlag;	  // flag for a random impact point       
 
-    // cks Implement also alpha and proton particles for the simulation of Juan Pablo Urrego
+  // cks delete    G4ParticleDefinition* muonMinusParticle;
+    // cks Alpha and proton particles implemented for the simulation of Juan Pablo Urrego
     G4ParticleDefinition* alphaParticle;
     G4ParticleDefinition* protonParticle;
     // csk

include/musrPrimaryGeneratorMessenger.hh

@@ -46,7 +46,8 @@ class musrPrimaryGeneratorMessenger: public G4UImessenger
     void SetNewValue(G4UIcommand*, G4String);
     
   private:
-    musrPrimaryGeneratorAction* musrAction; 
+    musrPrimaryGeneratorAction* musrAction;
+    G4UIcmdWithAString*        setPrimaryParticleCmd; 
     G4UIcmdWith3VectorAndUnit* setvertexCmd;
     G4UIcmdWith3VectorAndUnit* setvertexSigmaCmd;
     G4UIcmdWith3VectorAndUnit* setvertexBoundaryCmd;

musrSim.cc

@@ -139,8 +139,8 @@ int main(int argc,char** argv) {
   delete myRootOutput;
   delete myErrorMessage;
   delete myParameters;
+  // cks  runManager->SetVerboseLevel(2);   // This line can help debug crashes during the runManager delete
   delete runManager;
-  
   //  F04GlobalField* myGlobalField = F04GlobalField::getObject();
   //  if (myGlobalField!=NULL)  {delete myGlobalField;}
 

src/musrPhysicsList.cc

@@ -113,9 +113,9 @@ void musrPhysicsList::ConstructLeptons()
   //    MuonPlusDecayTable -> Insert(new G4MuonRadiativeDecayChannelWithSpin("mu+",0.014));
   G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
   //
-  //  G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
+  G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
-  //  MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00));
+  MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00));
-  //  G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
+  G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
   //csk
   //
   //TS: Using the muonium decay with and without spin

src/musrPrimaryGeneratorAction.cc

@@ -77,6 +77,7 @@ musrPrimaryGeneratorAction::musrPrimaryGeneratorAction(
   // default particle kinematic
   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
   G4ParticleDefinition* muonParticle= particleTable->FindParticle("mu+");
+  // cks delete muonMinusParticle= particleTable->FindParticle("mu-");
   // cks Implement also alpha and proton particles for the simulation of Juan Pablo Urrego
   alphaParticle= particleTable->FindParticle("alpha");
   protonParticle= particleTable->FindParticle("proton");
@@ -448,3 +449,20 @@ void musrPrimaryGeneratorAction::SetKEnergy(G4double val)  {
   //  G4cout<<"musrPrimaryGeneratorAction::SetKEnergy: Muon kinetic energy of "
   //	<<val<<" MeV requested ==> initial muon momentum set to "<<p0<<" MeV/c"<<G4endl;
 }
+
+//===============================================================================
+void musrPrimaryGeneratorAction::SetPrimaryParticule(G4String particleName) {
+  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
+  G4ParticleDefinition* particule= particleTable->FindParticle(particleName);
+  if (particule!=NULL) {
+    G4cout<<"   musrPrimaryGeneratorAction::SetPrimaryParticleMuMinus():  USING "
+	  <<particleName<<" AS PRIMARY PARTICLE!"<<G4endl; 
+    particleGun->SetParticleDefinition(particule);
+  }
+  else {
+    G4cout<<"\n   musrPrimaryGeneratorAction::SetPrimaryParticleMuMinus():  Particle "
+	  <<particleName<<" required as primary particle, but not found !!!";
+    G4cout<<"S T O P     F O R C E D" << G4endl;
+    exit(1);
+  }
+}

src/musrPrimaryGeneratorMessenger.cc

@@ -35,7 +35,10 @@
 musrPrimaryGeneratorMessenger::musrPrimaryGeneratorMessenger(musrPrimaryGeneratorAction* musrGun)
 :musrAction(musrGun)
 { 
- 
+  setPrimaryParticleCmd = new G4UIcmdWithAString("/gun/primaryparticle",this);
+  setPrimaryParticleCmd -> SetGuidance("Set the primary particle (default is mu+).");
+  setPrimaryParticleCmd ->AvailableForStates(G4State_PreInit,G4State_Idle);
+
   setvertexCmd = new G4UIcmdWith3VectorAndUnit("/gun/vertex",this);
   setvertexCmd->SetGuidance(" Set x0, y0, z0 of the generated muons (with unit)");
   setvertexCmd->SetParameterName("mes_x0","mes_y0","mes_z0",true,true);
@@ -137,6 +140,7 @@ musrPrimaryGeneratorMessenger::musrPrimaryGeneratorMessenger(musrPrimaryGenerato
 
 musrPrimaryGeneratorMessenger::~musrPrimaryGeneratorMessenger()
 {
+  delete setPrimaryParticleCmd;
   delete setvertexCmd;
   delete setvertexSigmaCmd;
   delete setvertexBoundaryCmd;
@@ -160,6 +164,8 @@ musrPrimaryGeneratorMessenger::~musrPrimaryGeneratorMessenger()
 
 void musrPrimaryGeneratorMessenger::SetNewValue(G4UIcommand * command,G4String newValue)
 { 
+  if( command == setPrimaryParticleCmd)
+    { musrAction->SetPrimaryParticule(newValue); }
   if( command == setvertexCmd)
    { musrAction->Setvertex(setvertexCmd->GetNew3VectorValue(newValue));}
   if( command == setvertexSigmaCmd)

src/musrScintSD.cc

@@ -245,7 +245,7 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
 	    det_edep_pos[j]           += aHit->GetEdep();
 	  }  else if (aParticleName=="gamma") {
 	    det_edep_gam[j]           += aHit->GetEdep();
-	  }  else if (aParticleName=="mu+") {
+	  }  else if ((aParticleName=="mu+")||(aParticleName=="mu-")) {
 	    det_edep_mup[j]           += aHit->GetEdep();
 	  }  else {
 	    char message[200];
@@ -290,7 +290,7 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
 	    det_edep_pos[nSignals]           += aHit->GetEdep();
 	  }  else if (aParticleName=="gamma") {
 	    det_edep_gam[nSignals]           += aHit->GetEdep();
-	  }  else if (aParticleName=="mu+") {
+	  }  else if ((aParticleName=="mu+")||(aParticleName=="mu-")) {
 	    det_edep_mup[nSignals]           += aHit->GetEdep();
 	  }  else {
 	    char message[200];