diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index 5962a05..d9b5be4 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index b1f5da9..a0d73f8 100644 --- a/doc/musrSim.tex +++ b/doc/musrSim.tex @@ -1,4 +1,5 @@ \documentclass[twoside]{dis04} +\usepackage{epsfig} %\def\runauthor{Kamil Sedl\'{a}k} \def\runauthor{PSI} % for the H1 collaboration} @@ -320,13 +321,19 @@ Three special volumes ``Target, M0, M1 and M2''. The following lines specify the \emph{x, y, z, Field\_x, Field\_y, Field\_z} values (non-compact format) or \emph{Field\_x, Field\_y, Field\_z} values (compact format).\\ - {\bf 3DBOpera} -- 3D magnetic field in the form of OPERA output. + {\bf 3DBOpera, 3DEOpera} -- 3D magnetic field in the form of OPERA output. It is expected that the \emph{length unit} is 1\,m, and the \emph{field normalisation factor} is 1. (Note that this default normalisation is different from 2DBOpera and 2DBOperaXY options). However, a different \emph{field normalisation factor} can be specified in the field map file using the keyword ``fieldNormalisation \emph{number}'' before the line started with 0.\\ + It is expected that the we first loop over the $z$ coordinate of the field map, + then (when $z$ changed from minimum to maximum) it is looped over $y$ coordinate, + and the highest-lever loop goes over $x$ coordinate. Hoever, if the order + of looping is reversed in the field map, it can be specified using the + keyword ``variableIncreasingOrder xyz'' placed in the field map before + the line started with 0.\\ The \emph{length unit} can be changed to 1\,cm by specifying ``[CENTIMETRE]'' after the ``0'' character in the field map file.\\ It is expected that the field map is defined in the full volume of the field. @@ -888,6 +895,7 @@ The list of variables that can be stored in the Root tree: the given ``save'' volume. Save volumes can therefore be made of vacuum. \item{\bf save\_detID[save\_n]} (array of Int\_t) -- ID number of the save volume. \item{\bf save\_particleID[save\_n]} (array of Int\_t) -- particle ID of the particle that entered the save volume. +\item{\bf save\_time[save\_n]} (array of Double\_t) -- time when the particle enetered in the volume (in $\mu$s). \item{\bf save\_x[save\_n], save\_y[save\_n], save\_z[save\_n]} (array of Double\_t) -- position of the particle where it entered the save volume (``GetPreStepPoint()'') (in mm). \item{\bf save\_px[save\_n], save\_py[save\_n], save\_pz[save\_n]} (array of Double\_t) -- momentum of the particle when it @@ -941,6 +949,24 @@ These are typically not interesting for the analysis of the results of the simul ``musrDetectorConstruction::SetColourOfLogicalVolume''. \end{description} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\section{Example 1 -- Electrons passing through two scintillator tiles (101.mac)} +One of the easiest example to illustrate the basic features of the musrSim (and/or Geant4) is to +shoot electrons into a scintillator block, and to observe the energy deposited inside it. +Figure~\ref{vis101} +% +\begin{figure}[tb]\centering +\epsfig{file=pict/vis_101_c.eps,width=8cm,%\linewidth,% +%bbllx=83pt,bblly=330pt,bburx=538pt,bbury=513pt, +clip=} +\caption{A simple simulation of an electron passing through two +scintillator tiles.} +\label{vis101} +\end{figure} +% +illustrates a simple geometry made of an electron source and two blocks +of scintilator tiles with the dimensions of $3 \times 3 \times 2\,$mm, +which is defined in the macro file ``101.mac''. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/doc/pict/vis_101_a.eps b/doc/pict/vis_101_a.eps new file mode 100644 index 0000000..202451e --- /dev/null +++ b/doc/pict/vis_101_a.eps @@ -0,0 +1,138 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%BoundingBox: 130 365 530 480 +% FILE: g4_09.eps + +% DEFAULT FONT: +/Times-Roman findfont 24 scalefont setfont + +% DEFAULT COLOR: +0 0 0 setrgbcolor + +% DEFAULT LINE: +0.283464567 setlinewidth +[1.41732283] 0 setdash + +% DEFINE PostScript Functions: +/PATH +{/POINT 0 def + /X 0 def + newpath + {X 1 eq{POINT 0 eq{moveto + /POINT 1 def + }{lineto }ifelse + /X 0 def + }{/X X 1 add def}ifelse + }forall +} def +/FILLPOLYGON{ PATH closepath fill }def +/DRAWFILLPOLYGON{ PATH closepath gsave stroke grestore fill }def +/DRAWPOLYGON{ PATH closepath stroke }def +/DRAWPATH { PATH stroke } def +/BOXPATH +{/X 0 def + {X 0 eq{/LLX exch def}if + X 1 eq{/LLY exch def}if + X 2 eq{/URX exch def}if + X 3 eq{/URY exch def}if + /X X 1 add def + }forall + LLX LLY moveto + URX LLY lineto + URX URY lineto + LLX URY lineto + closepath +}def + +% CLIPPING: +[56.6929134 56.6929134 538.582677 785.19685] BOXPATH clip + +% PARAMETERS: +% Generated by Fukui Renderer DAWN version 3.88a (Dev. indep. Mode) +% O.x = 0 +% O.y = 0 +% O.z = 9.5 +% camera_d = 6845.21 +% v_angle = 80 +% h_angle = 10 +% focal_d = 337.582961 +% e_2d = 1e-07 +% e_3d = 0.001 + +0 0.605994623 0.605994623 setrgbcolor +0 setlinewidth +[] 0 setdash +[ 385.804284 382.325119 432.58785 383.759675 432.597965 454.010008 385.810893 452.57891 ] DRAWFILLPOLYGON +0 0.807830738 0.807830738 setrgbcolor +[ 385.810893 452.57891 432.597965 454.010008 419.979122 465.993597 373.211292 464.563714 ] DRAWFILLPOLYGON +0 1 1 setrgbcolor +[ 385.804284 382.325119 385.810893 452.57891 373.211292 464.563714 373.205629 394.338877 ] DRAWFILLPOLYGON +0 0.605994623 0.605994623 setrgbcolor +[ 175.219473 375.867816 222.02441 377.303027 222.018741 447.568921 175.210295 446.137168 ] DRAWFILLPOLYGON +0 0.807830738 0.807830738 setrgbcolor +[ 175.210295 446.137168 222.018741 447.568921 209.486516 459.557983 162.697329 458.127447 ] DRAWFILLPOLYGON +0 1 1 setrgbcolor +[ 175.219473 375.867816 175.210295 446.137168 162.697329 458.127447 162.707441 387.88706 ] DRAWFILLPOLYGON +0 0 0 setrgbcolor +0.850393701 setlinewidth +[ 385.804284 382.325119 432.58785 383.759675] DRAWPATH +[ 432.58785 383.759675 432.597965 454.010008] DRAWPATH +[ 432.597965 454.010008 385.810893 452.57891] DRAWPATH +[ 385.810893 452.57891 385.804284 382.325119] DRAWPATH +[ 385.810893 452.57891 432.597965 454.010008] DRAWPATH +[ 432.597965 454.010008 419.979122 465.993597] DRAWPATH +[ 419.979122 465.993597 373.211292 464.563714] DRAWPATH +[ 373.211292 464.563714 385.810893 452.57891] DRAWPATH +[ 385.804284 382.325119 385.810893 452.57891] DRAWPATH +[ 385.810893 452.57891 373.211292 464.563714] DRAWPATH +[ 373.211292 464.563714 373.205629 394.338877] DRAWPATH +[ 373.205629 394.338877 385.804284 382.325119] DRAWPATH +[ 175.219473 375.867816 222.02441 377.303027] DRAWPATH +[ 222.02441 377.303027 222.018741 447.568921] DRAWPATH +[ 222.018741 447.568921 175.210295 446.137168] DRAWPATH +[ 175.210295 446.137168 175.219473 375.867816] DRAWPATH +[ 175.210295 446.137168 222.018741 447.568921] DRAWPATH +[ 222.018741 447.568921 209.486516 459.557983] DRAWPATH +[ 209.486516 459.557983 162.697329 458.127447] DRAWPATH +[ 162.697329 458.127447 175.210295 446.137168] DRAWPATH +[ 175.219473 375.867816 175.210295 446.137168] DRAWPATH +[ 175.210295 446.137168 162.697329 458.127447] DRAWPATH +[ 162.697329 458.127447 162.707441 387.88706] DRAWPATH +[ 162.707441 387.88706 175.219473 375.867816] DRAWPATH +1 0 1 setrgbcolor +[ 168.54679 406.309781 -48.7057921 160.988791] DRAWPATH +[ 168.54679 406.309781 -48.7057921 160.988791] DRAWPATH +[ 168.54679 406.309781 168.54679 406.309781] DRAWPATH +[ 168.54679 406.309781 168.54679 406.309781] DRAWPATH +[ 168.446039 408.225987 168.446039 408.225987] DRAWPATH +[ 168.446039 408.225987 168.446039 408.225987] DRAWPATH +[ 380.22515 425.545444 222.021947 407.830611] DRAWPATH +[ 380.22515 425.545444 222.021947 407.830611] DRAWPATH +[ 380.22515 425.545444 380.22515 425.545444] DRAWPATH +[ 380.22515 425.545444 380.22515 425.545444] DRAWPATH +[ 380.285088 426.103451 380.285088 426.103451] DRAWPATH +[ 380.285088 426.103451 380.285088 426.103451] DRAWPATH +[ 519.819541 427.747833 432.593799 425.077081] DRAWPATH +[ 519.819541 427.747833 432.593799 425.077081] DRAWPATH +1 0 0 setrgbcolor +newpath 426.282383 424.883832 0.5 0 360 arc fill +newpath 384.488285 425.992003 0.5 0 360 arc fill +newpath 215.875845 407.1424 0.5 0 360 arc fill +newpath 192.909649 410.155967 0.5 0 360 arc fill +newpath 170.403911 408.380447 0.5 0 360 arc fill +1 1 0 setrgbcolor +newpath 519.819541 427.747833 0.5 0 360 arc fill +newpath 426.282383 424.883832 0.5 0 360 arc fill +newpath 384.488285 425.992003 0.5 0 360 arc fill +newpath 380.22515 425.545444 0.5 0 360 arc fill +newpath 215.875845 407.1424 0.5 0 360 arc fill +newpath 192.909649 410.155967 0.5 0 360 arc fill +newpath 170.403911 408.380447 0.5 0 360 arc fill +newpath 168.54679 406.309781 0.5 0 360 arc fill +newpath -48.7057921 160.988791 0.5 0 360 arc fill +1 0 0 setrgbcolor +newpath 426.282383 424.883832 0.5 0 360 arc fill +newpath 384.488285 425.992003 0.5 0 360 arc fill +newpath 215.875845 407.1424 0.5 0 360 arc fill +newpath 192.909649 410.155967 0.5 0 360 arc fill +newpath 170.403911 408.380447 0.5 0 360 arc fill +showpage diff --git a/doc/pict/vis_101_b.eps b/doc/pict/vis_101_b.eps new file mode 100644 index 0000000..4dbeb51 --- /dev/null +++ b/doc/pict/vis_101_b.eps @@ -0,0 +1,138 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%BoundingBox: 60 365 505 480 +% FILE: g4_09.eps + +% DEFAULT FONT: +/Times-Roman findfont 24 scalefont setfont + +% DEFAULT COLOR: +0 0 0 setrgbcolor + +% DEFAULT LINE: +0.283464567 setlinewidth +[1.41732283] 0 setdash + +% DEFINE PostScript Functions: +/PATH +{/POINT 0 def + /X 0 def + newpath + {X 1 eq{POINT 0 eq{moveto + /POINT 1 def + }{lineto }ifelse + /X 0 def + }{/X X 1 add def}ifelse + }forall +} def +/FILLPOLYGON{ PATH closepath fill }def +/DRAWFILLPOLYGON{ PATH closepath gsave stroke grestore fill }def +/DRAWPOLYGON{ PATH closepath stroke }def +/DRAWPATH { PATH stroke } def +/BOXPATH +{/X 0 def + {X 0 eq{/LLX exch def}if + X 1 eq{/LLY exch def}if + X 2 eq{/URX exch def}if + X 3 eq{/URY exch def}if + /X X 1 add def + }forall + LLX LLY moveto + URX LLY lineto + URX URY lineto + LLX URY lineto + closepath +}def + +% CLIPPING: +[56.6929134 56.6929134 538.582677 785.19685] BOXPATH clip + +% PARAMETERS: +% Generated by Fukui Renderer DAWN version 3.88a (Dev. indep. Mode) +% O.x = 0 +% O.y = 0 +% O.z = 9.5 +% camera_d = 6845.21 +% v_angle = 80 +% h_angle = 170.5 +% focal_d = 340.001383 +% e_2d = 1e-07 +% e_3d = 0.001 + +0 0.5 0.5 setrgbcolor +0 setlinewidth +[] 0 setdash +[ 208.838521 453.027486 161.714149 454.395364 161.723832 383.538844 208.844847 382.167468 ] DRAWFILLPOLYGON +0 0.807830738 0.807830738 setrgbcolor +[ 221.510042 464.500028 174.405081 465.866742 161.714149 454.395364 208.838521 453.027486 ] DRAWFILLPOLYGON +0 1 1 setrgbcolor +[ 221.515463 393.669347 221.510042 464.500028 208.838521 453.027486 208.844847 382.167468 ] DRAWFILLPOLYGON +0 0.5 0.5 setrgbcolor +[ 420.957389 446.870317 373.811489 448.23882 373.806062 377.366557 420.948601 375.994555 ] DRAWFILLPOLYGON +0 0.807830738 0.807830738 setrgbcolor +[ 433.54151 458.3481 386.41504 459.715438 373.811489 448.23882 420.957389 446.870317 ] DRAWFILLPOLYGON +0 1 1 setrgbcolor +[ 433.53183 387.501687 433.54151 458.3481 420.957389 446.870317 420.948601 375.994555 ] DRAWFILLPOLYGON +0 0 0 setrgbcolor +0.850393701 setlinewidth +[ 208.838521 453.027486 161.714149 454.395364] DRAWPATH +[ 161.714149 454.395364 161.723832 383.538844] DRAWPATH +[ 161.723832 383.538844 208.844847 382.167468] DRAWPATH +[ 208.844847 382.167468 208.838521 453.027486] DRAWPATH +[ 221.510042 464.500028 174.405081 465.866742] DRAWPATH +[ 174.405081 465.866742 161.714149 454.395364] DRAWPATH +[ 161.714149 454.395364 208.838521 453.027486] DRAWPATH +[ 208.838521 453.027486 221.510042 464.500028] DRAWPATH +[ 221.515463 393.669347 221.510042 464.500028] DRAWPATH +[ 221.510042 464.500028 208.838521 453.027486] DRAWPATH +[ 208.838521 453.027486 208.844847 382.167468] DRAWPATH +[ 208.844847 382.167468 221.515463 393.669347] DRAWPATH +[ 420.957389 446.870317 373.811489 448.23882] DRAWPATH +[ 373.811489 448.23882 373.806062 377.366557] DRAWPATH +[ 373.806062 377.366557 420.948601 375.994555] DRAWPATH +[ 420.948601 375.994555 420.957389 446.870317] DRAWPATH +[ 433.54151 458.3481 386.41504 459.715438] DRAWPATH +[ 386.41504 459.715438 373.811489 448.23882] DRAWPATH +[ 373.811489 448.23882 420.957389 446.870317] DRAWPATH +[ 420.957389 446.870317 433.54151 458.3481] DRAWPATH +[ 433.53183 387.501687 433.54151 458.3481] DRAWPATH +[ 433.54151 458.3481 420.957389 446.870317] DRAWPATH +[ 420.957389 446.870317 420.948601 375.994555] DRAWPATH +[ 420.948601 375.994555 433.53183 387.501687] DRAWPATH +1 0 1 setrgbcolor +[ 426.72948 407.356942 426.72948 407.356942] DRAWPATH +[ 426.830321 405.613781 426.830321 405.613781] DRAWPATH +[ 426.72948 407.356942 426.72948 407.356942] DRAWPATH +[ 426.830321 405.613781 426.830321 405.613781] DRAWPATH +[ 426.830321 405.613781 676.816453 188.228361] DRAWPATH +[ 426.830321 405.613781 676.816453 188.228361] DRAWPATH +[ 215.900622 426.797536 373.808419 408.141437] DRAWPATH +[ 215.900622 426.797536 373.808419 408.141437] DRAWPATH +[ 215.900622 426.797536 215.900622 426.797536] DRAWPATH +[ 215.900622 426.797536 215.900622 426.797536] DRAWPATH +[ 215.96079 427.472432 215.96079 427.472432] DRAWPATH +[ 215.96079 427.472432 215.96079 427.472432] DRAWPATH +[ 73.8539467 427.447758 161.718181 424.894535] DRAWPATH +[ 73.8539467 427.447758 161.718181 424.894535] DRAWPATH +1 0 0 setrgbcolor +newpath 168.06558 424.710087 0.5 0 360 arc fill +newpath 211.584374 427.220024 0.5 0 360 arc fill +newpath 380.231921 407.382529 0.5 0 360 arc fill +newpath 402.593632 409.737525 0.5 0 360 arc fill +newpath 424.797809 407.547468 0.5 0 360 arc fill +1 1 0 setrgbcolor +newpath 73.8539467 427.447758 0.5 0 360 arc fill +newpath 168.06558 424.710087 0.5 0 360 arc fill +newpath 211.584374 427.220024 0.5 0 360 arc fill +newpath 215.900622 426.797536 0.5 0 360 arc fill +newpath 380.231921 407.382529 0.5 0 360 arc fill +newpath 402.593632 409.737525 0.5 0 360 arc fill +newpath 424.797809 407.547468 0.5 0 360 arc fill +newpath 426.830321 405.613781 0.5 0 360 arc fill +newpath 676.816453 188.228361 0.5 0 360 arc fill +1 0 0 setrgbcolor +newpath 168.06558 424.710087 0.5 0 360 arc fill +newpath 211.584374 427.220024 0.5 0 360 arc fill +newpath 380.231921 407.382529 0.5 0 360 arc fill +newpath 402.593632 409.737525 0.5 0 360 arc fill +newpath 424.797809 407.547468 0.5 0 360 arc fill +showpage diff --git a/include/musrRootOutput.hh b/include/musrRootOutput.hh index 3e8e898..c4e8d12 100644 --- a/include/musrRootOutput.hh +++ b/include/musrRootOutput.hh @@ -78,7 +78,7 @@ class musrRootOutput { G4double ekVertex, G4double xVertex, G4double yVertex, G4double zVertex, G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex, G4int particleID) ; - void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z, + void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z, G4double time, G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) ; void SetInitialMuonParameters(G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz, @@ -304,6 +304,7 @@ class musrRootOutput { G4int save_detID[save_nMax]; G4int save_particleID[save_nMax]; G4double save_ke[save_nMax]; + G4double save_time[save_nMax]; G4double save_x[save_nMax]; G4double save_y[save_nMax]; G4double save_z[save_nMax]; diff --git a/include/musrTabulatedElementField.hh b/include/musrTabulatedElementField.hh index 5c7a762..3950358 100644 --- a/include/musrTabulatedElementField.hh +++ b/include/musrTabulatedElementField.hh @@ -85,7 +85,8 @@ private: double maximumWidth, maximumHeight, maximumLength; int symmetryType; // this variable defines whether (and how) should be the field extended // in the case it is defined just in one octant of the Cartesian coordinates. - + char variableIncreasingOrder[100]; + int jx, jy, jz; void Invert(const char* indexToInvert); }; diff --git a/run/101.mac b/run/101.mac new file mode 100644 index 0000000..c5104ca --- /dev/null +++ b/run/101.mac @@ -0,0 +1,203 @@ +#----------------------------------------------------------------------- +# Macro file for the simulation of electron/positrons from the Sr decay +# passing through the scintiallator detectors (blocks). +# Unless specified otherwises, the default units are mm, ns, MeV, MeV/c. +# Lines starting with star "#" are comments. +################################################################################### +############################# G E O M E T R Y ################################### +# +# WORLD +/musr/command construct box World 10 10 100 G4_AIR 0 0 0 no_logical_volume norot dead -1 +# +# Sr SOURCE +#/musr/command construct sphere source 0 0.02 0 360 0 180 G4_Sr 0 0 0 log_World norot dead 301 +# +# SCINTILLATOR BLOCKS +/musr/command construct box scintFarAwayC1 1.5 1.5 1 G4_PLASTIC_SC_VINYLTOLUENE 0 0 5 log_World norot musr/ScintSD 10 +/musr/command construct box scintFarAwayC2 1.5 1.5 1 G4_PLASTIC_SC_VINYLTOLUENE 0 0 14 log_World norot musr/ScintSD 11 +# +#============================================================ +/musr/command visattributes log_World invisible +/musr/command visattributes log_source red +/musr/command visattributes G4_PLASTIC_SC_VINYLTOLUENE lightblue +################################################################################### +######################### P H Y S I C S P R O C E S S E S ################## +################################################################################### +# --- Low Energy (default) --- +/musr/command process addDiscreteProcess gamma G4LowEnergyPhotoElectric +/musr/command process addDiscreteProcess gamma G4LowEnergyCompton +/musr/command process addDiscreteProcess gamma G4LowEnergyGammaConversion +/musr/command process addDiscreteProcess gamma G4LowEnergyRayleigh +/musr/command process addProcess e- G4MultipleScattering -1 1 1 +#/musr/command process addDiscreteProcess e- G4CoulombScattering +/musr/command process addProcess e- G4LowEnergyIonisation -1 2 2 +/musr/command process addProcess e- G4LowEnergyBremsstrahlung -1 -1 3 +/musr/command process addProcess e+ G4MultipleScattering -1 1 1 +#/musr/command process addDiscreteProcess e+ G4CoulombScattering +/musr/command process addProcess e+ G4eIonisation -1 2 2 +/musr/command process addProcess e+ G4eBremsstrahlung -1 3 3 +/musr/command process addProcess e+ G4eplusAnnihilation 0 -1 4 +# +# --- High Energy --- +#/musr/command process addDiscreteProcess gamma G4PhotoElectricEffect +#/musr/command process addDiscreteProcess gamma G4ComptonScattering +#/musr/command process addDiscreteProcess gamma G4GammaConversion +#/musr/command process addProcess e- G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e- G4CoulombScattering +#/musr/command process addProcess e- G4eIonisation -1 2 2 +#/musr/command process addProcess e- G4eBremsstrahlung -1 3 3 +#/musr/command process addProcess e+ G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e+ G4CoulombScattering +#/musr/command process addProcess e+ G4eIonisation -1 2 2 +#/musr/command process addProcess e+ G4eBremsstrahlung -1 3 3 +#/musr/command process addProcess e+ G4eplusAnnihilation 0 -1 4 +# +# --- Penelope --- +#/musr/command process addDiscreteProcess gamma G4PenelopePhotoElectric +#/musr/command process addDiscreteProcess gamma G4PenelopeCompton +#/musr/command process addDiscreteProcess gamma G4PenelopeGammaConversion +#/musr/command process addDiscreteProcess gamma G4PenelopeRayleigh +#/musr/command process addProcess e- G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e- G4CoulombScattering +#/musr/command process addProcess e- G4PenelopeIonisation -1 2 2 +#/musr/command process addProcess e- G4PenelopeBremsstrahlung -1 -1 3 +#/musr/command process addProcess e+ G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e+ G4CoulombScattering +#/musr/command process addProcess e+ G4PenelopeIonisation, -1 2 2 +#/musr/command process addProcess e+ G4PenelopeBremsstrahlung, -1 -1 3 +#/musr/command process addProcess e+ G4PenelopeAnnihilation, 0 -1 4 +# +# --- Muons --- +#/musr/command process addProcess mu+ G4MultipleScattering -1 1 1 +##/musr/command process addProcess mu+ MultipleAndCoulombScattering -1 1 1 goulombRegion +##/musr/command process addDiscreteProcess mu+ G4CoulombScattering +#/musr/command process addProcess mu+ G4MuIonisation -1 2 2 +#/musr/command process addProcess mu+ G4MuBremsstrahlung -1 3 3 +#/musr/command process addProcess mu+ G4MuPairProduction -1 4 4 +#/musr/command process addProcess mu- G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess mu- G4CoulombScattering +#/musr/command process addProcess mu- G4MuIonisation -1 2 2 +#/musr/command process addProcess mu- G4MuBremsstrahlung -1 3 3 +#/musr/command process addProcess mu- G4MuPairProduction -1 4 4 +################################################################################### +################## S O M E O T H E R P A R A M E T E R S ################## +################################################################################### +# Store all events into the ROOT tree or just the interesting ones ? (true is default) +#/musr/command storeOnlyEventsWithHits false +/musr/command storeOnlyEventsWithHitInDetID 11 +# Set the minimum time separation between two subsequent signals in the same detector (in ns) +/musr/command signalSeparationTime 50 +# +/musr/run/howOftenToPrintEvent 10000 +/musr/run/randomOption 2 +################################################################################### +######################### R O O T O U T P U T ############################## +################################################################################### +#/musr/command rootOutput runID off +#/musr/command rootOutput eventID off +/musr/command rootOutput weight off +/musr/command rootOutput BFieldAtDecay off +/musr/command rootOutput muIniPosX off +/musr/command rootOutput muIniPosY off +/musr/command rootOutput muIniPosZ off +/musr/command rootOutput muIniMomX off +/musr/command rootOutput muIniMomY off +/musr/command rootOutput muIniMomZ off +/musr/command rootOutput muIniPolX off +/musr/command rootOutput muIniPolY off +/musr/command rootOutput muIniPolZ off +/musr/command rootOutput muIniTime off +/musr/command rootOutput muDecayDetID off +/musr/command rootOutput muDecayPosX off +/musr/command rootOutput muDecayPosY off +/musr/command rootOutput muDecayPosZ off +/musr/command rootOutput muDecayTime off +/musr/command rootOutput muDecayPolX off +/musr/command rootOutput muDecayPolY off +/musr/command rootOutput muDecayPolZ off +/musr/command rootOutput muTargetTime off +/musr/command rootOutput muTargetPolX off +/musr/command rootOutput muTargetPolY off +/musr/command rootOutput muTargetPolZ off +/musr/command rootOutput muM0Time off +/musr/command rootOutput muM0PolX off +/musr/command rootOutput muM0PolY off +/musr/command rootOutput muM0PolZ off +/musr/command rootOutput muM1Time off +/musr/command rootOutput muM1PolX off +/musr/command rootOutput muM1PolY off +/musr/command rootOutput muM1PolZ off +/musr/command rootOutput muM2Time off +/musr/command rootOutput muM2PolX off +/musr/command rootOutput muM2PolY off +/musr/command rootOutput muM2PolZ off +#/musr/command rootOutput posIniMomX off +#/musr/command rootOutput posIniMomY off +#/musr/command rootOutput posIniMomZ off +/musr/command rootOutput fieldNomVal off +#/musr/command rootOutput det_ID off +#/musr/command rootOutput det_edep off +/musr/command rootOutput det_edep_el off +/musr/command rootOutput det_edep_pos off +/musr/command rootOutput det_edep_gam off +/musr/command rootOutput det_edep_mup off +#/musr/command rootOutput det_nsteps off +#/musr/command rootOutput det_length off +#/musr/command rootOutput det_start off +#/musr/command rootOutput det_end off +#/musr/command rootOutput det_x off +#/musr/command rootOutput det_y off +#/musr/command rootOutput det_z off +#/musr/command rootOutput det_kine off +/musr/command rootOutput det_VrtxKine off +/musr/command rootOutput det_VrtxX off +/musr/command rootOutput det_VrtxY off +/musr/command rootOutput det_VrtxZ off +/musr/command rootOutput det_VrtxVolID off +/musr/command rootOutput det_VrtxProcID off +/musr/command rootOutput det_VrtxTrackID off +/musr/command rootOutput det_VrtxParticleID off +/musr/command rootOutput det_VvvKine off +/musr/command rootOutput det_VvvX off +/musr/command rootOutput det_VvvY off +/musr/command rootOutput det_VvvZ off +/musr/command rootOutput det_VvvVolID off +/musr/command rootOutput det_VvvProcID off +/musr/command rootOutput det_VvvTrackID off +/musr/command rootOutput det_VvvParticleID off +### Root variables that are not written out by default, but can be switched on: +#/musr/command rootOutput fieldIntegralBx on +#/musr/command rootOutput fieldIntegralBy on +#/musr/command rootOutput fieldIntegralBz on +#/musr/command rootOutput fieldIntegralBz1 on +#/musr/command rootOutput fieldIntegralBz2 on +#/musr/command rootOutput fieldIntegralBz3 on +################################################################################### +######################### V I S U A L I S A T I O N ############################## +################################################################################### +#/vis/disable +#/control/execute visFromToni.mac +/control/execute visDawn102.mac +#/control/execute visVRML.mac +################################################################################### +######################### P A R T I C L E G U N ################################# +################################################################################### +/gun/primaryparticle e- +/gun/vertex 0 0 0 mm +/gun/momentum 2.15 MeV +# sigma = 3% ==> sigma 27*0.03 = 0.81 +#/gun/momentumsmearing 0.3 MeV +#/gun/tiltsigma 5.15 5.15 0 deg +#/gun/pitch 10.0573 deg +# +#/gps/particle ion +#/gps/ion 39 86 +#/gps/ion 27 57 0 0 +#/gps/ion 38 90 0 0 +# /gps/position seems to be in cm !!!! +#/gps/position 0 0 0 +#/gps/energy 0 keV +#/gps/ang/maxtheta 2 deg +#/gps/ang/maxphi 2 deg +######################## B E A M O N ####################################### +/run/beamOn 10 diff --git a/run/102.mac b/run/102.mac new file mode 100644 index 0000000..e87030d --- /dev/null +++ b/run/102.mac @@ -0,0 +1,195 @@ +#----------------------------------------------------------------------- +# Macro file for the simulation of electron/positrons from the Sr decay +# passing through the scintiallator detectors (blocks). +# Unless specified otherwises, the default units are mm, ns, MeV, MeV/c. +# Lines starting with star "#" are comments. +################################################################################### +############################# G E O M E T R Y ################################### +# +# WORLD +/musr/command construct box World 10 10 100 G4_AIR 0 0 0 no_logical_volume norot dead -1 +# +# Sr SOURCE +/musr/command construct sphere source 0 0.02 0 360 0 180 G4_Sr 0 0 0 log_World norot dead 301 +# +# SCINTILLATOR BLOCKS +/musr/command construct box scintFarAwayC1 1.5 1.5 1 G4_PLASTIC_SC_VINYLTOLUENE 0 0 5 log_World norot musr/ScintSD 10 +/musr/command construct box scintFarAwayC2 1.5 1.5 1 G4_PLASTIC_SC_VINYLTOLUENE 0 0 14 log_World norot musr/ScintSD 11 +# +#============================================================ +/musr/command visattributes log_World invisible +/musr/command visattributes log_source red +/musr/command visattributes G4_PLASTIC_SC_VINYLTOLUENE lightblue +################################################################################### +######################### P H Y S I C S P R O C E S S E S ################## +################################################################################### +# --- Low Energy (default) --- +/musr/command process addDiscreteProcess gamma G4LowEnergyPhotoElectric +/musr/command process addDiscreteProcess gamma G4LowEnergyCompton +/musr/command process addDiscreteProcess gamma G4LowEnergyGammaConversion +/musr/command process addDiscreteProcess gamma G4LowEnergyRayleigh +/musr/command process addProcess e- G4MultipleScattering -1 1 1 +#/musr/command process addDiscreteProcess e- G4CoulombScattering +/musr/command process addProcess e- G4LowEnergyIonisation -1 2 2 +/musr/command process addProcess e- G4LowEnergyBremsstrahlung -1 -1 3 +/musr/command process addProcess e+ G4MultipleScattering -1 1 1 +#/musr/command process addDiscreteProcess e+ G4CoulombScattering +/musr/command process addProcess e+ G4eIonisation -1 2 2 +/musr/command process addProcess e+ G4eBremsstrahlung -1 3 3 +/musr/command process addProcess e+ G4eplusAnnihilation 0 -1 4 +# +# --- High Energy --- +#/musr/command process addDiscreteProcess gamma G4PhotoElectricEffect +#/musr/command process addDiscreteProcess gamma G4ComptonScattering +#/musr/command process addDiscreteProcess gamma G4GammaConversion +#/musr/command process addProcess e- G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e- G4CoulombScattering +#/musr/command process addProcess e- G4eIonisation -1 2 2 +#/musr/command process addProcess e- G4eBremsstrahlung -1 3 3 +#/musr/command process addProcess e+ G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e+ G4CoulombScattering +#/musr/command process addProcess e+ G4eIonisation -1 2 2 +#/musr/command process addProcess e+ G4eBremsstrahlung -1 3 3 +#/musr/command process addProcess e+ G4eplusAnnihilation 0 -1 4 +# +# --- Penelope --- +#/musr/command process addDiscreteProcess gamma G4PenelopePhotoElectric +#/musr/command process addDiscreteProcess gamma G4PenelopeCompton +#/musr/command process addDiscreteProcess gamma G4PenelopeGammaConversion +#/musr/command process addDiscreteProcess gamma G4PenelopeRayleigh +#/musr/command process addProcess e- G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e- G4CoulombScattering +#/musr/command process addProcess e- G4PenelopeIonisation -1 2 2 +#/musr/command process addProcess e- G4PenelopeBremsstrahlung -1 -1 3 +#/musr/command process addProcess e+ G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess e+ G4CoulombScattering +#/musr/command process addProcess e+ G4PenelopeIonisation, -1 2 2 +#/musr/command process addProcess e+ G4PenelopeBremsstrahlung, -1 -1 3 +#/musr/command process addProcess e+ G4PenelopeAnnihilation, 0 -1 4 +# +# --- Muons --- +#/musr/command process addProcess mu+ G4MultipleScattering -1 1 1 +##/musr/command process addProcess mu+ MultipleAndCoulombScattering -1 1 1 goulombRegion +##/musr/command process addDiscreteProcess mu+ G4CoulombScattering +#/musr/command process addProcess mu+ G4MuIonisation -1 2 2 +#/musr/command process addProcess mu+ G4MuBremsstrahlung -1 3 3 +#/musr/command process addProcess mu+ G4MuPairProduction -1 4 4 +#/musr/command process addProcess mu- G4MultipleScattering -1 1 1 +##/musr/command process addDiscreteProcess mu- G4CoulombScattering +#/musr/command process addProcess mu- G4MuIonisation -1 2 2 +#/musr/command process addProcess mu- G4MuBremsstrahlung -1 3 3 +#/musr/command process addProcess mu- G4MuPairProduction -1 4 4 +################################################################################### +################## S O M E O T H E R P A R A M E T E R S ################## +################################################################################### +# Store all events into the ROOT tree or just the interesting ones ? (true is default) +#/musr/command storeOnlyEventsWithHits false +/musr/command storeOnlyEventsWithHitInDetID 11 +# Set the minimum time separation between two subsequent signals in the same detector (in ns) +/musr/command signalSeparationTime 50 +# +/musr/run/howOftenToPrintEvent 10000 +/musr/run/randomOption 2 +################################################################################### +######################### R O O T O U T P U T ############################## +################################################################################### +#/musr/command rootOutput runID off +#/musr/command rootOutput eventID off +/musr/command rootOutput weight off +/musr/command rootOutput BFieldAtDecay off +/musr/command rootOutput muIniPosX off +/musr/command rootOutput muIniPosY off +/musr/command rootOutput muIniPosZ off +/musr/command rootOutput muIniMomX off +/musr/command rootOutput muIniMomY off +/musr/command rootOutput muIniMomZ off +/musr/command rootOutput muIniPolX off +/musr/command rootOutput muIniPolY off +/musr/command rootOutput muIniPolZ off +/musr/command rootOutput muIniTime off +/musr/command rootOutput muDecayDetID off +/musr/command rootOutput muDecayPosX off +/musr/command rootOutput muDecayPosY off +/musr/command rootOutput muDecayPosZ off +/musr/command rootOutput muDecayTime off +/musr/command rootOutput muDecayPolX off +/musr/command rootOutput muDecayPolY off +/musr/command rootOutput muDecayPolZ off +/musr/command rootOutput muTargetTime off +/musr/command rootOutput muTargetPolX off +/musr/command rootOutput muTargetPolY off +/musr/command rootOutput muTargetPolZ off +/musr/command rootOutput muM0Time off +/musr/command rootOutput muM0PolX off +/musr/command rootOutput muM0PolY off +/musr/command rootOutput muM0PolZ off +/musr/command rootOutput muM1Time off +/musr/command rootOutput muM1PolX off +/musr/command rootOutput muM1PolY off +/musr/command rootOutput muM1PolZ off +/musr/command rootOutput muM2Time off +/musr/command rootOutput muM2PolX off +/musr/command rootOutput muM2PolY off +/musr/command rootOutput muM2PolZ off +#/musr/command rootOutput posIniMomX off +#/musr/command rootOutput posIniMomY off +#/musr/command rootOutput posIniMomZ off +/musr/command rootOutput fieldNomVal off +#/musr/command rootOutput det_ID off +#/musr/command rootOutput det_edep off +/musr/command rootOutput det_edep_el off +/musr/command rootOutput det_edep_pos off +/musr/command rootOutput det_edep_gam off +/musr/command rootOutput det_edep_mup off +#/musr/command rootOutput det_nsteps off +#/musr/command rootOutput det_length off +#/musr/command rootOutput det_start off +#/musr/command rootOutput det_end off +#/musr/command rootOutput det_x off +#/musr/command rootOutput det_y off +#/musr/command rootOutput det_z off +#/musr/command rootOutput det_kine off +/musr/command rootOutput det_VrtxKine off +/musr/command rootOutput det_VrtxX off +/musr/command rootOutput det_VrtxY off +/musr/command rootOutput det_VrtxZ off +/musr/command rootOutput det_VrtxVolID off +/musr/command rootOutput det_VrtxProcID off +/musr/command rootOutput det_VrtxTrackID off +/musr/command rootOutput det_VrtxParticleID off +/musr/command rootOutput det_VvvKine off +/musr/command rootOutput det_VvvX off +/musr/command rootOutput det_VvvY off +/musr/command rootOutput det_VvvZ off +/musr/command rootOutput det_VvvVolID off +/musr/command rootOutput det_VvvProcID off +/musr/command rootOutput det_VvvTrackID off +/musr/command rootOutput det_VvvParticleID off +### Root variables that are not written out by default, but can be switched on: +#/musr/command rootOutput fieldIntegralBx on +#/musr/command rootOutput fieldIntegralBy on +#/musr/command rootOutput fieldIntegralBz on +#/musr/command rootOutput fieldIntegralBz1 on +#/musr/command rootOutput fieldIntegralBz2 on +#/musr/command rootOutput fieldIntegralBz3 on +################################################################################### +######################### V I S U A L I S A T I O N ############################## +################################################################################### +#/vis/disable +#/control/execute visFromToni.mac +/control/execute visDawn1.mac +#/control/execute visVRML.mac +################################################################################### +######################### P A R T I C L E G U N ################################# +################################################################################### +#/gps/particle ion +#/gps/ion 39 86 +#/gps/ion 27 57 0 0 +/gps/ion 38 90 0 0 +# /gps/position seems to be in cm !!!! +/gps/position 0 0 0 +/gps/energy 0 keV +#/gps/ang/maxtheta 2 deg +#/gps/ang/maxphi 2 deg +######################## B E A M O N ####################################### +/run/beamOn 1000000 diff --git a/run/visDawn101.mac b/run/visDawn101.mac new file mode 100644 index 0000000..a314f71 --- /dev/null +++ b/run/visDawn101.mac @@ -0,0 +1,70 @@ +# This is a macro file for visualizing G4 events. +# It can either be included in another macro or called with /control/exec vis.mac + +# Create an OpenGL driver (i.e. a scene handler and viewer) +# Some useful choices: VRML2FILE, OGLSX, OGLIX, DAWNFILE, etc. +#/vis/open VRML2FILE +#*/vis/open OGLIX 600x600-0+0 +/vis/open DAWNFILE + +# To calculate volumes and masses uncomment the next two lines +#*/vis/open ATree +#*/vis/ASCIITree/verbose 4 + + +# Create a new empty scene and attach it to handler +/vis/scene/create + +# Add world volume, trajectories and hits to the scene +/vis/scene/add/volume +/vis/scene/add/trajectories +/vis/scene/add/hits +/vis/sceneHandler/attach + +# Configure the viewer (optional) +#/vis/viewer/set/viewpointThetaPhi 235 -45 +/vis/viewer/set/viewpointThetaPhi 80 190 +#/vis/viewer/set/lightsThetaPhi 120 60 +#/vis/viewer/set/hiddenEdge true +/vis/viewer/set/style surface +/vis/viewer/zoom 0.8 +# Style: s - surface, w - wireframe +# Note: "surface style" and "hiddenEdge true" remove transparency! +# Other viewpoints (25 55) (235 -45) (125 35) + + +# Store trajectory information for visualisation (set to 0 if too many tracks cause core dump) +/tracking/storeTrajectory 1 + +#At the end of each event (default behaviour) +#/vis/scene/endOfEventAction refresh +#At the end of run of X events - Data from X events will be superimposed +#cks +#/vis/scene/endOfEventAction accumulate +#At the end of Y runs - Data from Y runs will be superimposed +#/vis/scene/endOfRunAction accumulate + +# Coloured trajectories for an easier particle identification: +# PDG IDs and colours: e- 11 red, e+ -11 blue, nu_e 12 yellow, +# mu+ -13 magenta, anti_nu_mu -14 green, gamma 22 grey +# +#/vis/modeling/trajectories/create/drawByCharge +#/vis/modeling/trajectories/drawByCharge-0/set 1 cyan + +/vis/modeling/trajectories/create/drawByParticleID +#*/vis/modeling/trajectories/drawByParticleID-0/set gamma grey +#/vis/modeling/trajectories/drawByParticleID-0/setRGBA gamma 1 1 1 0 +/vis/modeling/trajectories/drawByParticleID-0/setRGBA mu+ 1 0 0 1 +/vis/modeling/trajectories/drawByParticleID-0/setRGBA e+ 0 0 1 1 +/vis/modeling/trajectories/drawByParticleID-0/setRGBA gamma 0 1 0 1 +/vis/modeling/trajectories/drawByParticleID-0/setRGBA e- 1 0 1 1 +/vis/modeling/trajectories/drawByParticleID-0/setRGBA nu_e 1 1 1 0 1 +/vis/modeling/trajectories/drawByParticleID-0/setRGBA anti_nu_mu 1 1 1 0.5 +#/vis/modeling/trajectories/drawByParticleID-0/set nu_e white +#/vis/modeling/trajectories/drawByParticleID-0/set anti_nu_mu white + +# Verbosity of hits +#/hits/verbose 2 + +# Output just the detector geometry +/vis/viewer/flush diff --git a/src/G4EqEMFieldWithSpin.cc_for_Geant4.9.2p02_only b/src/G4EqEMFieldWithSpin.cc_for_Geant4.9.2p02_only new file mode 100644 index 0000000..41a6f3b --- /dev/null +++ b/src/G4EqEMFieldWithSpin.cc_for_Geant4.9.2p02_only @@ -0,0 +1,149 @@ +// +// ******************************************************************** +// * 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: G4EqEMFieldWithSpin.cc,v 1.8 2009/11/06 22:31:35 gum Exp $ +// GEANT4 tag $Name: geant4-09-03 $ +// +// +// This is the standard right-hand side for equation of motion. +// +// 30.08.2007 Chris Gong, Peter Gumplinger +// 14.02.2009 Kevin Lynch +// 06.11.2009 Hiromi Iinuma see: +// http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/161.html +// +// ------------------------------------------------------------------- + +#include "G4EqEMFieldWithSpin.hh" +#include "G4ElectroMagneticField.hh" +#include "G4ThreeVector.hh" +#include "globals.hh" + +G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField ) + : G4EquationOfMotion( emField ) +{ + anomaly = 0.0011659208; +} + +G4EqEMFieldWithSpin::~G4EqEMFieldWithSpin() +{ +} + +void +G4EqEMFieldWithSpin::SetChargeMomentumMass(G4double particleCharge, // e+ units + G4double MomentumXc, + G4double particleMass) +{ + fElectroMagCof = eplus*particleCharge*c_light ; + fMassCof = particleMass*particleMass ; + + omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss); + + ParticleCharge = particleCharge; + + E = std::sqrt(sqr(MomentumXc)+sqr(particleMass)); + beta = MomentumXc/E; + gamma = E/particleMass; + +} + +void +G4EqEMFieldWithSpin::EvaluateRhsGivenB(const G4double y[], + const G4double Field[], + G4double dydx[] ) const +{ + + // Components of y: + // 0-2 dr/ds, + // 3-5 dp/ds - momentum derivatives + // 9-11 dSpin/ds = (1/beta) dSpin/dt - spin derivatives + + // The BMT equation, following J.D.Jackson, Classical + // Electrodynamics, Second Edition, + // dS/dt = (e/mc) S \cross + // [ (g/2-1 +1/\gamma) B + // -(g/2-1)\gamma/(\gamma+1) (\beta \cdot B)\beta + // -(g/2-\gamma/(\gamma+1) \beta \cross E ] + // where + // S = \vec{s}, where S^2 = 1 + // B = \vec{B} + // \beta = \vec{\beta} = \beta \vec{u} with u^2 = 1 + // E = \vec{E} + + G4double pSquared = y[3]*y[3] + y[4]*y[4] + y[5]*y[5] ; + + G4double Energy = std::sqrt( pSquared + fMassCof ); + G4double cof2 = Energy/c_light ; + + G4double pModuleInverse = 1.0/std::sqrt(pSquared) ; + + G4double inverse_velocity = Energy * pModuleInverse / c_light; + + G4double cof1 = fElectroMagCof*pModuleInverse ; + + dydx[0] = y[3]*pModuleInverse ; + dydx[1] = y[4]*pModuleInverse ; + dydx[2] = y[5]*pModuleInverse ; + + dydx[3] = cof1*(cof2*Field[3] + (y[4]*Field[2] - y[5]*Field[1])) ; + + dydx[4] = cof1*(cof2*Field[4] + (y[5]*Field[0] - y[3]*Field[2])) ; + + dydx[5] = cof1*(cof2*Field[5] + (y[3]*Field[1] - y[4]*Field[0])) ; + + dydx[6] = dydx[8] = 0.;//not used + + // Lab Time of flight + dydx[7] = inverse_velocity; + + G4ThreeVector BField(Field[0],Field[1],Field[2]); + G4ThreeVector EField(Field[3],Field[4],Field[5]); + + EField /= c_light; + + G4ThreeVector u(y[3], y[4], y[5]); + u *= pModuleInverse; + + G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u); + G4double ucb = (anomaly+1./gamma)/beta; + G4double uce = anomaly + 1./(gamma+1.); + + G4ThreeVector Spin(y[9],y[10],y[11]); + + G4ThreeVector dSpin + = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u)) + // from Jackson + // -uce*Spin.cross(u.cross(EField)) ); + // but this form has one less operation + - uce*(u*(Spin*EField) - EField*(Spin*u)) ); + + dydx[ 9] = dSpin.x(); + dydx[10] = dSpin.y(); + dydx[11] = dSpin.z(); + + return ; +} + diff --git a/src/musrRootOutput.cc b/src/musrRootOutput.cc index 145dbae..515a70f 100644 --- a/src/musrRootOutput.cc +++ b/src/musrRootOutput.cc @@ -249,6 +249,7 @@ void musrRootOutput::BeginOfRunAction() { rootTree->Branch("save_detID",&save_detID,"save_detID[save_n]/I"); rootTree->Branch("save_particleID",&save_particleID,"save_particleID[save_n]/I"); rootTree->Branch("save_ke",&save_ke,"save_ke[save_n]/D"); + rootTree->Branch("save_time",&save_time,"save_time[save_n]/D"); rootTree->Branch("save_x",&save_x,"save_x[save_n]/D"); rootTree->Branch("save_y",&save_y,"save_y[save_n]/D"); rootTree->Branch("save_z",&save_z,"save_z[save_n]/D"); @@ -390,7 +391,8 @@ G4int musrRootOutput::ConvertProcessToID(std::string processName) { //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void musrRootOutput::SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, - G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) { + G4double x, G4double y, G4double z, G4double time, + G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) { if (save_n>=save_nMax) { char message[200]; sprintf(message,"musrRootOutput.cc::SetSaveDetectorInfo(): more \"save\" hits then allowed: save_nMax=%i",save_nMax); @@ -403,6 +405,7 @@ void musrRootOutput::SetSaveDetectorInfo (G4int ID, G4int particleID, G4double k save_x[save_n]=x/mm; save_y[save_n]=y/mm; save_z[save_n]=z/mm; + save_time[save_n]=time/microsecond; save_px[save_n]=px/MeV; save_py[save_n]=py/MeV; save_pz[save_n]=pz/MeV; diff --git a/src/musrSteppingAction.cc b/src/musrSteppingAction.cc index 8ac57ca..aae160c 100644 --- a/src/musrSteppingAction.cc +++ b/src/musrSteppingAction.cc @@ -205,13 +205,14 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) { G4double x_save=aStep->GetPreStepPoint()->GetPosition().x(); G4double y_save=aStep->GetPreStepPoint()->GetPosition().y(); G4double z_save=aStep->GetPreStepPoint()->GetPosition().z(); + G4double time_save=aStep->GetPreStepPoint()->GetGlobalTime(); G4double px_save=aStep->GetPreStepPoint()->GetMomentum().x(); G4double py_save=aStep->GetPreStepPoint()->GetMomentum().y(); G4double pz_save=aStep->GetPreStepPoint()->GetMomentum().z(); G4double polx_save=aStep->GetPreStepPoint()->GetPolarization().x(); G4double poly_save=aStep->GetPreStepPoint()->GetPolarization().y(); G4double polz_save=aStep->GetPreStepPoint()->GetPolarization().z(); - myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,px_save,py_save,pz_save,polx_save,poly_save,polz_save); + myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,time_save,px_save,py_save,pz_save,polx_save,poly_save,polz_save); } } } @@ -266,6 +267,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) { CoordinateForFieldIntegral[2] = position_tmp.z(); CoordinateForFieldIntegral[3] = aTrack->GetGlobalTime(); F04GlobalField::getObject() -> GetFieldValue(CoordinateForFieldIntegral,FieldForFieldIntegral); + // G4cout<<"B=("<GetStepLength(); // BxIntegral += stepLength; // ByIntegral += stepLength; diff --git a/src/musrTabulatedElementField.cc b/src/musrTabulatedElementField.cc index 4ecc46c..157ed4c 100644 --- a/src/musrTabulatedElementField.cc +++ b/src/musrTabulatedElementField.cc @@ -31,6 +31,7 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons { // The following posibilities of the format of the field map are distinguieshed: // 3DBOpera = 3D ... 3D field , magnetic, Opera format (Kamil-like) + // 3DEOpera = 3D ... 3D field , electric, Opera format (Kamil-like) // 3DE ... 3D field , electric, Toni-like (3DE WAS TESTED) // 3DB ... 3D field , magnetic, Toni-like // 2DBOpera = 2D .... 2D field, magnetic, Opera format (Kamil-like), X and Z components (2D WAS TESTED) @@ -54,6 +55,7 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons } symmetryType=0; + strcpy(variableIncreasingOrder,"zyx"); jx=-1; jy=0; jz=0; fldDim = (fieldTableType[0]=='3') ? 3:2; if (fldDim==2) ny=1; @@ -108,7 +110,7 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons file.ignore(256, '\n'); } while (file.peek() == '%'); } - else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)) { + else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)||(strcmp(fieldTableType,"3DEOpera")==0)) { // OPERA format of the input file: // Read table dimensions lenUnit = 1*m; @@ -142,6 +144,10 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons G4cout<<" of the Cartesian coord. system according to symmetryType = "<> xval >> yval >> zval >> bx >> by >> bz; } } - else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)) { + else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)||(strcmp(fieldTableType,"3DEOpera")==0)) { file >> xval >> yval >> zval >> bx >> by >> bz >> permeability; + // G4cout<< xval <<" "<< yval <<" "<< zval <<" "<< bx <<" "<< by <<" "<< bz <> xval >> zval >> bx >> bz; @@ -287,9 +293,34 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons minimumy = yval * lenUnit; minimumz = zval * lenUnit; } - xField[ix][iy][iz] = bx*fieldNormalisation; - yField[ix][iy][iz] = by*fieldNormalisation; - zField[ix][iy][iz] = bz*fieldNormalisation; + if (strcmp(variableIncreasingOrder,"xyz")==0) { + // The order of how the x,y,z variables increase in the field map file is non-standard. + // Recalculate the indexes of the fieldmap array for this special case: + jx++; + if (jx==nx) { + jx=0; jy++; + if (jy==ny) { + jy=0; jz++; + if (jz==nz) { // this should never happen - problem! + musrErrorMessage::GetInstance()->musrError(FATAL,"musrTabulatedElementField: too many lines found in the field map file !",false); + } + } + } + xField[jx][jy][jz] = bx*fieldNormalisation; + yField[jx][jy][jz] = by*fieldNormalisation; + zField[jx][jy][jz] = bz*fieldNormalisation; + // if ((jx==nx-1)&&(jy==ny-1)&&(jz==nz-1)) { + // G4cout<<" jx="<