Implementation of Spin Rotator.

This commit is contained in:
shiroka
2008-12-22 17:53:30 +00:00
parent 4c7938ed6d
commit 76c7f4062e
111 changed files with 409662 additions and 0 deletions

View File

@ -0,0 +1,596 @@
# Macro file for sr1.cc - Construct detector, set fields and other parameters.
# Last modified by T. Shiroka: 17.03.2008
# How to run: sr1 10xx.mac (append "idle" for prompt after running)
# sr1 10xx.mac > fname.txt (stores output on a txt file)
###############################################################################################
# #
# Specify the geometry parameters in this file (all dimensions in mm) #
# a. Lines starting with hash marks "#" are comments #
# b Lines starting with #* are temporary comments. Remove/modify to change the configuration #
# c. Lines starting with /sr1/command are commands for the executable program #
# d. Lines starting with /vis, /gun, etc. are common macro commands #
# e. Beam-line components are ordered from exp. area (MCP2) to trigger detector (TD) #
#---------------------------------------------------------------------------------------------#
# Syntax example (following /sr1/command): #
# construct solid_type volume_name parameters_defining_solid material position mothers_name #
# (mothers_name starts with log_) #
###############################################################################################
# For the meaning of the acronyms see also the original G3 file ugeom.F at:
# http://savannah.psi.ch/viewcvs/trunk/simulation/geant3/src/lemsr/ugeom.F?root=nemu%2Flem&rev=2964&view=markup
################################################################################################################
# -- ROTATION MATRICES --
################################################################################################################
# 3 parameters -> Define Euler angles (the 4th par. is set to zero).
# 4 parameters -> Define axis + rotation.
# HEP computations ordinarily use the active rotation viewpoint (object is rotated NOT axes).
# Therefore, rotations about an axis imply ACTIVE COUNTER-CLOCKWISE rotation in this package.
# Rotation around a specified axis means counter-clockwise rot. around the positive direction of the axis.
# Define rotations for the field maps of Trigger and Ring Anode:
/sr1/command rotation rotTrig 0 1 0 -45
/sr1/command rotation rotRAnR 0 0 1 -90
/sr1/command rotation rotRAnL 0 0 1 90
/sr1/command rotation rotRAnD 0 0 1 180
################################################################################################################
# -- LEM GEOMETRY --
################################################################################################################
# WORLD = Laboratory reference frame
/sr1/command construct box World 250 250 2250 G4_Galactic 0 0 0 no_logical_volume norot dead -1
# MINIMUM WORD HALF LENGTH 1250 mm!
#/sr1/command construct box World 2000 2000 4000 G4_Galactic 0 0 0 no_logical_volume norot dead -1
# World visual attributes (optional)
/sr1/command visattributes log_World invisible
#===============================================================================================================
# Sc - Scintillators: U - up, D - down, L - left, R - right (with respect to muon's view - momentum direction)
# 8 Scintillators in two concentric rings - Inner and Outer (see also the convention for the Ring Anode)
#===============================================================================================================
## Inner Scintillators - I
#*/sr1/command construct tubs ScIU 90 95 130 45 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 011 nofield
#*/sr1/command construct tubs ScIR 90 95 130 135 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 012 nofield
#*/sr1/command construct tubs ScID 90 95 130 225 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 013 nofield
#*/sr1/command construct tubs ScIL 90 95 130 315 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 014 nofield
## Outer Scintillators - O
/sr1/command construct tubs ScOU 96 101 130 45 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 021 nofield
/sr1/command construct tubs ScOR 96 101 130 135 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 022 nofield
/sr1/command construct tubs ScOD 96 101 130 225 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 023 nofield
/sr1/command construct tubs ScOL 96 101 130 315 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 024 nofield
# Visual attributes (optional)
#*/sr1/command visattributes log_ScOU SCINT_style
#*/sr1/command visattributes log_ScOD dSCINT_style
#*/sr1/command visattributes log_ScOL darkred
#===============================================================================================================
# Experimental Area - Can host EITHER the Cryostat OR the MCP2 (For the tests we usually use the MCP)
# Delimited by the F160 and F100 (blank end) flanges
#===============================================================================================================
# MCP - Multi-Channel Plate 2 Chamber; V - Vacuum, S - Solid # Note: VERY IMPORTANT: mcpv_z = -92.5 mm!
# OLD way of assigning a field
#/sr1/command construct tubs MCPV 0 76.5 254.5 0 360 G4_Galactic 0 0 -92.5 log_World norot dead 100 MCPSfield
/sr1/command construct tubs MCPV 0 76.5 254.5 0 360 G4_Galactic 0 0 -92.5 log_World norot dead 100 nofield
/sr1/command construct tubs MCPS 76.5 79.5 162.0 0 360 Steel 0 0 0 log_World norot dead 101 nofield
# F - Flanges: F160, F100, F200 (used only when the Asymmetry check is OFF)
# F160 - 160 CF flange upstream of MCP2 tube
# F100 (Blank end flange) # OLD Value was 162.0
/sr1/command construct tubs F160 79.5 101.25 11 0 360 Steel 0 0 -151.0 log_World norot dead 901 nofield
/sr1/command construct tubs F100 0 76.5 10 0 360 Steel 0 0 172.0 log_World norot dead 902 nofield
# NOTE: Original F100 referred to MCPV (as shown below) was moved to World.
#/sr1/command construct tubs F100 0 76.5 10 0 360 Steel 0 0 264.5 log_MCPV norot dead 902 nofield
# Experimental Area visual attributes (optional)
/sr1/command visattributes log_MCPV invisible
/sr1/command visattributes log_MCPS invisible
/sr1/command visattributes log_F160 blue_style
/sr1/command visattributes log_F100 blue_style
#===============================================================================================================
# MCP - Micro Channel Plate Detector MCP2 (Used as an alternative to cryostat) # mcpv_z = -92.5 mm!
#===============================================================================================================
# MCPM1 - MCP Macor ring 1
# MCPD - electron multiplying glass disk (also known as target) # Sensitive surface at z = 14 mm wrt. World
# MCPM2 - MCP Macor ring 2
#*/sr1/command construct tubs MCPM1 24 32.5 0.75 0 360 Macor 0 0 105.75 log_MCPV norot dead 201 nofield
# Use it either as (DMCP-sr1/ScintSD) - no info on mu+ polariz., or as (target-dead) with info on mu+ polariz.
#*/sr1/command construct tubs target 0 25.0 1.50 0 360 MCPglass 0 0 108.0 log_MCPV norot dead 032 nofield
/sr1/command construct tubs MCPM2 24 32.5 0.75 0 360 Macor 0 0 110.25 log_MCPV norot dead 203 nofield
# NOTE: To intercept ALL the incoming muons, comment the DMCP and MCPM1 lines above and uncomment this one:
#*/sr1/command construct tubs DMCP 0 76.5 1.5 0 360 MCPglass 0 0 108 log_MCPV norot sr1/ScintSD 202 nofield
/sr1/command construct tubs target 0 76.5 1.5 0 360 MCPglass 0 0 108 log_MCPV norot sr1/ScintSD 202 nofield
# MCSR - Stainless Steel Ring for MCP2 mounting (modelled as a box with a circular hole)
# MCVR - "Vacuum Ring" (circular hole)
/sr1/command construct box MCSR 36.5 36.5 1 Steel 0 0 112.5 log_MCPV norot dead 204 nofield
/sr1/command construct tubs MCVR 0 27.5 1 0 360 G4_Galactic 0 0 0 log_MCSR norot dead 205 nofield
# MCPA = MCP Anode (modelled as a box with two symmetrically subtracted "vacuum" disks)
# ANVA1 - Anode "Vacuum" 1 - Part of MCP Anode
# ANVA2 - Anode "Vacuum" 2 - Part of MCP Anode
/sr1/command construct box MCPA 36.5 36.5 4 Steel 0 0 123.5 log_MCPV norot dead 206 nofield
/sr1/command construct tubs ANVA1 0 27.5 1.5 0 360 G4_Galactic 0 0 -2.5 log_MCPA norot dead 207 nofield
/sr1/command construct tubs ANVA2 0 27.5 1.5 0 360 G4_Galactic 0 0 2.5 log_MCPA norot dead 208 nofield
# MCSS - MCP Stainless Steel Support Ring
/sr1/command construct tubs MCSS 40 48 2.5 0 360 Steel 0 0 156.3 log_MCPV norot dead 209 nofield
# MCP2 visual attributes (optional)
#/sr1/command visattributes log_DMCP MCP_style
#*/sr1/command visattributes log_target MCP_style
#*/sr1/command visattributes log_MCPM1 MACOR_style
/sr1/command visattributes log_MCPM2 MACOR_style
#===============================================================================================================
# CRY - Cryostat - Used as an ALTERNATIVE to MCP2 - Uncomment lines with #*. (Offset = 0.0 cm)
#===============================================================================================================
# SAH - SAmple Holder components (Cu plate) Cu or Al plates 1. Cu plate (sample holder) on Cold finger, 0.5cm
# SAPH - SAPpHire plate mounted between 1st and 2nd Cu plates, 6 mm thick, 60 mm diameter.
# SAH3 is ignored because currently NOT use.
#*/sr1/command construct tubs SAH1 0 35 2.5 0 360 G4_Cu 0 0 119.0 log_MCPV norot dead 251 nofield
#*/sr1/command construct tubs SAH2 0 35 2 0 360 G4_Cu 0 0 108.5 log_MCPV norot dead 252 nofield
#/sr1/command construct tubs SAH3 20 35 0.5 0 360 G4_Cu 0 0 106.0 log_MCPV norot dead 253 nofield
#*/sr1/command construct tubs SAPH 0 30 3 0 360 G4_ALUMINUM_OXIDE 0 0 113.5 log_MCPV norot dead 254 nofield
# Other components of the CRYostat (dimensions and position of CRY4 are only approx. because unknown)
# COFI - COld FInger
# CRY1 - End plate of cryostat (7 mm thick, 30 mm diameter)
# CRY2 - Heat exchanger (assuming a 10 mm opening - Original dimensions not known.) # OLD pos. 160.0
# CRY3 - Mounting ring for He-shield
# CRY4 - 2 mm thick plate for mounting ring. This is just to close the downstream side.
# CRSH - Lateral He-shield
# CRSH2- Frontal He-shield Ring
#*/sr1/command construct tubs COFI 0 27.5 5 0 360 G4_Cu 0 0 126.5 log_MCPV norot dead 261 nofield
#*/sr1/command construct tubs CRY1 0 15 3.5 0 360 G4_Cu 0 0 135.0 log_MCPV norot dead 262 nofield
#*/sr1/command construct tubs CRY2 5 15 25 0 360 G4_Cu 0 0 163.5 log_MCPV norot dead 263 nofield
#*/sr1/command construct tubs CRY3 38 47 5.5 0 360 G4_Cu 0 0 143.5 log_MCPV norot dead 264 nofield
#*/sr1/command construct tubs CRY4 15 38 1 0 360 G4_Cu 0 0 143.5 log_MCPV norot dead 265 nofield
#*/sr1/command construct tubs CRSH 47 48 45 0 360 G4_Cu 0 0 108.5 log_MCPV norot dead 266 nofield
#*/sr1/command construct tubs CRSH2 30 48 0.5 0 360 G4_Cu 0 0 63.0 log_MCPV norot dead 267 nofield
# Electrical Field Guard Rings (distance between the guard rings: 16 mm)
#*/sr1/command construct tubs Guard1 29 38 1.5 0 360 G4_Cu 0 0 76.0 log_MCPV norot dead 271 nofield
#*/sr1/command construct tubs Guard2 29 38 1.5 0 360 G4_Cu 0 0 92.0 log_MCPV norot dead 272 nofield
# Cryostat visual attributes (optional)
#*/sr1/command visattributes log_SAH1 oxsteel
#*/sr1/command visattributes log_SAH2 oxsteel
#*/sr1/command visattributes log_SAPH MACOR_style
#/sr1/command visattributes log_SAH3 oxsteel
#/sr1/command visattributes log_CRSH invisible
#===============================================================================================================
# RA - Ring Anode, M - middle part (closer to Ground Anode), E - end part (farther from the Ground Anode)
# U - up, D - down, L - left, R - right (with respect to muon's view - momentum direction)
# Note: 3.0 mm HALF gap at 45.1469 mm half radius => delta_ang = asin(3.0/45.1469)*180/pi = 3.81 deg.
# Note: delta_ang = 3.1744 deg. for 2.5 mm HG. The angular extension goes e.g. from (45 + da) to (90 - 2*da).
# Note: Ring Anode - Ground Anode distance was 15 mm => CHANGED to 12 mm! (Positions: 11.5 -> 8.5, -33.5 -> -36.5)
#===============================================================================================================
# RA_Ez = -10.35+2.25 = -8.1 cm; RA_Mz= -10.35 - 2.25 = -12.6 cm; RA_Gz= -25.45+3.75 = -21.7 cm; mcpv_z = -9.25 cm
/sr1/command construct cones RA_EU 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV norot dead 301 nofield
/sr1/command construct cones RA_MU 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV norot dead 302 nofield
/sr1/command construct cones RA_ER 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV rotRAnR dead 303 nofield
/sr1/command construct cones RA_MR 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV rotRAnR dead 304 nofield
/sr1/command construct cones RA_ED 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV rotRAnD dead 305 nofield
/sr1/command construct cones RA_MD 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV rotRAnD dead 306 nofield
/sr1/command construct cones RA_EL 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV rotRAnL dead 307 nofield
/sr1/command construct cones RA_ML 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV rotRAnL dead 308 nofield
# Dummy, thin cylindres used for applying the SAME RA field-map (ROTATED by 90 deg.) to different anodes.
# NOTE: EM field cannot be applied to non simply connected bodies, as e.g. rings, cones, tori, etc.!
/sr1/command construct tubs RA_U 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.50 log_MCPV norot dead 322 nofield
/sr1/command construct tubs RA_R 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.52 log_MCPV rotRAnR dead 324 nofield
/sr1/command construct tubs RA_D 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.54 log_MCPV rotRAnD dead 326 nofield
/sr1/command construct tubs RA_L 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.56 log_MCPV rotRAnL dead 328 nofield
# RA_G - Ring Anode Ground Cylinder
/sr1/command construct tubs RA_G 58 62.5 58.0 0 360 G4_Cu 0 0 -129.0 log_MCPV norot dead 351 nofield
# Ring Anodes visual attributes (optional)
/sr1/command visattributes log_RA_EU oxsteel
/sr1/command visattributes log_RA_MR oxsteel
/sr1/command visattributes log_RA_G Grid_style
# Alternative placement using World as a mother volume (mcpv_z = -92.5 mm). Check latter. These values refer to a 5 mm GAP!
#/sr1/command construct cones RA_EU 45.1469 62.5 33.5 39.0 22.5 48.1711 83.6578 Steel 0 0 -81 log_World norot dead 301 nofield
#/sr1/command construct cones RA_MU 56.7937 62.5 45.147 62.5 22.5 48.1711 83.6578 Steel 0 0 -126 log_World norot dead 302 nofield
#===============================================================================================================
# Gate Valve Area - Hosts the Ground Anode (upstream part of the Ring Anode) - Delimited by L3F1 and F200 flanges
#===============================================================================================================
# GATS - 200 CF flange upstream of MCP2 tube covering the whole length of the gate valve chamber. For simplicity, we
# choose the INNER diameter of the GATe valve Steel tube the same as the OUTER diameter of F200 (200 CF flange)
/sr1/command construct tubs GATS 103.25 126.5 92.5 0 360 Steel 0 0 -254.5 log_World norot dead 371 nofield
# Vacuum "Ring" (to avoid intersections with MCPV) - Not needed if world is already filled with vacuum.
#*/sr1/command construct tubs GATV 76.5 103.25 92.5 0 360 G4_Galactic 0 0 -254.5 log_World norot dead 370 nofield
# F200 - 200 CF flange upstream of MCP2 tube to connect to gate valve chamber
/sr1/command construct tubs F200 76.5 103.25 12 0 360 Steel 0 0 -174.0 log_World norot dead 372 nofield
# NOTE: When using GATV comment the F200 above and uncomment the following (change mother to log_GATV).
#*/sr1/command construct tubs F200 76.5 103.25 12 0 360 Steel 0 0 80.5 log_GATV norot dead 372 nofield
# Gate Valve Area visual attributes (optional)
/sr1/command visattributes log_GATS SCINT_style
/sr1/command visattributes log_F200 blue_style
#===============================================================================================================
# L3 - 3rd Einzel Lens # L3z = -56.7 cm. ATT: DUMMY FIELD change to electric L3Efield!
#===============================================================================================================
# Lens Gap = 12.0 mm => G/D = 12/130 ~ 0.1 (Lens Gap = gap between Ground and Anode, D - Diameter)
# L3 envelope (Tube + Flanges)
# L3VA - Lens 3 Vacuum
# L3ST - Lens 3 Steel tube (inner dia: 200 mm, outer dia: 206 mm, length: 720 mm)
# L3F1 - Lens 3 Flange 1, z = L3z + 208 mm
# L3F2 - Lens 3 Flange 2, z = L3z - 208 mm
/sr1/command construct tubs L3VA 0 100 220 0 360 G4_Galactic 0 0 -567 log_World norot dead 400 nofield
/sr1/command construct tubs L3ST 100 103 220 0 360 Steel 0 0 -567 log_World norot dead 401 nofield
/sr1/command construct tubs L3F1 103 126.5 12 0 360 Steel 0 0 -359 log_World norot dead 402 nofield
/sr1/command construct tubs L3F2 103 126.5 12 0 360 Steel 0 0 -775 log_World norot dead 403 nofield
# GPn - Ground Potential Electrodes
# n = 1-4 and 5-8 - components of the Ground Electrodes
# GP1 - Ground Electrode (inner dia: 130 mm, outer dia: 134 mm, length: 133 mm)
# GP2 - outer electrode surface (LN2 cooling vessel)
# GP3 - first ring cap
# GP4 - second ring cap
# n = 1-4 - Ground Electrode 1 (further from TD). See above for the meaning of acronyms.
/sr1/command construct tubs L3GP1 65 67 66.5 0 360 Steel 0 0 133.5 log_L3VA norot dead 421 nofield
/sr1/command construct tubs L3GP2 81 83 66.5 0 360 Steel 0 0 133.5 log_L3VA norot dead 422 nofield
/sr1/command construct tubs L3GP3 67 81 4 0 360 Steel 0 0 196.0 log_L3VA norot dead 423 nofield
/sr1/command construct tubs L3GP4 67 81 4 0 360 Steel 0 0 71.0 log_L3VA norot dead 424 nofield
# n = 5-8 - Ground Electrode 2 (closer to TD). See above for the meaning of acronyms.
/sr1/command construct tubs L3GP5 65 67 66.5 0 360 Steel 0 0 -133.5 log_L3VA norot dead 431 nofield
/sr1/command construct tubs L3GP6 81 83 66.5 0 360 Steel 0 0 -133.5 log_L3VA norot dead 432 nofield
/sr1/command construct tubs L3GP7 67 81 4 0 360 Steel 0 0 -196.0 log_L3VA norot dead 433 nofield
/sr1/command construct tubs L3GP8 67 81 4 0 360 Steel 0 0 -71.0 log_L3VA norot dead 434 nofield
# HP - High Potential Electrode (Central Anode - usually at +8.7 kV, for a 15 keV muon beam)
/sr1/command construct tubs L3HP 65 83 55 0 360 Steel 0 0 0 log_L3VA norot dead 451 nofield
# Lens 3 visual attributes (optional)
/sr1/command visattributes log_L3VA invisible
/sr1/command visattributes log_L3ST invisible
/sr1/command visattributes log_L3HP darkred
#===============================================================================================================
# IP - Intermediate Piece (between Trigger Detector and Einzel Lens 3)
# Original name was CGate - B-field Coil Compensation Gate?! # CompGatez = -86.55 cm; // L3z - 22 - 7.85 cm
#===============================================================================================================
# IPV (but also others) are just empty volumes "filled" with vacuum. Sometimes used to apply EM fields to restricted areas.
/sr1/command construct tubs IPV 0 100 78.5 0 360 G4_Galactic 0 0 -865.5 log_World norot dead 500 nofield
# IPCF - Intermediate Piece Central Flange (same as L3 Flanges)
# IPST - Intermediate Piece Steel Tube (same diameter as L3 Steel Tube)
/sr1/command construct tubs IPCF 103 126.5 12.0 0 360 Steel 0 0 -865.5 log_World norot dead 501 nofield
/sr1/command construct tubs IPST 100 103 78.5 0 360 Steel 0 0 -865.5 log_World norot dead 502 nofield
# IP visual attributes (optional)
/sr1/command visattributes log_IPV invisible
/sr1/command visattributes log_IPCF blue_style
/sr1/command visattributes log_IPST gray
#===============================================================================================================
# Trigger - Trigger Detector # Triggerz = -1092 mm
#===============================================================================================================
# Trigger tube and relative vacuum
/sr1/command construct tubs TriggerV 0 100 148 0 360 G4_Galactic 0 0 -1092 log_World norot dead 600 nofield
/sr1/command construct tubs Trigger 100 103 148 0 360 Steel 0 0 -1092 log_World norot dead 601 nofield
# TF - Trigger tube flanges
/sr1/command construct tubs TF1 103 126.5 12 0 360 Steel 0 0 -956 log_World norot dead 611 nofield
/sr1/command construct tubs TF2 103 126.5 12 0 360 Steel 0 0 -1228 log_World norot dead 612 nofield
# Carbon Foil (default thickness 0.000005147 mm, see below).
# USE THE NAME CFoil or coulombCFoil, otherwise sr1MuFormation won't work!
/sr1/command construct box CFoil 60 60 0.000005147 G4_GRAPHITE 0 0 -45 log_TriggerV norot dead 621 nofield
###/sr1/command construct box coulombCFoil 60 60 0.000005147 G4_GRAPHITE 0 0 -45 log_TriggerV norot dead 621 nofield
# Notes: NIST tables use G4_GRAPHITE with 1.7 g/cm3 and 78 eV ioniz. energy.
# An area density of 1.75 ug/cm2 implies a CF thickn. = (1.75*1.e-6/1.70)*cm = 1.029e-5 mm - Total thickness
# If necessary, use Graphite as defined in sr1DetectorConstruction.cc and set any density.
# Electrical Field areas in the Trigger Detector
# En = Electrical Field n: TnFieldMgr (n = 1-3)
# Original TriggE2: [4.*sqrt(2), 4.5, 0.7/sqrt(2)] cm -> changed due to overlaps with E1 and E3
/sr1/command construct box TriggE1 45 45 4 G4_Galactic 0 0 -38.5 log_TriggerV norot dead 631 nofield
/sr1/command construct box TriggE2 45 45 4.9497 G4_Galactic 0 0 2.25 log_TriggerV rotTrig dead 632 nofield
/sr1/command construct box TriggE3 45 45 4 G4_Galactic 0 0 43.0 log_TriggerV norot dead 633
# Beam spot (just for having a visual idea!)
/sr1/command construct tubs BSpot 0 20 1 0 360 G4_Galactic 0 0 -48.0 log_TriggerV norot dead 650 nofield
# Trigger visual attributes (optional)
/sr1/command visattributes log_TriggerV invisible
/sr1/command visattributes log_Trigger invisible
#/sr1/command visattributes log_CFoil MACOR_style
/sr1/command visattributes log_BSpot fblue_style
# One can set visible attrib. also on a MATERIAL basis, rather than on log_VOL.
# E.g. /sr1/command visattributes Steel red
################################################################################################################
# -- Setting the ELECTRIC and MAGNETIC fields --
################################################################################################################
# Use ABSOLUTE coordinates to specify the field position (i.e. with respect to GLOBAL WORLD)!
# Default field units: Magnetic - T, Electric - kV/mm (or kV for E-field maps).
# NOTE: Applying a field to an invisible log_vol makes is visible!
### Electric field at TRIGGER Detector TD: Three different uniform fields
/sr1/command globalfield Trigg1_field 0. 0. -1130.5 uniform log_TriggE1 0 0 0 0 0 -0.02375
/sr1/command globalfield Trigg2_field 0. 0. -1089.75 uniform log_TriggE2 0 0 0 0 0 0.041416
/sr1/command globalfield Trigg3_field 0. 0. -1049.0 uniform log_TriggE3 0 0 0 0 0 -0.49375
### Electric field at Einzel LENS 3 - from folded 2D axial field map (f - folded, i.e. symmetric)
# Typically V = +8.7 kV for a muon beam at 15 keV. Use either 2DE L3_Erz.map or
/sr1/command globalfield Lens3_field 0. 0. -567. fromfile 2DEf L3_Erz4.map log_L3VA 8.7
# To change the field in regular steps (e.g. for testing) use (f_min f_max step_no), e.g.:
#*/sr1/command globalfield Lens3_field 0. 0. -567. fromfile 2DEf L3_Erz4.map log_L3VA 7 11 5
### Electric field at RING ANODE - from 3DE field map
# Typically set at +11.0 kV for a muon beam at 15 keV
# To create an arbitrary configuration, switch on all fields and set different potentials.
#######/sr1/command globalfield RngAnU_field 0. 0. -167.00 fromfile 3DE EM_3D_extc.map log_RA_U 11.0
/sr1/command globalfield RngAnU_field 0. 0. -143.00 fromfile 3DE EM_3D_ext_gridf.map log_RA_U 11.0
/sr1/command globalfield RngAnR_field 0. 0. -143.02 fromfile 3DE EM_3D_ext_gridf.map log_RA_R 11.0
/sr1/command globalfield RngAnD_field 0. 0. -143.04 fromfile 3DE EM_3D_ext_gridf.map log_RA_D 11.0
/sr1/command globalfield RngAnL_field 0. 0. -143.06 fromfile 3DE EM_3D_ext_gridf.map log_RA_L 11.0
### LAST FIELD OK: EM_3D_ext2f.map
# EXTENDED MAPS (from -28 to +20 mm) EM_3D_extc.map or EM_3D_extf.map (coord + field or field only)
# "Best" results with EM_RA_3D.map, even though this is not a precise map (poor meshing).
#RA_1kV_upf_raw.map
#RA_1kV_upf.map # EM_extendf2.map give rise to "strange spots" -> sr1_1047_RA13.eps
### Electric field at SAMPLE space. Three possible field maps: Sample, G1 and G2 (closest to sample)
# To create an arbitrary configuration, switch on all fields and set different potentials.
# Field extension along z: 50 mm. Center of MCPV at z = -92.5 mm. Sample face at: z = 108.5 (4 mm thick - SAH2).
# Position of the field: (-92.5 - 50/2) + (108.5 - 4/2) = -11 mm wrt. WORLD! =>
#*/sr1/command globalfield Sample_field 0. 0. -11.0 fromfile 2DE sample_Erz.map log_MCPV 9.0
#*/sr1/command globalfield Guard2_field 0. 0. -11.0 fromfile 2DE guard2_Erz.map log_MCPV 6.0
#*/sr1/command globalfield Guard1_field 0. 0. -11.0 fromfile 2DE guard1_Erz.map log_MCPV 3.0
### Magnetic field at SAMPLE space (use either a uniform field or a field map).
# Use either DMCP or MCPV to apply a CONSTANT field strictly to the sample or also to the surroundings, resp.!
#*/sr1/command globalfield Magnet_field 0. 0. 14.5 uniform log_DMCP 0 0.005 0 0 0 0
# Use field map to set field to 20 G TF
# Extended map has -100/100 cm z extension, field center at +70, original -85/85, center at +85
#####/sr1/command globalfield Magnet_field 0. 0. -836.0 fromfile 2DB sample_Brz.map log_IPV 0.002
/sr1/command globalfield Magnet_field 0. 0. -686.0 fromfile 2DB sample_Brz_ext.map log_L3VA 0.002
# Set parameters for particle tracking in an EM field
/sr1/command globalfield setparameter SetLargestAcceptableStep 5
/sr1/command globalfield setparameter SetMinimumEpsilonStep 5e-5
/sr1/command globalfield setparameter SetMaximumEpsilonStep 0.001
/sr1/command globalfield setparameter SetDeltaOneStep 0.1
/sr1/command globalfield setparameter SetDeltaIntersection 0.01
/sr1/command globalfield printparameters
################################################################################################################
# -- Testing the ELECTRIC and MAGNETIC fields (OPTIONAL) --
################################################################################################################
# FIELD CHECKS at different beam transport components (from trigg. to sample)
# All distances in mm and in GLOBAL coordinates (preferably at field center)
# The test points below are just some examples. Any point can be checked.
# Trigger 1
/sr1/command globalfield printFieldValueAtPoint 0. 0. -1130.5
# Trigger 2
/sr1/command globalfield printFieldValueAtPoint 0. 0. -1089.75
# Trigger 3
/sr1/command globalfield printFieldValueAtPoint 0. 0. -1049.0
# Einzel Lens 3 - L3 (center at -567.0, but max field at rel. +/-62 mm)
/sr1/command globalfield printFieldValueAtPoint 0. 0. -507.0
# Ring Anode - RA (center at -167.0, but max field at rel. -16/+132 mm)
/sr1/command globalfield printFieldValueAtPoint 0. 0. -35.0
# Sample space (center at -11.0, but max field at rel. +24 mm)
/sr1/command globalfield printFieldValueAtPoint 0. 0. 13.0
# Check magnetic field at sample space
#*/sr1/command globalfield printFieldValueAtPoint 10. 0. 13.0
# Check magnetic field at the lower field end limit
#*/sr1/command globalfield printFieldValueAtPoint 0. 20. -1680.
################################################################################################################
# -- Setting simulation PARAMETERS --
################################################################################################################
# Set processes from: lowenergy, penelope, coulombAndMultiple (default), coulomb (only for CFoil).
/sr1/command typeofprocesses lowenergy
#*/sr1/command typeofprocesses penelope
#*/sr1/command includeMuoniumProcesses false
# Set the overall range cut
#*/run/setCut 1 mm
# Set the range cut on particular volumes (in mm)
#*/sr1/command SetUserLimits log_target 0.01
#*/sr1/command SetUserLimits log_targetscint 0.01
#*/sr1/command SetUserLimits log_cryostatscint 0.01
# Set particle energy cuts on particular volumes (in eV)
/sr1/command SetUserMinEkine log_World 0.1
# Store ALL the events in a ROOT tree or just the interesting ones? (default is true)
#*/sr1/command storeOnlyEventsWithHits false
# Set the minimum time separation between two subsequent signals in the same detector (in ns)
/sr1/command signalSeparationTime 0.1
# Override runID number
#*/sr1/run/runID 21
# Set the frequency of event printing
/sr1/run/howOftenToPrintEvent 1000
# RANDOM option choices: (specify the random number generator initialisation)
# 0 ... no initialisation (default)
# 1 ... use actual computer time to initialise now
# 2 ... use event number to initialise at the beginning of each event
# 3 ... read in the random no. initial values for each event from a file
/sr1/run/randomOption 2
# VISUALIZATION options
# To enable or disable visualization uncomment one of these lines
# To modify visualization options edit the file vis.mac
#/vis/disable
/control/execute vis.mac
################################################################################################################
# -- Setting PARTICLE GUN parameters --
################################################################################################################
# Default momentum direction: 001, i.e. 0z.
# Default muon spin direction: 100, i.e. 0x.
# Default particle type: mu+ (can be changed to Mu)
# Set particle type
#*/gun/particle Mu
/gun/particle mu+
# Set beam vertex (to start after C Foil use z = -1130)
/gun/vertex 0. 0. -1170. mm
# A point-like uniform beam
/gun/vertexsigma -0.1 -0.1 0 mm
# Set beam transverse spread (default GAUSSIAN spread)
# If FWHM = 10 mm ==> sigma = 10/2.354 = 4.2481 mm (last 0 is a dummy value)
# Negative sigma values => random FLAT RECTANGULAR distribution (area 2x.2y)
# Use vertexboundary with (vb < sigma_xy) to obtain a CIRCULAR beam spot
# /gun/vertexsigma 0 0 0 mm ==> Very SLOW with mag. field ON and centered beam
#*/gun/vertexsigma 42.5 42.5 0 mm
/gun/vertexsigma -20 -20 0 mm
/gun/vertexboundary 30 -1e6 1e6 mm
# /gun/vertexboundary: rMaxAllowed, zMinAllowed, zMaxAllowed # Beam AND gating
#*/gun/vertexboundary 7 -1314.4 -1305 mm
# Without restrictions in z, but only on r:
#*/gun/vertexboundary 3 -1e6 1e6 mm
# Set beam momentum (USE only as an ALTERNATIVE to setting energy!)
# /gun/momentum 0 0 29.79 MeV
#*/gun/momentum 0 0 1.8 MeV
# Energy loss at p = 1.2 MeV/c (E = 6.8 keV) => 1.23 +/- 0.2 keV
# Energy loss at p = 1.8 MeV/c (E = 15.3 keV) => 1.25 +/- 0.3 keV
# 1.2 MeV/c -> 6.8 keV, 1.8 MeV/c -> 15.3 keV
# muon rest mass = 105.658 MeV/c2
# With Trigger ON use 15 keV, with Trigger OFF use 11.27 keV (15 - 3.73 = 11.27 keV)
/gun/kenergy 15 keV
# Set beam momentum direction
/gun/direction 0.0 0.0 1.0
# Set muon spin direction
/gun/muonPolarizVector 1 0 0
# Other useful test parameters:
#
# FWHM = 3% ==> sigma = 29.79*0.03/2.354 = 0.37965 MeV/c
#*/gun/momentumsmearing 0.37965 MeV
#---/gun/momentumboundary: pMinAllowed, pMaxAllowed, dummy
#*/gun/momentumboundary 20 40 0 MeV
#---/gun/tilt: xangle, yangle, dummy
#*/gun/tilt 0 0.5 0 deg
#---/gun/tiltsigma: xangleSigma, yangleSigma, dummy (1 degree at 1 m => 17 mm)
#*/gun/tiltsigma 0.2 0.2 0 deg
#*/gun/pitch 0.5 deg
#---/gun/decaytimelimits: decayMin, decayMax, decayTime
#*/gun/decaytimelimits 10400 10420 2197.03 ns
# BEAM ON
#*/run/beamOn 1000000
/run/beamOn 2
#/run/beamOn 10000

View File

@ -0,0 +1,612 @@
# Macro file for sr1.cc - Construct detector, set fields and other parameters.
# Last modified by T. Shiroka: 17.03.2008
# How to run: sr1 10xx.mac (append "idle" for prompt after running)
# sr1 10xx.mac > fname.txt (stores output on a txt file)
###############################################################################################
# #
# Specify the geometry parameters in this file (all dimensions in mm) #
# a. Lines starting with hash marks "#" are comments #
# b Lines starting with #* are temporary comments. Remove/modify to change the configuration #
# c. Lines starting with /sr1/command are commands for the executable program #
# d. Lines starting with /vis, /gun, etc. are common macro commands #
# e. Beam-line components are ordered from exp. area (MCP2) to trigger detector (TD) #
#---------------------------------------------------------------------------------------------#
# Syntax example (following /sr1/command): #
# construct solid_type volume_name parameters_defining_solid material position mothers_name #
# (mothers_name starts with log_) #
###############################################################################################
# For the meaning of the acronyms see also the original G3 file ugeom.F at:
# http://savannah.psi.ch/viewcvs/trunk/simulation/geant3/src/lemsr/ugeom.F?root=nemu%2Flem&rev=2964&view=markup
################################################################################################################
# -- ROTATION MATRICES --
################################################################################################################
# 3 parameters -> Define Euler angles (the 4th par. is set to zero).
# 4 parameters -> Define axis + rotation.
# HEP computations ordinarily use the active rotation viewpoint (object is rotated NOT axes).
# Therefore, rotations about an axis imply ACTIVE COUNTER-CLOCKWISE rotation in this package.
# Rotation around a specified axis means counter-clockwise rot. around the positive direction of the axis.
# Define rotations for the field maps of Trigger and Ring Anode:
/sr1/command rotation rotTrig 0 1 0 -45
/sr1/command rotation rotRAnR 0 0 1 -90
/sr1/command rotation rotRAnL 0 0 1 90
/sr1/command rotation rotRAnD 0 0 1 180
################################################################################################################
# -- LEM GEOMETRY --
################################################################################################################
# WORLD = Laboratory reference frame
/sr1/command construct box World 250 250 2250 G4_Galactic 0 0 0 no_logical_volume norot dead -1
# MINIMUM WORD HALF LENGTH 1250 mm!
#/sr1/command construct box World 2000 2000 4000 G4_Galactic 0 0 0 no_logical_volume norot dead -1
# World visual attributes (optional)
/sr1/command visattributes log_World invisible
#===============================================================================================================
# Sc - Scintillators: U - up, D - down, L - left, R - right (with respect to muon's view - momentum direction)
# 8 Scintillators in two concentric rings - Inner and Outer (see also the convention for the Ring Anode)
#===============================================================================================================
## Inner Scintillators - I
#*/sr1/command construct tubs ScIU 90 95 130 45 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 011 nofield
#*/sr1/command construct tubs ScIR 90 95 130 135 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 012 nofield
#*/sr1/command construct tubs ScID 90 95 130 225 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 013 nofield
#*/sr1/command construct tubs ScIL 90 95 130 315 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 014 nofield
## Outer Scintillators - O
/sr1/command construct tubs ScOU 96 101 130 45 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 021 nofield
/sr1/command construct tubs ScOR 96 101 130 135 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 022 nofield
/sr1/command construct tubs ScOD 96 101 130 225 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 023 nofield
/sr1/command construct tubs ScOL 96 101 130 315 90 G4_PLASTIC_SC_VINYLTOLUENE 0 0 0 log_World norot sr1/ScintSD 024 nofield
# Visual attributes (optional)
#*/sr1/command visattributes log_ScOU SCINT_style
#*/sr1/command visattributes log_ScOD dSCINT_style
#*/sr1/command visattributes log_ScOL darkred
#===============================================================================================================
# Experimental Area - Can host EITHER the Cryostat OR the MCP2 (For the tests we usually use the MCP)
# Delimited by the F160 and F100 (blank end) flanges
#===============================================================================================================
# MCP - Multi-Channel Plate 2 Chamber; V - Vacuum, S - Solid # Note: VERY IMPORTANT: mcpv_z = -92.5 mm!
# OLD way of assigning a field
#/sr1/command construct tubs MCPV 0 76.5 254.5 0 360 G4_Galactic 0 0 -92.5 log_World norot dead 100 MCPSfield
/sr1/command construct tubs MCPV 0 76.5 254.5 0 360 G4_Galactic 0 0 -92.5 log_World norot dead 100 nofield
/sr1/command construct tubs MCPS 76.5 79.5 162.0 0 360 Steel 0 0 0 log_World norot dead 101 nofield
# F - Flanges: F160, F100, F200 (used only when the Asymmetry check is OFF)
# F160 - 160 CF flange upstream of MCP2 tube
# F100 (Blank end flange) # OLD Value was 162.0
/sr1/command construct tubs F160 79.5 101.25 11 0 360 Steel 0 0 -151.0 log_World norot dead 901 nofield
/sr1/command construct tubs F100 0 76.5 10 0 360 Steel 0 0 172.0 log_World norot dead 902 nofield
# NOTE: Original F100 referred to MCPV (as shown below) was moved to World.
#/sr1/command construct tubs F100 0 76.5 10 0 360 Steel 0 0 264.5 log_MCPV norot dead 902 nofield
# Experimental Area visual attributes (optional)
/sr1/command visattributes log_MCPV invisible
/sr1/command visattributes log_MCPS invisible
/sr1/command visattributes log_F160 blue_style
/sr1/command visattributes log_F100 blue_style
#===============================================================================================================
# MCP - Micro Channel Plate Detector MCP2 (Used as an alternative to cryostat) # mcpv_z = -92.5 mm!
#===============================================================================================================
# MCPM1 - MCP Macor ring 1
# MCPD - electron multiplying glass disk (also known as target) # Sensitive surface at z = 14 mm wrt. World
# MCPM2 - MCP Macor ring 2
#*/sr1/command construct tubs MCPM1 24 32.5 0.75 0 360 Macor 0 0 105.75 log_MCPV norot dead 201 nofield
# Use it either as (DMCP-sr1/ScintSD) - no info on mu+ polariz., or as (target-dead) with info on mu+ polariz.
#*/sr1/command construct tubs target 0 25.0 1.50 0 360 MCPglass 0 0 108.0 log_MCPV norot dead 032 nofield
/sr1/command construct tubs MCPM2 24 32.5 0.75 0 360 Macor 0 0 110.25 log_MCPV norot dead 203 nofield
# NOTE: To intercept ALL the incoming muons, comment the DMCP and MCPM1 lines above and uncomment this one:
#*/sr1/command construct tubs DMCP 0 76.5 1.5 0 360 MCPglass 0 0 108 log_MCPV norot sr1/ScintSD 202 nofield
/sr1/command construct tubs target 0 76.5 1.5 0 360 MCPglass 0 0 108 log_MCPV norot sr1/ScintSD 202 nofield
# MCSR - Stainless Steel Ring for MCP2 mounting (modelled as a box with a circular hole)
# MCVR - "Vacuum Ring" (circular hole)
/sr1/command construct box MCSR 36.5 36.5 1 Steel 0 0 112.5 log_MCPV norot dead 204 nofield
/sr1/command construct tubs MCVR 0 27.5 1 0 360 G4_Galactic 0 0 0 log_MCSR norot dead 205 nofield
# MCPA = MCP Anode (modelled as a box with two symmetrically subtracted "vacuum" disks)
# ANVA1 - Anode "Vacuum" 1 - Part of MCP Anode
# ANVA2 - Anode "Vacuum" 2 - Part of MCP Anode
/sr1/command construct box MCPA 36.5 36.5 4 Steel 0 0 123.5 log_MCPV norot dead 206 nofield
/sr1/command construct tubs ANVA1 0 27.5 1.5 0 360 G4_Galactic 0 0 -2.5 log_MCPA norot dead 207 nofield
/sr1/command construct tubs ANVA2 0 27.5 1.5 0 360 G4_Galactic 0 0 2.5 log_MCPA norot dead 208 nofield
# MCSS - MCP Stainless Steel Support Ring
/sr1/command construct tubs MCSS 40 48 2.5 0 360 Steel 0 0 156.3 log_MCPV norot dead 209 nofield
# MCP2 visual attributes (optional)
#/sr1/command visattributes log_DMCP MCP_style
#*/sr1/command visattributes log_target MCP_style
#*/sr1/command visattributes log_MCPM1 MACOR_style
/sr1/command visattributes log_MCPM2 MACOR_style
#===============================================================================================================
# CRY - Cryostat - Used as an ALTERNATIVE to MCP2 - Uncomment lines with #*. (Offset = 0.0 cm)
#===============================================================================================================
# SAH - SAmple Holder components (Cu plate) Cu or Al plates 1. Cu plate (sample holder) on Cold finger, 0.5cm
# SAPH - SAPpHire plate mounted between 1st and 2nd Cu plates, 6 mm thick, 60 mm diameter.
# SAH3 is ignored because currently NOT use.
#*/sr1/command construct tubs SAH1 0 35 2.5 0 360 G4_Cu 0 0 119.0 log_MCPV norot dead 251 nofield
#*/sr1/command construct tubs SAH2 0 35 2 0 360 G4_Cu 0 0 108.5 log_MCPV norot dead 252 nofield
#/sr1/command construct tubs SAH3 20 35 0.5 0 360 G4_Cu 0 0 106.0 log_MCPV norot dead 253 nofield
#*/sr1/command construct tubs SAPH 0 30 3 0 360 G4_ALUMINUM_OXIDE 0 0 113.5 log_MCPV norot dead 254 nofield
# Other components of the CRYostat (dimensions and position of CRY4 are only approx. because unknown)
# COFI - COld FInger
# CRY1 - End plate of cryostat (7 mm thick, 30 mm diameter)
# CRY2 - Heat exchanger (assuming a 10 mm opening - Original dimensions not known.) # OLD pos. 160.0
# CRY3 - Mounting ring for He-shield
# CRY4 - 2 mm thick plate for mounting ring. This is just to close the downstream side.
# CRSH - Lateral He-shield
# CRSH2- Frontal He-shield Ring
#*/sr1/command construct tubs COFI 0 27.5 5 0 360 G4_Cu 0 0 126.5 log_MCPV norot dead 261 nofield
#*/sr1/command construct tubs CRY1 0 15 3.5 0 360 G4_Cu 0 0 135.0 log_MCPV norot dead 262 nofield
#*/sr1/command construct tubs CRY2 5 15 25 0 360 G4_Cu 0 0 163.5 log_MCPV norot dead 263 nofield
#*/sr1/command construct tubs CRY3 38 47 5.5 0 360 G4_Cu 0 0 143.5 log_MCPV norot dead 264 nofield
#*/sr1/command construct tubs CRY4 15 38 1 0 360 G4_Cu 0 0 143.5 log_MCPV norot dead 265 nofield
#*/sr1/command construct tubs CRSH 47 48 45 0 360 G4_Cu 0 0 108.5 log_MCPV norot dead 266 nofield
#*/sr1/command construct tubs CRSH2 30 48 0.5 0 360 G4_Cu 0 0 63.0 log_MCPV norot dead 267 nofield
# Electrical Field Guard Rings (distance between the guard rings: 16 mm)
#*/sr1/command construct tubs Guard1 29 38 1.5 0 360 G4_Cu 0 0 76.0 log_MCPV norot dead 271 nofield
#*/sr1/command construct tubs Guard2 29 38 1.5 0 360 G4_Cu 0 0 92.0 log_MCPV norot dead 272 nofield
# Cryostat visual attributes (optional)
#*/sr1/command visattributes log_SAH1 oxsteel
#*/sr1/command visattributes log_SAH2 oxsteel
#*/sr1/command visattributes log_SAPH MACOR_style
#/sr1/command visattributes log_SAH3 oxsteel
#/sr1/command visattributes log_CRSH invisible
#===============================================================================================================
# RA - Ring Anode, M - middle part (closer to Ground Anode), E - end part (farther from the Ground Anode)
# U - up, D - down, L - left, R - right (with respect to muon's view - momentum direction)
# Note: 3.0 mm HALF gap at 45.1469 mm half radius => delta_ang = asin(3.0/45.1469)*180/pi = 3.81 deg.
# Note: delta_ang = 3.1744 deg. for 2.5 mm HG. The angular extension goes e.g. from (45 + da) to (90 - 2*da).
# Note: Ring Anode - Ground Anode distance was 15 mm => CHANGED to 12 mm! (Positions: 11.5 -> 8.5, -33.5 -> -36.5)
#===============================================================================================================
# RA_Ez = -10.35+2.25 = -8.1 cm; RA_Mz= -10.35 - 2.25 = -12.6 cm; RA_Gz= -25.45+3.75 = -21.7 cm; mcpv_z = -9.25 cm
/sr1/command construct cones RA_EU 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV norot dead 301 nofield
/sr1/command construct cones RA_MU 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV norot dead 302 nofield
/sr1/command construct cones RA_ER 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV rotRAnR dead 303 nofield
/sr1/command construct cones RA_MR 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV rotRAnR dead 304 nofield
/sr1/command construct cones RA_ED 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV rotRAnD dead 305 nofield
/sr1/command construct cones RA_MD 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV rotRAnD dead 306 nofield
/sr1/command construct cones RA_EL 45.1469 62.5 33.5 39.0 22.5 48.81 82.38 Steel 0 0 8.5 log_MCPV rotRAnL dead 307 nofield
/sr1/command construct cones RA_ML 56.7937 62.5 45.147 62.5 22.5 48.81 82.38 Steel 0 0 -36.5 log_MCPV rotRAnL dead 308 nofield
# Dummy, thin cylindres used for applying the SAME RA field-map (ROTATED by 90 deg.) to different anodes.
# NOTE: EM field cannot be applied to non simply connected bodies, as e.g. rings, cones, tori, etc.!
/sr1/command construct tubs RA_U 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.50 log_MCPV norot dead 322 nofield
/sr1/command construct tubs RA_R 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.52 log_MCPV rotRAnR dead 324 nofield
/sr1/command construct tubs RA_D 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.54 log_MCPV rotRAnD dead 326 nofield
/sr1/command construct tubs RA_L 0 0.01 0.005 0 360 G4_Galactic 0 0 -50.56 log_MCPV rotRAnL dead 328 nofield
# RA_G - Ring Anode Ground Cylinder
/sr1/command construct tubs RA_G 58 62.5 58.0 0 360 G4_Cu 0 0 -129.0 log_MCPV norot dead 351 nofield
# Ring Anodes visual attributes (optional)
/sr1/command visattributes log_RA_EU oxsteel
/sr1/command visattributes log_RA_MR oxsteel
/sr1/command visattributes log_RA_G Grid_style
# Alternative placement using World as a mother volume (mcpv_z = -92.5 mm). Check latter. These values refer to a 5 mm GAP!
#/sr1/command construct cones RA_EU 45.1469 62.5 33.5 39.0 22.5 48.1711 83.6578 Steel 0 0 -81 log_World norot dead 301 nofield
#/sr1/command construct cones RA_MU 56.7937 62.5 45.147 62.5 22.5 48.1711 83.6578 Steel 0 0 -126 log_World norot dead 302 nofield
#===============================================================================================================
# Gate Valve Area - Hosts the Ground Anode (upstream part of the Ring Anode) - Delimited by L3F1 and F200 flanges
#===============================================================================================================
# GATS - 200 CF flange upstream of MCP2 tube covering the whole length of the gate valve chamber. For simplicity, we
# choose the INNER diameter of the GATe valve Steel tube the same as the OUTER diameter of F200 (200 CF flange)
/sr1/command construct tubs GATS 103.25 126.5 92.5 0 360 Steel 0 0 -254.5 log_World norot dead 371 nofield
# Vacuum "Ring" (to avoid intersections with MCPV) - Not needed if world is already filled with vacuum.
#*/sr1/command construct tubs GATV 76.5 103.25 92.5 0 360 G4_Galactic 0 0 -254.5 log_World norot dead 370 nofield
# F200 - 200 CF flange upstream of MCP2 tube to connect to gate valve chamber
/sr1/command construct tubs F200 76.5 103.25 12 0 360 Steel 0 0 -174.0 log_World norot dead 372 nofield
# NOTE: When using GATV comment the F200 above and uncomment the following (change mother to log_GATV).
#*/sr1/command construct tubs F200 76.5 103.25 12 0 360 Steel 0 0 80.5 log_GATV norot dead 372 nofield
# Gate Valve Area visual attributes (optional)
/sr1/command visattributes log_GATS SCINT_style
/sr1/command visattributes log_F200 blue_style
#===============================================================================================================
# L3 - 3rd Einzel Lens # L3z = -56.7 cm. ATT: DUMMY FIELD change to electric L3Efield!
#===============================================================================================================
# Lens Gap = 12.0 mm => G/D = 12/130 ~ 0.1 (Lens Gap = gap between Ground and Anode, D - Diameter)
# L3 envelope (Tube + Flanges)
# L3VA - Lens 3 Vacuum
# L3ST - Lens 3 Steel tube (inner dia: 200 mm, outer dia: 206 mm, length: 720 mm)
# L3F1 - Lens 3 Flange 1, z = L3z + 208 mm
# L3F2 - Lens 3 Flange 2, z = L3z - 208 mm
/sr1/command construct tubs L3VA 0 100 220 0 360 G4_Galactic 0 0 -567 log_World norot dead 400 nofield
/sr1/command construct tubs L3ST 100 103 220 0 360 Steel 0 0 -567 log_World norot dead 401 nofield
/sr1/command construct tubs L3F1 103 126.5 12 0 360 Steel 0 0 -359 log_World norot dead 402 nofield
/sr1/command construct tubs L3F2 103 126.5 12 0 360 Steel 0 0 -775 log_World norot dead 403 nofield
# GPn - Ground Potential Electrodes
# n = 1-4 and 5-8 - components of the Ground Electrodes
# GP1 - Ground Electrode (inner dia: 130 mm, outer dia: 134 mm, length: 133 mm)
# GP2 - outer electrode surface (LN2 cooling vessel)
# GP3 - first ring cap
# GP4 - second ring cap
# n = 1-4 - Ground Electrode 1 (further from TD). See above for the meaning of acronyms.
/sr1/command construct tubs L3GP1 65 67 66.5 0 360 Steel 0 0 133.5 log_L3VA norot dead 421 nofield
/sr1/command construct tubs L3GP2 81 83 66.5 0 360 Steel 0 0 133.5 log_L3VA norot dead 422 nofield
/sr1/command construct tubs L3GP3 67 81 4 0 360 Steel 0 0 196.0 log_L3VA norot dead 423 nofield
/sr1/command construct tubs L3GP4 67 81 4 0 360 Steel 0 0 71.0 log_L3VA norot dead 424 nofield
# n = 5-8 - Ground Electrode 2 (closer to TD). See above for the meaning of acronyms.
/sr1/command construct tubs L3GP5 65 67 66.5 0 360 Steel 0 0 -133.5 log_L3VA norot dead 431 nofield
/sr1/command construct tubs L3GP6 81 83 66.5 0 360 Steel 0 0 -133.5 log_L3VA norot dead 432 nofield
/sr1/command construct tubs L3GP7 67 81 4 0 360 Steel 0 0 -196.0 log_L3VA norot dead 433 nofield
/sr1/command construct tubs L3GP8 67 81 4 0 360 Steel 0 0 -71.0 log_L3VA norot dead 434 nofield
# HP - High Potential Electrode (Central Anode - usually at +8.7 kV, for a 15 keV muon beam)
/sr1/command construct tubs L3HP 65 83 55 0 360 Steel 0 0 0 log_L3VA norot dead 451 nofield
# Lens 3 visual attributes (optional)
/sr1/command visattributes log_L3VA invisible
/sr1/command visattributes log_L3ST invisible
/sr1/command visattributes log_L3HP darkred
#===============================================================================================================
# IP - Intermediate Piece (between Trigger Detector and Einzel Lens 3)
# Original name was CGate - B-field Coil Compensation Gate?! # CompGatez = -86.55 cm; // L3z - 22 - 7.85 cm
#===============================================================================================================
# IPV (but also others) are just empty volumes "filled" with vacuum. Sometimes used to apply EM fields to restricted areas.
/sr1/command construct tubs IPV 0 100 78.5 0 360 G4_Galactic 0 0 -865.5 log_World norot dead 500 nofield
# IPCF - Intermediate Piece Central Flange (same as L3 Flanges)
# IPST - Intermediate Piece Steel Tube (same diameter as L3 Steel Tube)
/sr1/command construct tubs IPCF 103 126.5 12.0 0 360 Steel 0 0 -865.5 log_World norot dead 501 nofield
/sr1/command construct tubs IPST 100 103 78.5 0 360 Steel 0 0 -865.5 log_World norot dead 502 nofield
# IP visual attributes (optional)
/sr1/command visattributes log_IPV invisible
/sr1/command visattributes log_IPCF blue_style
/sr1/command visattributes log_IPST gray
#===============================================================================================================
# Trigger - Trigger Detector # Triggerz = -1092 mm
#===============================================================================================================
# Trigger tube and relative vacuum
/sr1/command construct tubs TriggerV 0 100 148 0 360 G4_Galactic 0 0 -1092 log_World norot dead 600 nofield
/sr1/command construct tubs Trigger 100 103 148 0 360 Steel 0 0 -1092 log_World norot dead 601 nofield
# TF - Trigger tube flanges
/sr1/command construct tubs TF1 103 126.5 12 0 360 Steel 0 0 -956 log_World norot dead 611 nofield
/sr1/command construct tubs TF2 103 126.5 12 0 360 Steel 0 0 -1228 log_World norot dead 612 nofield
# Carbon Foil (default HALF-thickness 0.000005147 mm, see below => CFoil thick = 10.3 nm).
# USE THE NAME CFoil or coulombCFoil, otherwise sr1MuFormation won't work!
####/sr1/command construct box CFoil 60 60 0.000005147 G4_GRAPHITE 0 0 -45 log_TriggerV norot dead 621 nofield
/sr1/command construct box CFoil 60 60 0.000006471 G4_GRAPHITE 0 0 -45 log_TriggerV norot dead 621 nofield
####/sr1/command construct box coulombCFoil 60 60 0.000005147 G4_GRAPHITE 0 0 -45 log_TriggerV norot dead 621 nofield
# Notes: NIST tables use G4_GRAPHITE with 1.7 g/cm3 and 78 eV ioniz. energy.
# An area density of 2.20 ug/cm2 implies a CF thickn. = (2.20*1.e-6/1.70)*cm = 1.294e-5 mm - Total thickness
# An area density of 1.75 ug/cm2 implies a CF thickn. = (1.75*1.e-6/1.70)*cm = 1.029e-5 mm - Total thickness
# If necessary, use Graphite as defined in sr1DetectorConstruction.cc and set any density.
# Dummy plane to intercept outgoing muons from the Carbon foil.
/sr1/command construct box saveCFoil 60 60 5e-4 G4_Galactic 0 0 -44.995 log_TriggerV norot dead 623 nofield
# Electrical Field areas in the Trigger Detector
# En = Electrical Field n: TnFieldMgr (n = 1-3)
# Original TriggE2: [4.*sqrt(2), 4.5, 0.7/sqrt(2)] cm -> changed due to overlaps with E1 and E3
/sr1/command construct box TriggE1 45 45 4 G4_Galactic 0 0 -38.5 log_TriggerV norot dead 631 nofield
/sr1/command construct box TriggE2 45 45 4.9497 G4_Galactic 0 0 2.25 log_TriggerV rotTrig dead 632 nofield
/sr1/command construct box TriggE3 45 45 4 G4_Galactic 0 0 43.0 log_TriggerV norot dead 633
# Beam spot (just for having a visual idea!)
/sr1/command construct tubs BSpot 0 20 1 0 360 G4_Galactic 0 0 -48.0 log_TriggerV norot dead 650 nofield
# Trigger visual attributes (optional)
/sr1/command visattributes log_TriggerV invisible
/sr1/command visattributes log_Trigger invisible
#/sr1/command visattributes log_CFoil MACOR_style
/sr1/command visattributes log_BSpot fblue_style
# One can set visible attrib. also on a MATERIAL basis, rather than on log_VOL.
# E.g. /sr1/command visattributes Steel red
################################################################################################################
# -- Setting the ELECTRIC and MAGNETIC fields --
################################################################################################################
# Use ABSOLUTE coordinates to specify the field position (i.e. with respect to GLOBAL WORLD)!
# Default field units: Magnetic - T, Electric - kV/mm (or kV for E-field maps).
# NOTE: Applying a field to an invisible log_vol makes is visible!
### Electric field at TRIGGER Detector TD: Three different uniform fields
/sr1/command globalfield Trigg1_field 0. 0. -1130.5 uniform log_TriggE1 0 0 0 0 0 -0.02375
/sr1/command globalfield Trigg2_field 0. 0. -1089.75 uniform log_TriggE2 0 0 0 0 0 0.041416
/sr1/command globalfield Trigg3_field 0. 0. -1049.0 uniform log_TriggE3 0 0 0 0 0 -0.49375
### Electric field at Einzel LENS 3 - from folded 2D axial field map (f - folded, i.e. symmetric)
# Typically V = +8.7 kV for a muon beam at 15 keV. Use either 2DE L3_Erz.map or
#
# ATTENTION: The electric field is ANTI-symmetric: DO NOT use folded field map L3_Erz4.map!!
#
/sr1/command globalfield Lens3_field 0. 0. -567. fromfile 2DE L3_Erz.map log_L3VA 6.78
# To change the field in regular steps (e.g. for testing) use (f_min f_max step_no), e.g.:
#*/sr1/command globalfield Lens3_field 0. 0. -567. fromfile 2DE L3_Erz.map log_L3VA 7 11 5
### Electric field at RING ANODE - from 3DE field map
# Typically set at +11.0 kV for a muon beam at 15 keV
# To create an arbitrary configuration, switch on all fields and set different potentials.
#######/sr1/command globalfield RngAnU_field 0. 0. -167.00 fromfile 3DE EM_3D_extc.map log_RA_U 11.0
/sr1/command globalfield RngAnU_field 0. 0. -143.00 fromfile 3DE EM_3D_ext_gridf.map log_RA_U 8.6
/sr1/command globalfield RngAnR_field 0. 0. -143.02 fromfile 3DE EM_3D_ext_gridf.map log_RA_R 8.6
/sr1/command globalfield RngAnD_field 0. 0. -143.04 fromfile 3DE EM_3D_ext_gridf.map log_RA_D 8.6
/sr1/command globalfield RngAnL_field 0. 0. -143.06 fromfile 3DE EM_3D_ext_gridf.map log_RA_L 8.6
### LAST FIELD OK: EM_3D_ext2f.map
# EXTENDED MAPS (from -28 to +20 mm) EM_3D_extc.map or EM_3D_extf.map (coord + field or field only)
# "Best" results with EM_RA_3D.map, even though this is not a precise map (poor meshing).
#RA_1kV_upf_raw.map
#RA_1kV_upf.map # EM_extendf2.map give rise to "strange spots" -> sr1_1047_RA13.eps
### Electric field at SAMPLE space. Three possible field maps: Sample, G1 and G2 (closest to sample)
# To create an arbitrary configuration, switch on all fields and set different potentials.
# Field extension along z: 50 mm. Center of MCPV at z = -92.5 mm. Sample face at: z = 108.5 (4 mm thick - SAH2).
# Position of the field: (-92.5 - 50/2) + (108.5 - 4/2) = -11 mm wrt. WORLD! =>
#*/sr1/command globalfield Sample_field 0. 0. -11.0 fromfile 2DE sample_Erz.map log_MCPV 9.0
#*/sr1/command globalfield Guard2_field 0. 0. -11.0 fromfile 2DE guard2_Erz.map log_MCPV 6.0
#*/sr1/command globalfield Guard1_field 0. 0. -11.0 fromfile 2DE guard1_Erz.map log_MCPV 3.0
### Magnetic field at SAMPLE space (use either a uniform field or a field map).
# Use either DMCP or MCPV to apply a CONSTANT field strictly to the sample or also to the surroundings, resp.!
#*/sr1/command globalfield Magnet_field 0. 0. 14.5 uniform log_DMCP 0 0.005 0 0 0 0
# Use field map to set field to 20 G TF
# Extended map has -100/100 cm z extension, field center at +70, original -85/85, center at +85
#######/sr1/command globalfield Magnet_field 0. 0. -836.0 fromfile 2DB sample_Brz.map log_IPV 0.002
#*/sr1/command globalfield Magnet_field 0. 0. -686.0 fromfile 2DB sample_Brz_ext.map log_L3VA 0.002
# Set parameters for particle tracking in an EM field
/sr1/command globalfield setparameter SetLargestAcceptableStep 5
/sr1/command globalfield setparameter SetMinimumEpsilonStep 5e-5
/sr1/command globalfield setparameter SetMaximumEpsilonStep 0.001
/sr1/command globalfield setparameter SetDeltaOneStep 0.1
/sr1/command globalfield setparameter SetDeltaIntersection 0.01
/sr1/command globalfield printparameters
################################################################################################################
# -- Testing the ELECTRIC and MAGNETIC fields (OPTIONAL) --
################################################################################################################
# FIELD CHECKS at different beam transport components (from trigg. to sample)
# All distances in mm and in GLOBAL coordinates (preferably at field center)
# The test points below are just some examples. Any point can be checked.
# Trigger 1
/sr1/command globalfield printFieldValueAtPoint 0. 0. -1130.5
# Trigger 2
/sr1/command globalfield printFieldValueAtPoint 0. 0. -1089.75
# Trigger 3
/sr1/command globalfield printFieldValueAtPoint 0. 0. -1049.0
# Einzel Lens 3 - L3 (center at -567.0, but max field at rel. +/-62 mm)
/sr1/command globalfield printFieldValueAtPoint 0. 0. -507.0
# Ring Anode - RA (center at -167.0, but max field at rel. -16/+132 mm)
/sr1/command globalfield printFieldValueAtPoint 0. 0. -35.0
# Sample space (center at -11.0, but max field at rel. +24 mm)
/sr1/command globalfield printFieldValueAtPoint 0. 0. 13.0
# Check magnetic field at sample space
#*/sr1/command globalfield printFieldValueAtPoint 10. 0. 13.0
# Check magnetic field at the lower field end limit
#*/sr1/command globalfield printFieldValueAtPoint 0. 20. -1680.
################################################################################################################
# -- Setting simulation PARAMETERS --
################################################################################################################
# Set processes from: lowenergy, penelope, coulombAndMultiple (default, Coul. only for CFoil), coulomb (for all, very slow).
/sr1/command typeofprocesses coulombAndMultiple
#*/sr1/command typeofprocesses penelope
#*/sr1/command includeMuoniumProcesses false
# Set the overall range cut (default 0.1 mm)
#*/run/setCut 1 mm
# Set the range cut on particular volumes (in mm)
/sr1/command SetUserLimits log_CFoil 1e-8
#*/sr1/command SetUserLimits log_target 0.01
#*/sr1/command SetUserLimits log_targetscint 0.01
#*/sr1/command SetUserLimits log_cryostatscint 0.01
# Set particle energy cuts on particular volumes (in eV)
/sr1/command SetUserMinEkine log_World 0.1
# Store ALL the events in a ROOT tree or just the interesting ones? (default is true)
#*/sr1/command storeOnlyEventsWithHits false
# Set the minimum time separation between two subsequent signals in the same detector (in ns)
/sr1/command signalSeparationTime 0.1
# Override runID number
#*/sr1/run/runID 21
# Set the frequency of event printing
/sr1/run/howOftenToPrintEvent 1000
# RANDOM option choices: (specify the random number generator initialisation)
# 0 ... no initialisation (default)
# 1 ... use actual computer time to initialise now
# 2 ... use event number to initialise at the beginning of each event
# 3 ... read in the random no. initial values for each event from a file
/sr1/run/randomOption 2
# VISUALIZATION options
# To enable or disable visualization uncomment one of these lines
# To modify visualization options edit the file vis.mac
/vis/disable
#*/control/execute vis.mac
################################################################################################################
# -- Setting PARTICLE GUN parameters --
################################################################################################################
# Default momentum direction: 001, i.e. 0z.
# Default muon spin direction: 100, i.e. 0x.
# Default particle type: mu+ (can be changed to Mu)
# Set particle type
#*/gun/particle Mu
/gun/particle mu+
# Set beam vertex (CFoil at -1137. To start after C Foil use z = -1130, before -1170)
/gun/vertex 0. 0. -1138. mm
# A point-like uniform beam
/gun/vertexsigma -0.1 -0.1 0 mm
# Set beam transverse spread (default GAUSSIAN spread)
# If FWHM = 10 mm ==> sigma = 10/2.354 = 4.2481 mm (last 0 is a dummy value)
# Negative sigma values => random FLAT RECTANGULAR distribution (area 2x.2y)
# Use vertexboundary with (vb < sigma_xy) to obtain a CIRCULAR beam spot
# /gun/vertexsigma 0 0 0 mm ==> Very SLOW with mag. field ON and centered beam
#*/gun/vertexsigma 42.5 42.5 0 mm
/gun/vertexsigma -20 -20 0 mm
/gun/vertexboundary 20 -1e6 1e6 mm
# /gun/vertexboundary: rMaxAllowed, zMinAllowed, zMaxAllowed # Beam AND gating
#*/gun/vertexboundary 7 -1314.4 -1305 mm
# Without restrictions in z, but only on r:
#*/gun/vertexboundary 3 -1e6 1e6 mm
# Set beam momentum (USE only as an ALTERNATIVE to setting energy!)
# /gun/momentum 0 0 29.79 MeV
#*/gun/momentum 0 0 1.8 MeV
# Energy loss at p = 1.2 MeV/c (E = 6.8 keV) => 1.23 +/- 0.2 keV
# Energy loss at p = 1.8 MeV/c (E = 15.3 keV) => 1.25 +/- 0.3 keV
# 1.2 MeV/c -> 6.8 keV, 1.8 MeV/c -> 15.3 keV
# muon rest mass = 105.658 MeV/c2
# With Trigger ON use 15 keV, with Trigger OFF use 11.27 keV (15 - 3.73 = 11.27 keV)
/gun/kenergy 12.0 keV
# Set beam momentum direction
/gun/direction 0.0 0.0 1.0
# Set muon spin direction
/gun/muonPolarizVector 1 0 0
# Other useful test parameters:
#
# FWHM = 3% ==> sigma = 29.79*0.03/2.354 = 0.37965 MeV/c
#*/gun/momentumsmearing 0.37965 MeV
#---/gun/momentumboundary: pMinAllowed, pMaxAllowed, dummy
#*/gun/momentumboundary 20 40 0 MeV
#---/gun/tilt: xangle, yangle, dummy
#*/gun/tilt 0 0.5 0 deg
#---/gun/tiltsigma: xangleSigma, yangleSigma, dummy (1 degree at 1 m => 17 mm)
#*/gun/tiltsigma 0.2 0.2 0 deg
#*/gun/pitch 0.5 deg
#---/gun/decaytimelimits: decayMin, decayMax, decayTime
#*/gun/decaytimelimits 10400 10420 2197.03 ns
# Selectively inactivate or activate sensitive detectors
#*/hits/inactivate /sr1/ScintSD
# Only for code debugging!
#/tracking/verbose 1
# BEAM ON
#*/run/beamOn 1000000
#*/run/beamOn 2
/run/beamOn 10000

View File

@ -0,0 +1,289 @@
# Macro file for sr1.cc - Construct detector, set fields and other parameters.
# Last modified by T. Shiroka: 31.10.2008
# How to run: sr1 10xx.mac (append "idle" for prompt after running)
# sr1 10xx.mac > fname.txt (stores output on a txt file)
###############################################################################################
# #
# Specify the geometry parameters in this file (all dimensions in mm) #
# a. Lines starting with hash marks "#" are comments #
# b Lines starting with #* are temporary comments. Remove/modify to change the configuration #
# c. Lines starting with /sr1/command are commands for the executable program #
# d. Lines starting with /vis, /gun, etc. are common macro commands #
# e. Beam-line components are ordered from exp. area (MCP2) to trigger detector (TD) #
#---------------------------------------------------------------------------------------------#
# Syntax example (following /sr1/command): #
# construct solid_type volume_name parameters_defining_solid material position mothers_name #
# (mothers_name starts with log_) #
###############################################################################################
# For the meaning of the acronyms see also the original G3 file ugeom.F at:
# http://savannah.psi.ch/viewcvs/trunk/simulation/geant3/src/lemsr/ugeom.F?root=nemu%2Flem&rev=2964&view=markup
################################################################################################################
# -- ROTATION MATRICES --
################################################################################################################
# 3 parameters -> Define Euler angles (the 4th par. is set to zero).
# 4 parameters -> Define axis + rotation.
# HEP computations ordinarily use the active rotation viewpoint (object is rotated NOT axes).
# Therefore, rotations about an axis imply ACTIVE COUNTER-CLOCKWISE rotation in this package.
# Rotation around a specified axis means counter-clockwise rot. around the positive direction of the axis.
# Define rotations for the field maps of Trigger and Ring Anode:
/sr1/command rotation rotPole 1 0 0 90
/sr1/command rotation rotTrig 0 1 0 -90
/sr1/command rotation rotRAnR 0 0 1 -90
/sr1/command rotation rotRAnL 0 0 1 90
/sr1/command rotation rotRAnD 0 0 1 180
################################################################################################################
# -- SPIN ROTATOR GEOMETRY --
################################################################################################################
# WORLD = Laboratory reference frame
/sr1/command construct box World 250 250 500 G4_Galactic 0 0 0 no_logical_volume norot dead -1
# World visual attributes (optional)
/sr1/command visattributes log_World invisible
# Electromagnet with iron yoke
# U - up, D - down, L - left, R - right (with respect to muon's view - momentum direction)
/sr1/command construct box YokeU 200 25 50 G4_Fe 0 200 0 log_World norot dead 201 nofield
/sr1/command construct box YokeD 200 25 50 G4_Fe 0 -200 0 log_World norot dead 204 nofield
/sr1/command construct box YokeL 25 175 50 G4_Fe 175 0 0 log_World norot dead 203 nofield
/sr1/command construct box YokeR 25 175 50 G4_Fe -175 0 0 log_World norot dead 202 nofield
/sr1/command construct tubs NPole 0 100 25 0 360 G4_Fe 0 100 0 log_World rotPole dead 211 nofield
/sr1/command construct tubs NYoke 0 50 25 0 360 G4_Fe 0 150 0 log_World rotPole dead 212 nofield
/sr1/command construct tubs NCoil 50 60 25 0 360 G4_Cu 0 150 0 log_World rotPole dead 213 nofield
/sr1/command construct tubs SPole 0 100 25 0 360 G4_Fe 0 -100 0 log_World rotPole dead 221 nofield
/sr1/command construct tubs SYoke 0 50 25 0 360 G4_Fe 0 -150 0 log_World rotPole dead 222 nofield
/sr1/command construct tubs SCoil 50 60 25 0 360 G4_Cu 0 -150 0 log_World rotPole dead 223 nofield
# Capacitor
/sr1/command construct box Cap_p 0.5 50 150 Brass 50 0 0 log_World norot dead 301 nofield
/sr1/command construct box Cap_n 0.5 50 150 Brass -50 0 0 log_World norot dead 302 nofield
#*/sr1/command construct box Uniform 49 49 150 G4_Galactic 0 0 0 log_World norot dead 303 nofield
# Visualize Electric and Magnetic field volumes. COMMENT during normal run, otherwise the root file is empty!
#/sr1/command construct box Evol 50 50 300 G4_Galactic 0 0 0 log_World norot dead 311 nofield
#/sr1/command construct box Bvol 70 70 500 G4_Galactic 0 0 0 log_World norot dead 312 nofield
# Dummy, thin cylindres used for applying the field maps.
# NOTE: EM field cannot be applied to non simply connected bodies, as e.g. rings, cones, tori, etc.!
/sr1/command construct tubs EField 0 0.01 0.005 0 360 G4_Galactic 0 0 0.0 log_World norot dead 322 nofield
/sr1/command construct tubs BField 0 0.01 0.005 0 360 G4_Galactic 0 0 0.1 log_World norot dead 324 nofield
# Beam spot (just for having a visual idea!). Assume r = 3*sigma to intercept 99.73% of all events, with sigma = 6.83 mm
#/sr1/command construct tubs BSpot 0 30 1 0 360 G4_Galactic 0 0 -490 log_World norot dead 650 nofield
/sr1/command construct tubs saveBSpot 0 20.5 1 0 360 G4_Galactic 0 0 -490 log_World norot sr1/ScintSD 650 nofield
# Dummy plane or target to intercept muon events.
# Use it either as (DMCP-sr1/ScintSD) - with no info on mu+ polariz., or as (target-dead) with info on mu+ polariz.
#/sr1/command construct box saveCFoil 65 60 5e-4 G4_Galactic 0 0 -199 log_World norot dead 623 nofield
/sr1/command construct tubs target 0 65 1.5 0 360 G4_Galactic 0 0 490 log_World norot dead 101 nofield
#/sr1/command construct tubs DMCP 0 76.5 1.5 0 360 MCPglass 0 0 108 log_MCPV norot dead 202 nofield
# Visual attributes (optional)
/sr1/command visattributes log_NPole red
/sr1/command visattributes log_SPole lightblue
/sr1/command visattributes log_NCoil Grid_style
/sr1/command visattributes log_SCoil Grid_style
/sr1/command visattributes log_YokeR oxsteel
/sr1/command visattributes log_Cap_p darkred
/sr1/command visattributes log_Cap_n blue_style
#/sr1/command visattributes log_Uniform invisible
################################################################################################################
# -- Setting the ELECTRIC and MAGNETIC fields --
################################################################################################################
# Use ABSOLUTE coordinates to specify the field position (i.e. with respect to GLOBAL WORLD)!
# Default field units: Magnetic - T, Electric - kV/mm (or kV for E-field maps).
# NOTE: Applying a field to an invisible log_vol makes is visible!
### Uniform Electric and Magnetic field at Spin Rotator:
#*/sr1/command globalfield UniMag_field 0. 0. 0. uniform log_Uniform 0 -0.036 0 0 0 0
#*/sr1/command globalfield UniEle_field 0. 0. 0. uniform log_Uniform 0 0 0 -0.21 0 0
#*/sr1/command globalfield UniEM_field 0. 0. 0. uniform log_Uniform 0 -0.036 0 -0.21 0 0
###15 cm gap /sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_capacitor.map log_EField 22.1
###1000-0 V /sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_capacitor.map log_EField 14.74
#*/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_const.map log_EField 21.0
## Symmetric field map +/- 500 V.
/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_cap_symm.map log_EField 21.0
/sr1/command globalfield Mag_field 0. 0. 0.1 fromfile 3DB B3D_iron2.map log_BField 0.036
# Set parameters for particle tracking in an EM field
/sr1/command globalfield setparameter SetLargestAcceptableStep 5
/sr1/command globalfield setparameter SetMinimumEpsilonStep 5e-5
/sr1/command globalfield setparameter SetMaximumEpsilonStep 0.001
/sr1/command globalfield setparameter SetDeltaOneStep 0.1
/sr1/command globalfield setparameter SetDeltaIntersection 0.01
/sr1/command globalfield printparameters
# TESTING EM FIELD
/sr1/command globalfield printFieldValueAtPoint 0. 0. 0.
################################################################################################################
# -- Setting simulation PARAMETERS --
################################################################################################################
# Set processes from: lowenergy, penelope, coulombAndMultiple (default, Coul. only for CFoil), coulomb (for all, very slow).
#**/sr1/command typeofprocesses coulombAndMultiple
#*/sr1/command typeofprocesses penelope
#*/sr1/command includeMuoniumProcesses false
# Set the overall range cut (default 0.1 mm)
#*/run/setCut 1 mm
# Set the range cut on particular volumes (in mm)
#*/sr1/command SetUserLimits log_target 0.01
#*/sr1/command SetUserLimits log_targetscint 0.01
#*/sr1/command SetUserLimits log_cryostatscint 0.01
# Set particle energy cuts on particular volumes (in eV)
#*/sr1/command SetUserMinEkine log_World 0.1
# Store ALL the events in a ROOT tree or just the interesting ones? (default is true)
#*/sr1/command storeOnlyEventsWithHits false
# Set the minimum time separation between two subsequent signals in the same detector (in ns)
/sr1/command signalSeparationTime 0.1
# Override runID number
#*/sr1/run/runID 21
# Set the frequency of event printing
/sr1/run/howOftenToPrintEvent 1000
# RANDOM option choices: (specify the random number generator initialisation)
# 0 ... no initialisation (default)
# 1 ... use actual computer time to initialise now # Pseudo-random numbers
# 2 ... use event number to initialise at the beginning of each event # Reproducible numbers
# 3 ... read in the random no. initial values for each event from a file
/sr1/run/randomOption 1
# VISUALIZATION options
# To enable or disable visualization uncomment one of these lines
# To modify visualization options edit the file vis.mac
/vis/disable
#*/control/execute vis.mac
#alias g4='export G4WORKDIR="/home/sedlak/bin_4.9.1"; source /home/geant4/4.9.1/env.sh;
#export G4VRMLFILE_VIEWER="vrmlview"; echo "On this machine the G4VRMLFILE_VIEWER=$G4VRMLFILE_VIEWER"'
################################################################################################################
# -- Setting PARTICLE GUN parameters --
################################################################################################################
# Default momentum direction: 001, i.e. 0z.
# Default muon spin direction: 100, i.e. 0x.
# Default particle type: mu+ (can be changed to Mu or proton)
# Set particle type
#*/gun/particle Mu
#*/gun/particle proton
/gun/particle mu+
# Set the position of the beam vertex
/gun/vertex 0. 0. -495. mm
# A point-like uniform beam
#*/gun/vertexsigma -0.1 -0.1 0 mm
# Set beam transverse spread (default GAUSSIAN spread)
# If FWHM = 10 mm ==> sigma = 10/2.354 = 4.2481 mm (last 0 is a dummy value)
# Negative sigma values => random FLAT RECTANGULAR distribution (area 2x.2y)
# Use vertexboundary with (vb < sigma_xy) to obtain a CIRCULAR beam spot
# /gun/vertexsigma 0 0 0 mm ==> Very SLOW with mag. field ON and centered beam
## The LEMUSR beam has a sigma of 6.83 mm => FWHM = 2.355*6.83 = 16.08 mm
/gun/vertexsigma 6.83 6.83 0 mm
#/gun/vertexsigma -20 -20 0 mm
#/gun/vertexboundary 20 -1e6 1e6 mm
# /gun/vertexboundary: rMaxAllowed, zMinAllowed, zMaxAllowed # Beam AND gating
#*/gun/vertexboundary 7 -1314.4 -1305 mm
# Without restrictions in z, but only on r:
#*/gun/vertexboundary 3 -1e6 1e6 mm
# Set beam momentum (USE only as an ALTERNATIVE to setting energy!)
# /gun/momentum 0 0 29.79 MeV
#*/gun/momentum 0 0 1.8 MeV
# Energy loss at p = 1.2 MeV/c (E = 6.8 keV) => 1.23 +/- 0.2 keV
# Energy loss at p = 1.8 MeV/c (E = 15.3 keV) => 1.25 +/- 0.3 keV
# 1.2 MeV/c -> 6.8 keV, 1.8 MeV/c -> 15.3 keV, 2.056 MeV/c -> 20 keV.
# muon rest mass = 105.658 MeV/c2
# With Trigger ON use 15 keV, with Trigger OFF use 11.27 keV (15 - 3.73 = 11.27 keV)
#/gun/kenergy 3.4 MeV
/gun/kenergy 16.27 keV
# Set beam momentum direction
/gun/direction 0.0 0.0 1.0
# Set muon spin direction
/gun/muonPolarizVector 1 0 0
# Other useful test parameters:
#
# FWHM = 3% ==> sigma = 29.79*0.03/2.354 = 0.37965 MeV/c
#*/gun/momentumsmearing 0.37965 MeV
#---/gun/momentumboundary: pMinAllowed, pMaxAllowed, dummy
#*/gun/momentumboundary 20 40 0 MeV
#---/gun/tilt: xangle, yangle, dummy
#*/gun/tilt 0 0.5 0 deg
#---/gun/tiltsigma: xangleSigma, yangleSigma, dummy (1 degree at 1 m => 17 mm)
#*/gun/tiltsigma 0.2 0.2 0 deg
#*/gun/pitch 0.5 deg
#---/gun/decaytimelimits: decayMin, decayMax, decayTime
#*/gun/decaytimelimits 10400 10420 2197.03 ns
# Selectively inactivate or activate sensitive detectors
#*/hits/inactivate /sr1/ScintSD
# Only for code debugging!
#/tracking/verbose 1
# BEAM ON
#*/run/beamOn 1000000
#/run/beamOn 2
#/run/beamOn 1000000
/run/beamOn 5000

View File

@ -0,0 +1,367 @@
# Macro file for sr1.cc - Construct detector, set fields and other parameters.
# Last modified by T. Shiroka: 31.10.2008
# How to run: sr1 10xx.mac (append "idle" for prompt after running)
# sr1 10xx.mac > fname.txt (stores output on a txt file)
###############################################################################################
# #
# Specify the geometry parameters in this file (all dimensions in mm) #
# a. Lines starting with hash marks "#" are comments #
# b Lines starting with #* are temporary comments. Remove/modify to change the configuration #
# c. Lines starting with /sr1/command are commands for the executable program #
# d. Lines starting with /vis, /gun, etc. are common macro commands #
# e. Beam-line components are ordered from exp. area (MCP2) to trigger detector (TD) #
#---------------------------------------------------------------------------------------------#
# Syntax example (following /sr1/command): #
# construct solid_type volume_name parameters_defining_solid material position mothers_name #
# (mothers_name starts with log_) #
###############################################################################################
# For the meaning of the acronyms see also the original G3 file ugeom.F at:
# http://savannah.psi.ch/viewcvs/trunk/simulation/geant3/src/lemsr/ugeom.F?root=nemu%2Flem&rev=2964&view=markup
################################################################################################################
# -- ROTATION MATRICES --
################################################################################################################
# 3 parameters -> Define Euler angles (the 4th par. is set to zero).
# 4 parameters -> Define axis + rotation.
# HEP computations ordinarily use the active rotation viewpoint (object is rotated NOT axes).
# Therefore, rotations about an axis imply ACTIVE COUNTER-CLOCKWISE rotation in this package.
# Rotation around a specified axis means counter-clockwise rot. around the positive direction of the axis.
# Define rotations for the field maps of Trigger and Ring Anode:
/sr1/command rotation rotPole 1 0 0 90
/sr1/command rotation rotTrig 0 1 0 -90
/sr1/command rotation rotRAnR 0 0 1 -90
/sr1/command rotation rotRAnL 0 0 1 90
/sr1/command rotation rotRAnD 0 0 1 180
################################################################################################################
# -- SPIN ROTATOR GEOMETRY --
################################################################################################################
# WORLD = Laboratory reference frame
/sr1/command construct box World 350 250 900 G4_Galactic 0 0 0 no_logical_volume norot dead -1
# World visual attributes (optional)
/sr1/command visattributes log_World invisible
# Electromagnet with iron yoke
# U - up, D - down, L - left, R - right (with respect to muon's view - momentum direction)
/sr1/command construct box NPole 200 20 150 G4_Fe 0 95 0 log_World norot dead 201 nofield
/sr1/command construct box NYoke 150 30 50 G4_Fe 0 145 0 log_World norot dead 202 nofield
/sr1/command construct box NConn 75 17.5 50 G4_Fe 225 157.5 0 log_World norot dead 203 nofield
/sr1/command construct box YokeL 25 175 50 G4_Fe 325 0 0 log_World norot dead 210 nofield
/sr1/command construct box SPole 200 20 150 G4_Fe 0 -95 0 log_World norot dead 205 nofield
/sr1/command construct box SYoke 150 30 50 G4_Fe 0 -145 0 log_World norot dead 206 nofield
/sr1/command construct box SConn 75 17.5 50 G4_Fe 225 -157.5 0 log_World norot dead 207 nofield
# Front shield
/sr1/command construct box FShU 225 50 5 G4_Fe 0 125 -175 log_World norot dead 301 nofield
/sr1/command construct box FShD 225 50 5 G4_Fe 0 -125 -175 log_World norot dead 302 nofield
/sr1/command construct box FShL 75 75 5 G4_Fe 150 0 -175 log_World norot dead 303 nofield
/sr1/command construct box FShR 75 75 5 G4_Fe -150 0 -175 log_World norot dead 304 nofield
# Back shield
/sr1/command construct box BShU 225 50 5 G4_Fe 0 125 175 log_World norot dead 321 nofield
/sr1/command construct box BShD 225 50 5 G4_Fe 0 -125 175 log_World norot dead 322 nofield
/sr1/command construct box BShL 75 75 5 G4_Fe 150 0 175 log_World norot dead 323 nofield
/sr1/command construct box BShR 75 75 5 G4_Fe -150 0 175 log_World norot dead 324 nofield
# Top coils
/sr1/command construct box TCoil1 170 12.5 10 G4_Cu 0 127.5 -60 log_World norot dead 401 nofield
/sr1/command construct box TCoil2 170 12.5 10 G4_Cu 0 127.5 60 log_World norot dead 402 nofield
/sr1/command construct box TCoil3 10 12.5 50 G4_Cu 160 127.5 0 log_World norot dead 403 nofield
/sr1/command construct box TCoil4 10 12.5 50 G4_Cu -160 127.5 0 log_World norot dead 404 nofield
# Bottom coils
/sr1/command construct box BCoil1 170 12.5 10 G4_Cu 0 -127.5 -60 log_World norot dead 421 nofield
/sr1/command construct box BCoil2 170 12.5 10 G4_Cu 0 -127.5 60 log_World norot dead 422 nofield
/sr1/command construct box BCoil3 10 12.5 50 G4_Cu 160 -127.5 0 log_World norot dead 423 nofield
/sr1/command construct box BCoil4 10 12.5 50 G4_Cu -160 -127.5 0 log_World norot dead 424 nofield
# Capacitor
/sr1/command construct box Cap_p 0.5 50 150 Brass 50 0 0 log_World norot dead 501 nofield
/sr1/command construct box Cap_n 0.5 50 150 Brass -50 0 0 log_World norot dead 502 nofield
#*/sr1/command construct box Uniform 49 49 150 G4_Galactic 0 0 0 log_World norot dead 503 nofield
# Lens 1 - 1st Einzel Lens # L1z = -(38 + 25) = -63 cm.
# Gap = 10.0 mm => G/D = 10/80 = 0.125 (Lens Gap = gap between Ground and Anode, D - Diameter)
# L1ENV - Lens 1 envelope - for easy positioning of lens parts (outer dia: 100 mm, length: 300 mm)
/sr1/command construct tubs L1ENV 0 50 150 0 360 G4_Galactic 0 0 -630 log_World norot dead 600 nofield
# GPn - Ground Potential Electrodes. (n = 1-2, inner dia: 80 mm, outer dia: 84 mm, length: 100 mm)
# n = 1 - Ground Electrode 1 (further from SR).
/sr1/command construct tubs L1GP1 40 44 50 0 360 Steel 0 0 -100 log_L1ENV norot dead 651 nofield
# n = 2 - Ground Electrode 2 (closer to SR).
/sr1/command construct tubs L1GP2 40 44 50 0 360 Steel 0 0 100 log_L1ENV norot dead 653 nofield
# HP - High Potential Electrode (Central Anode - usually at +8.5 kV, for a 15 keV muon beam)
/sr1/command construct tubs L1HP 40 44 40 0 360 Steel 0 0 0 log_L1ENV norot dead 652 nofield
# Lens 1 visual attributes (optional)
/sr1/command visattributes log_L1ENV invisible
/sr1/command visattributes log_L1HP darkred
# Dummy, thin cylindres used for applying the field maps.
# NOTE: EM field cannot be applied to non simply connected bodies, as e.g. rings, cones, tori, etc.!
/sr1/command construct tubs EField 0 0.01 0.005 0 360 G4_Galactic 0 0 0.0 log_World norot dead 702 nofield
/sr1/command construct tubs BField 0 0.01 0.005 0 360 G4_Galactic 0 0 0.1 log_World norot dead 704 nofield
/sr1/command construct tubs LField 0 0.01 0.005 0 360 G4_Galactic 0 0 0.0 log_L1ENV norot dead 706 nofield
# Visualize Electric and Magnetic field volumes. COMMENT during normal run, otherwise the root file is empty!
#*/sr1/command construct box Evol 50 50 300 G4_Galactic 0 0 0 log_World norot dead 752 nofield
#*/sr1/command construct box Bvol 70 70 500 G4_Galactic 0 0 0 log_World norot dead 754 nofield
#*/sr1/command construct tubs Lvol 0 38 180 0 360 G4_Galactic 0 0 -630 log_World norot dead 756 nofield
# Beam spot (just for having a visual idea!).
# Assume r = 3*sigma to intercept 99.73% of all events, with sigma = 6.83 mm
#/sr1/command construct tubs BSpot 0 30 1 0 360 G4_Galactic 0 0 -490 log_World norot dead 800 nofield
#*/sr1/command construct tubs saveBSpot 0 20.5 1 0 360 G4_Galactic 0 0 -490 log_World norot sr1/ScintSD 800 nofield
/sr1/command construct tubs saveBSpot 0 20.5 1 0 360 G4_Galactic 0 0 -835 log_World norot sr1/ScintSD 800 nofield
# Dummy plane or target to intercept muon events.
# Use it either as (DMCP-sr1/ScintSD) - with no info on mu+ polariz., or as (target-dead) with info on mu+ polariz.
#/sr1/command construct box saveCFoil 65 60 5e-4 G4_Galactic 0 0 -199 log_World norot dead 851 nofield
/sr1/command construct tubs target 0 65 1.5 0 360 G4_Galactic 0 0 490 log_World norot dead 851 nofield
/sr1/command construct tubs saveIT 0 40 1 0 360 G4_Galactic 0 0 -300 log_World norot sr1/ScintSD 900 nofield
# Visual attributes (optional)
/sr1/command visattributes log_NPole red
/sr1/command visattributes log_SPole lightblue
#/sr1/command visattributes log_NCoil Grid_style
#/sr1/command visattributes log_SCoil Grid_style
/sr1/command visattributes log_YokeL oxsteel
/sr1/command visattributes log_TCoil1 Grid_style
/sr1/command visattributes log_TCoil2 Grid_style
/sr1/command visattributes log_TCoil3 Grid_style
/sr1/command visattributes log_TCoil4 Grid_style
/sr1/command visattributes log_BCoil1 Grid_style
/sr1/command visattributes log_BCoil2 Grid_style
/sr1/command visattributes log_BCoil3 Grid_style
/sr1/command visattributes log_BCoil4 Grid_style
/sr1/command visattributes log_Cap_p darkred
/sr1/command visattributes log_Cap_n blue_style
#/sr1/command visattributes log_Uniform invisible
################################################################################################################
# -- Setting the ELECTRIC and MAGNETIC fields --
################################################################################################################
# Use ABSOLUTE coordinates to specify the field position (i.e. with respect to GLOBAL WORLD)!
# Default field units: Magnetic - T, Electric - kV/mm (or kV for E-field maps).
# NOTE: Applying a field to an invisible log_vol makes is visible!
### Uniform Electric and Magnetic field at Spin Rotator:
#*/sr1/command globalfield UniMag_field 0. 0. 0. uniform log_Uniform 0 -0.036 0 0 0 0
#*/sr1/command globalfield UniEle_field 0. 0. 0. uniform log_Uniform 0 0 0 -0.21 0 0
#*/sr1/command globalfield UniEM_field 0. 0. 0. uniform log_Uniform 0 -0.036 0 -0.21 0 0
###15 cm gap /sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_capacitor.map log_EField 22.1
###1000-0 V /sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_capacitor.map log_EField 14.74
# E_0 = 20 keV, E = 21 kV, B = 360 G
# E_0 = 15 keV, E = 18.1865 kV, B = 311.7691 G (= sqrt(15/20))
# CONSTANT Electric and Magnetic fields (15 keV 15.7 0.0312 T; 20 keV 21.0 0.036 T) (16.6 if with B map)
#/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_const.map log_EField 16.6
#/sr1/command globalfield Mag_field 0. 0. 0.1 fromfile 3DB B3D_const.map log_BField 0.03117691
# FIELD MAPS (using symmetric E-field map +/- 500 V). E3D_cap_symm.map or E3D_cap_symm_edge.map
# E0 = 20 keV
#*/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_cap_symm_edge.map log_EField 21.0
#*/sr1/command globalfield Mag_field 0. 0. 0.1 fromfile 3DB B3D_iron2.map log_BField 0.036
# E0 = 15 keV
#/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_cap_symm_edge.map log_EField 15.8
/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_edge_rod.map log_EField 15.9
#*/sr1/command globalfield Ele_field 0. 0. 0. fromfile 3DE E3D_edge_3rods.map log_EField 15.2
/sr1/command globalfield Mag_field 0. 0. 0.1 fromfile 3DB B3D_iron2.map log_BField 0.0312
# LENS 1. Best L1 settings: 8.5 kV for 15 keV muons, 11.33 kV for 20 keV muons
/sr1/command globalfield L1E_field 0. 0. -630 fromfile 2DE L1_Erz.map log_L1ENV 0.
# Set parameters for particle tracking in an EM field
/sr1/command globalfield setparameter SetLargestAcceptableStep 5
/sr1/command globalfield setparameter SetMinimumEpsilonStep 5e-5
/sr1/command globalfield setparameter SetMaximumEpsilonStep 0.001
/sr1/command globalfield setparameter SetDeltaOneStep 0.1
/sr1/command globalfield setparameter SetDeltaIntersection 0.01
/sr1/command globalfield printparameters
# TESTING EM FIELD
/sr1/command globalfield printFieldValueAtPoint 0. 0. 0.
/sr1/command globalfield printFieldValueAtPoint 0. 35. -670.
################################################################################################################
# -- Setting simulation PARAMETERS --
################################################################################################################
# Set processes from: lowenergy, penelope, coulombAndMultiple (default, Coul. only for CFoil), coulomb (for all, very slow).
#**/sr1/command typeofprocesses coulombAndMultiple
#*/sr1/command typeofprocesses penelope
#*/sr1/command includeMuoniumProcesses false
# Set the overall range cut (default 0.1 mm)
#*/run/setCut 1 mm
# Set the range cut on particular volumes (in mm)
#*/sr1/command SetUserLimits log_target 0.01
#*/sr1/command SetUserLimits log_targetscint 0.01
#*/sr1/command SetUserLimits log_cryostatscint 0.01
# Set particle energy cuts on particular volumes (in eV)
#*/sr1/command SetUserMinEkine log_World 0.1
# Store ALL the events in a ROOT tree or just the interesting ones? (default is true)
#*/sr1/command storeOnlyEventsWithHits false
# Set the minimum time separation between two subsequent signals in the same detector (in ns)
/sr1/command signalSeparationTime 0.1
# Override runID number
#*/sr1/run/runID 21
# Set the frequency of event printing
/sr1/run/howOftenToPrintEvent 1000
# RANDOM option choices: (specify the random number generator initialisation)
# 0 ... no initialisation (default)
# 1 ... use actual computer time to initialise now # Pseudo-random numbers
# 2 ... use event number to initialise at the beginning of each event # Reproducible numbers
# 3 ... read in the random no. initial values for each event from a file
/sr1/run/randomOption 2
# VISUALIZATION options
# To enable or disable visualization uncomment one of these lines
# To modify visualization options edit the file vis.mac
/vis/disable
#*/control/execute vis.mac
#alias g4='export G4WORKDIR="/home/sedlak/bin_4.9.1"; source /home/geant4/4.9.1/env.sh;
#export G4VRMLFILE_VIEWER="vrmlview"; echo "On this machine the G4VRMLFILE_VIEWER=$G4VRMLFILE_VIEWER"'
################################################################################################################
# -- Setting PARTICLE GUN parameters --
################################################################################################################
# Default momentum direction: 001, i.e. 0z.
# Default muon spin direction: 100, i.e. 0x.
# Default particle type: mu+ (can be changed to Mu or proton)
# Set particle type
#*/gun/particle Mu
#*/gun/particle proton
/gun/particle mu+
# Set the position of the beam vertex
# For use only with SR: d = -495 mm, when L1 is included d = -835 mm.
#/gun/vertex 0. 0. -495. mm
/gun/vertex 0. 0. -835. mm
# A point-like uniform beam
#*/gun/vertexsigma -0.1 -0.1 0 mm
# Set beam transverse spread (default GAUSSIAN spread)
# If FWHM = 10 mm ==> sigma = 10/2.354 = 4.2481 mm (last 0 is a dummy value)
# Negative sigma values => random FLAT RECTANGULAR distribution (area 2x.2y)
# Use vertexboundary with (vb < sigma_xy) to obtain a CIRCULAR beam spot
# /gun/vertexsigma 0 0 0 mm ==> Very SLOW with mag. field ON and centered beam
## The LEMUSR beam has a sigma of 6.83 mm => FWHM = 2.355*6.83 = 16.08 mm
/gun/vertexsigma 6.83 6.83 0 mm
#/gun/vertexsigma -20 -20 0 mm
#/gun/vertexboundary 20 -1e6 1e6 mm
# /gun/vertexboundary: rMaxAllowed, zMinAllowed, zMaxAllowed # Beam AND gating
#*/gun/vertexboundary 7 -1314.4 -1305 mm
# Without restrictions in z, but only on r:
#*/gun/vertexboundary 3 -1e6 1e6 mm
# Set beam momentum (USE only as an ALTERNATIVE to setting energy!)
# /gun/momentum 0 0 29.79 MeV
#*/gun/momentum 0 0 1.8 MeV
# Energy loss at p = 1.2 MeV/c (E = 6.8 keV) => 1.23 +/- 0.2 keV
# Energy loss at p = 1.8 MeV/c (E = 15.3 keV) => 1.25 +/- 0.3 keV
# 1.2 MeV/c -> 6.8 keV, 1.8 MeV/c -> 15.3 keV, 2.056 MeV/c -> 20 keV.
# muon rest mass = 105.658 MeV/c2
#/gun/kenergy 3.4 MeV
/gun/kenergy 15.0 keV
# Set beam momentum direction
/gun/direction 0.0 0.0 1.0
# Set muon spin direction
/gun/muonPolarizVector 1 0 0
# Other useful test parameters:
#
# FWHM = 3% ==> sigma = 29.79*0.03/2.354 = 0.37965 MeV/c
#*/gun/momentumsmearing 0.37965 MeV
#---/gun/momentumboundary: pMinAllowed, pMaxAllowed, dummy
#*/gun/momentumboundary 20 40 0 MeV
#---/gun/tilt: xangle, yangle, dummy
#*/gun/tilt 0 0.5 0 deg
#---/gun/tiltsigma: xangleSigma, yangleSigma, dummy (1 degree at 1 m => 17 mm)
/gun/tiltsigma 1.4 1.4 0 deg
#*/gun/pitch 0.5 deg
#---/gun/decaytimelimits: decayMin, decayMax, decayTime
#*/gun/decaytimelimits 10400 10420 2197.03 ns
# Selectively inactivate or activate sensitive detectors
#*/hits/inactivate /sr1/ScintSD
# Only for code debugging!
#/tracking/verbose 1
# For running a batch see the files v1.mac and changefield.mac in /home/l_shiroka/geant4/musim_alc2
# BEAM ON
#*/run/beamOn 1000000
#*/run/beamOn 2
#/run/beamOn 1000000
/run/beamOn 10000

View File

@ -0,0 +1,14 @@
2 2 2 cm 1
% nx ny nz length_unit norm_fact
% Bx(T) By(T) Bz (T)
% A norm_fact 1 means that the magnetic field is expressed in Tesla.
%
% x(cm) y(cm) z(cm) Bx_emqv(V/m) By_emqv(V/m) Bz_emqv(V/m)
-4.9 -4.9 -15.0 0 -1 0
-4.9 -4.9 15.0 0 -1 0
-4.9 4.9 -15.0 0 -1 0
-4.9 4.9 15.0 0 -1 0
4.9 -4.9 -15.0 0 -1 0
4.9 -4.9 15.0 0 -1 0
4.9 4.9 -15.0 0 -1 0
4.9 4.9 15.0 0 -1 0

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,16 @@
2 2 2 cm 1e-4
%nx ny nz length_unit norm_fact
% Ex (V/m) Ey(V/m) Ez (V/m)
% A norm_fact 1e-6 converts V/m (FL) to kV/mm (G4).
% This FemLab model uses cm, therefore norm_fact = norm_fact*1e2.
%
% Exported from capacitor.fl
% x(cm) y(cm) z(cm) Ex_emqv(V/m) Ey_emqv(V/m) Ez_emqv(V/m)
-4.9 -4.9 -15.0 -100 0 0
-4.9 -4.9 15.0 -100 0 0
-4.9 4.9 -15.0 -100 0 0
-4.9 4.9 15.0 -100 0 0
4.9 -4.9 -15.0 -100 0 0
4.9 -4.9 15.0 -100 0 0
4.9 4.9 -15.0 -100 0 0
4.9 4.9 15.0 -100 0 0

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,191 @@
#define analysis_cxx
#include "analysis.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
void analysis::Loop()
{
// In a ROOT session, you can do:
// Root > .L analysis.C
// Root > analysis 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 <- Very Important!
// Root > .ls // Shows what histos are available for plotting
//
// 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 as:
// METHOD1:
// fChain->SetBranchStatus("*",0); // disable all branches
// fChain->SetBranchStatus("branchname",1); // activate branchname
// METHOD2: replace line
// fChain->GetEntry(jentry); //read all branches
// b_branchname->GetEntry(ientry); //read only this branch
//
//
// To create TEMPLATE files for analysis proceed as follows:
// 1. Load the root file containing all the variables, e.g.:
// TFile* f=new TFile("data/lem4_1049.root");
// 2. Quickly check the Tree and Variable names:
// f->ls(); // Shows tree name (i.e. t1)
// t1->Print(); // Shows variable names
// 3. Create the template files .C and .h (t1 is the Tree name)
// t1->MakeClass("myAnalysis")
//
//
// How to write data to an EXTERNAL file (use cout instead of outdata for printing on screen):
// See also: www.bgsu.edu/departments/compsci/docs/index.html
//
// ofstream outdata;
// outdata.open("export.dat");
// for (Int_t i=0; i<hEdepositCF->GetNbinsX(); i++) {outdata<<hEdepositCF->GetBinContent(i)<<endl;}
// outdata.close();
//
if (fChain == 0) return;
//DEFINE HISTOGRAMS
TH1F* hEdeposited = new TH1F("hEdeposited","Energy spectrum; E [MeV]", 250,0.,2.5);
TH1F* hEdepositCF = new TH1F("hEdepositCF","Energy spectrum; E [MeV]", 250,0.,3.0);
TH1F* hEdepTrig = new TH1F("hEdepTrig", "Radioactive electron kin. energy",250,0.,2.5);
TH1F* hDettUp = new TH1F("hDettUp", "Muon arrival times Up (#mus)", 250,0.,10);
TH1F* hDettDown = new TH1F("hDettDown", "Muon arrival times Down (#mus)", 250,0.,10);
TH1F* hDettRight = new TH1F("hDettRight", "Muon arrival times Right (#mus)", 250,0.,10);
TH1F* hDettLeft = new TH1F("hDettLeft", "Muon arrival times Left (#mus)", 250,0.,10);
TH1F* hEdepoTest = new TH1F("hEdepoTest", "Number of events in coincidence", 100,0.,1.);
hEdeposited->Sumw2();
hEdepositCF->Sumw2();
hEdepTrig->Sumw2();
Long64_t nentries = fChain->GetEntriesFast();
//nentries=1000;
Long64_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
// FILL IN HISTOGRAMS
// ATTENTION: change 11.27 to the starting muon energy:
hEdepositCF->Fill(11.27 + 3.73-save_ke[0]/1e3);
//for (Int_t i=0; i<det_n; i++)
// { if (det_ID[i]==623) { hEdeposited->Fill(save_ke[i]);}
// }
// Coincidence between Up scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==21) {
for (Int_t j=0; j<det_n; j++)
if (det_ID[j]==11) {
hDettUp->Fill(det_time_start[j]);
///hEdeposited->Fill(det_edep[i]);
// if (electronEscaped[i]) {
// hEdepTrig->Fill(det_edep[i]);
// }
}
}
}
// Coincidence between Down scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==23) {
for (Int_t j=0; j<det_n; j++)
if (det_ID[j]==13) {
hDettDown->Fill(det_time_start[j]);
}
}
}
// Coincidence between Right scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==22) {
for (Int_t j=0; j<det_n; j++)
if (det_ID[j]==12) {
hDettRight->Fill(det_time_start[j]);
}
}
}
// Coincidence between Left scintillators
for (Int_t i=0; i<det_n; i++) {
if (det_ID[i]==24) {
for (Int_t j=0; j<det_n; j++)
if (det_ID[j]==14) {
hDettLeft->Fill(det_time_start[j]);
}
}
}
// Testing event coincidence with SciFi
Bool_t twentyone= false;
Bool_t twozero7 = false;
Bool_t twozero8 = false;
Int_t ikk=0;
for (Int_t i=0; i<det_n; i++) {
if ((det_ID[i]== 21)&&(det_edep[i]>0.0)) {twentyone = true; ikk = i;}
if ((det_ID[i]==207)&&(det_edep[i]>0.0)) twozero7 = true;
if ((det_ID[i]==208)&&(det_edep[i]>0.0)) twozero8 = true;
}
//if (twentyone) hEdepoTest->Fill(det_edep[ikk]); // Fiber hits
//if (twozero7&& twozero8) hEdepoTest->Fill(det_edep[ikk]); // Coincidence
if (twentyone&&twozero7&& twozero8) hEdepoTest->Fill(det_edep[ikk]); // Fibre hits + coincidence
}
// PLOT HISTOGRAMS
TCanvas* c1= new TCanvas("c1","canvas 1");
// hEdeposited->Scale(1/(number_of_generated_events*hEdeposited->GetBinWidth(1)));
TH1F *htemp1 = (TH1F*) hDettUp->Clone(); htemp1->SetName("htemp1");
TH1F *htemp2 = (TH1F*) hDettUp->Clone(); htemp2->SetName("htemp2");
TH1F *htemp3 = (TH1F*) hDettUp->Clone(); htemp3->SetName("htemp3");
TH1F *htemp4 = (TH1F*) hDettUp->Clone(); htemp4->SetName("htemp4");
TH1F *htemp5 = (TH1F*) hDettUp->Clone(); htemp5->SetName("htemp5");
htemp5->SetTitle("Muon decay asymmetry U-D/R-L; Time (#mus); Asymmetry");
TH1F *htemp6 = (TH1F*) hDettUp->Clone(); htemp6->SetName("htemp6");
//htemp6->SetTitle("Muon decay asymmetry U-D; Time (#mus); Asymmetry");
htemp1->Add(hDettUp,hDettDown,1.,-1.);
htemp2->Add(hDettUp,hDettDown,1., 1.);
htemp5->Divide(htemp1,htemp2,1.,1.);
htemp3->Add(hDettLeft,hDettRight,1.,-1.);
htemp4->Add(hDettLeft,hDettRight,1., 1.);
htemp6->Divide(htemp3,htemp4,1.,1.);
htemp5->SetLineWidth(2);
htemp5->SetLineColor(2);
htemp5->SetAxisRange(-0.8,0.8,"Y");
htemp5->Draw("hist");
htemp6->SetLineWidth(2);
htemp6->SetLineColor(4);
htemp6->SetAxisRange(-0.8,0.8,"Y");
htemp6->Draw("same");
hEdepositCF->SetLineWidth(2);
hEdepositCF->SetLineColor(2);
hEdepositCF->Draw("hist");
hEdepositCF->Fit("gaus");
TCanvas* c2= new TCanvas("c2","canvas 2");
hEdepoTest->Draw();
// sprintf(filenamePrint,"musrSrClass_%s.ps",fileNr);
// c1->Print(filenamePrint);
}

View File

@ -0,0 +1,304 @@
//////////////////////////////////////////////////////////
// This class has been automatically generated on
// Fri Oct 26 16:55:39 2007 by ROOT version 5.08/00
// from TTree t1/a simple Tree with simple variables
// found on file: musr_200003.root
//////////////////////////////////////////////////////////
#ifndef analysis_h
#define analysis_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
class analysis {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
Int_t fCurrent; //!current Tree number in a TChain
// Declaration of leave types
Int_t runID;
Int_t eventID;
Double_t BFieldAtDecay_Bx;
Double_t BFieldAtDecay_By;
Double_t BFieldAtDecay_Bz;
Double_t BFieldAtDecay_B3;
Double_t BFieldAtDecay_B4;
Double_t BFieldAtDecay_B5;
Double_t muIniPosX;
Double_t muIniPosY;
Double_t muIniPosZ;
Double_t muIniMomX;
Double_t muIniMomY;
Double_t muIniMomZ;
Double_t muIniPolX;
Double_t muIniPolY;
Double_t muIniPolZ;
Int_t muDecayDetID;
Double_t muDecayPosX;
Double_t muDecayPosY;
Double_t muDecayPosZ;
Double_t muDecayTime;
Double_t muDecayPolX;
Double_t muDecayPolY;
Double_t muDecayPolZ;
Double_t muTargetTime;
Double_t muTargetPolX;
Double_t muTargetPolY;
Double_t muTargetPolZ;
Double_t fieldValue;
Int_t det_n;
Int_t det_ID[50]; //[det_n]
Double_t det_edep[50]; //[det_n]
Int_t det_nsteps[50]; //[det_n]
Double_t det_length[50]; //[det_n]
Double_t det_time_start[50]; //[det_n]
Double_t det_time_end[50]; //[det_n]
Double_t det_x[50]; //[det_n]
Double_t det_y[50]; //[det_n]
Double_t det_z[50]; //[det_n]
Int_t save_n;
Int_t save_detID[10]; //[save_n]
Int_t save_particleID[10]; //[save_n]
Double_t save_ke[10]; //[save_n]
Double_t save_x[10]; //[save_n]
Double_t save_y[10]; //[save_n]
Double_t save_z[10]; //[save_n]
Double_t save_px[10]; //[save_n]
Double_t save_py[10]; //[save_n]
Double_t save_pz[10]; //[save_n]
// List of branches
TBranch *b_runID; //!
TBranch *b_eventID; //!
TBranch *b_BFieldAtDecay; //!
TBranch *b_muIniPosX; //!
TBranch *b_muIniPosY; //!
TBranch *b_muIniPosZ; //!
TBranch *b_muIniMomX; //!
TBranch *b_muIniMomY; //!
TBranch *b_muIniMomZ; //!
TBranch *b_muIniPolX; //!
TBranch *b_muIniPolY; //!
TBranch *b_muIniPolZ; //!
TBranch *b_muDecayDetID; //!
TBranch *b_muDecayPosX; //!
TBranch *b_muDecayPosY; //!
TBranch *b_muDecayPosZ; //!
TBranch *b_muDecayTime; //!
TBranch *b_muDecayPolX; //!
TBranch *b_muDecayPolY; //!
TBranch *b_muDecayPolZ; //!
TBranch *b_muTargetTime; //!
TBranch *b_muTargetPolX; //!
TBranch *b_muTargetPolY; //!
TBranch *b_muTargetPolZ; //!
TBranch *b_fieldValue; //!
TBranch *b_det_n; //!
TBranch *b_det_ID; //!
TBranch *b_det_edep; //!
TBranch *b_det_nsteps; //!
TBranch *b_det_length; //!
TBranch *b_det_time_start; //!
TBranch *b_det_time_end; //!
TBranch *b_det_x; //!
TBranch *b_det_y; //!
TBranch *b_det_z; //!
TBranch *b_save_n; //!
TBranch *b_save_detID; //!
TBranch *b_save_particleID; //!
TBranch *b_save_ke; //!
TBranch *b_save_x; //!
TBranch *b_save_y; //!
TBranch *b_save_z; //!
TBranch *b_save_px; //!
TBranch *b_save_py; //!
TBranch *b_save_pz; //!
analysis(TTree *tree=0);
virtual ~analysis();
virtual Int_t Cut(Long64_t entry);
virtual Int_t GetEntry(Long64_t entry);
virtual Long64_t LoadTree(Long64_t entry);
virtual void Init(TTree *tree);
virtual void Loop();
virtual Bool_t Notify();
virtual void Show(Long64_t entry = -1);
};
#endif
#ifdef analysis_cxx
analysis::analysis(TTree *tree)
{
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.
if (tree == 0) {
TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("data/sr_1091.root");
if (!f) {
//f = new TFile("musr_200003.root");
f = new TFile("data/sr1_1091.root"); //1049
}
tree = (TTree*)gDirectory->Get("t1");
}
Init(tree);
}
analysis::~analysis()
{
if (!fChain) return;
delete fChain->GetCurrentFile();
}
Int_t analysis::GetEntry(Long64_t entry)
{
// Read contents of entry.
if (!fChain) return 0;
return fChain->GetEntry(entry);
}
Long64_t analysis::LoadTree(Long64_t entry)
{
// Set the environment to read one entry
if (!fChain) return -5;
Long64_t centry = fChain->LoadTree(entry);
if (centry < 0) return centry;
if (fChain->IsA() != TChain::Class()) return centry;
TChain *chain = (TChain*)fChain;
if (chain->GetTreeNumber() != fCurrent) {
fCurrent = chain->GetTreeNumber();
Notify();
}
return centry;
}
void analysis::Init(TTree *tree)
{
// The Init() function is called when the selector needs to initialize
// a new tree or chain. Typically here the branch addresses of the tree
// will be set. It is normaly not necessary to make changes to the
// generated code, but the routine can be extended by the user if needed.
// Init() will be called many times when running with PROOF.
// Set branch addresses
if (tree == 0) return;
fChain = tree;
fCurrent = -1;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("runID",&runID);
fChain->SetBranchAddress("eventID",&eventID);
fChain->SetBranchAddress("BFieldAtDecay",&BFieldAtDecay_Bx);
fChain->SetBranchAddress("muIniPosX",&muIniPosX);
fChain->SetBranchAddress("muIniPosY",&muIniPosY);
fChain->SetBranchAddress("muIniPosZ",&muIniPosZ);
fChain->SetBranchAddress("muIniMomX",&muIniMomX);
fChain->SetBranchAddress("muIniMomY",&muIniMomY);
fChain->SetBranchAddress("muIniMomZ",&muIniMomZ);
fChain->SetBranchAddress("muIniPolX",&muIniPolX);
fChain->SetBranchAddress("muIniPolY",&muIniPolY);
fChain->SetBranchAddress("muIniPolZ",&muIniPolZ);
fChain->SetBranchAddress("muDecayDetID",&muDecayDetID);
fChain->SetBranchAddress("muDecayPosX",&muDecayPosX);
fChain->SetBranchAddress("muDecayPosY",&muDecayPosY);
fChain->SetBranchAddress("muDecayPosZ",&muDecayPosZ);
fChain->SetBranchAddress("muDecayTime",&muDecayTime);
fChain->SetBranchAddress("muDecayPolX",&muDecayPolX);
fChain->SetBranchAddress("muDecayPolY",&muDecayPolY);
fChain->SetBranchAddress("muDecayPolZ",&muDecayPolZ);
fChain->SetBranchAddress("muTargetTime",&muTargetTime);
fChain->SetBranchAddress("muTargetPolX",&muTargetPolX);
fChain->SetBranchAddress("muTargetPolY",&muTargetPolY);
fChain->SetBranchAddress("muTargetPolZ",&muTargetPolZ);
fChain->SetBranchAddress("fieldValue",&fieldValue);
fChain->SetBranchAddress("det_n",&det_n);
fChain->SetBranchAddress("det_ID",det_ID);
fChain->SetBranchAddress("det_edep",det_edep);
fChain->SetBranchAddress("det_nsteps",det_nsteps);
fChain->SetBranchAddress("det_length",det_length);
fChain->SetBranchAddress("det_time_start",det_time_start);
fChain->SetBranchAddress("det_time_end",det_time_end);
fChain->SetBranchAddress("det_x",det_x);
fChain->SetBranchAddress("det_y",det_y);
fChain->SetBranchAddress("det_z",det_z);
fChain->SetBranchAddress("save_n", &save_n, &b_save_n);
fChain->SetBranchAddress("save_detID", save_detID, &b_save_detID);
fChain->SetBranchAddress("save_particleID", save_particleID, &b_save_particleID);
fChain->SetBranchAddress("save_ke", save_ke, &b_save_ke);
fChain->SetBranchAddress("save_x", save_x, &b_save_x);
fChain->SetBranchAddress("save_y", save_y, &b_save_y);
fChain->SetBranchAddress("save_z", save_z, &b_save_z);
fChain->SetBranchAddress("save_px", save_px, &b_save_px);
fChain->SetBranchAddress("save_py", save_py, &b_save_py);
fChain->SetBranchAddress("save_pz", save_pz, &b_save_pz);
Notify();
}
Bool_t analysis::Notify()
{
// The Notify() function is called when a new file is opened. This
// can be either for a new TTree in a TChain or when when a new TTree
// is started when using PROOF. Typically here the branch pointers
// will be retrieved. It is normaly not necessary to make changes
// to the generated code, but the routine can be extended by the
// user if needed.
// Get branch pointers
b_runID = fChain->GetBranch("runID");
b_eventID = fChain->GetBranch("eventID");
b_BFieldAtDecay = fChain->GetBranch("BFieldAtDecay");
b_muIniPosX = fChain->GetBranch("muIniPosX");
b_muIniPosY = fChain->GetBranch("muIniPosY");
b_muIniPosZ = fChain->GetBranch("muIniPosZ");
b_muIniMomX = fChain->GetBranch("muIniMomX");
b_muIniMomY = fChain->GetBranch("muIniMomY");
b_muIniMomZ = fChain->GetBranch("muIniMomZ");
b_muIniPolX = fChain->GetBranch("muIniPolX");
b_muIniPolY = fChain->GetBranch("muIniPolY");
b_muIniPolZ = fChain->GetBranch("muIniPolZ");
b_muDecayDetID = fChain->GetBranch("muDecayDetID");
b_muDecayPosX = fChain->GetBranch("muDecayPosX");
b_muDecayPosY = fChain->GetBranch("muDecayPosY");
b_muDecayPosZ = fChain->GetBranch("muDecayPosZ");
b_muDecayTime = fChain->GetBranch("muDecayTime");
b_muDecayPolX = fChain->GetBranch("muDecayPolX");
b_muDecayPolY = fChain->GetBranch("muDecayPolY");
b_muDecayPolZ = fChain->GetBranch("muDecayPolZ");
b_muTargetTime = fChain->GetBranch("muTargetTime");
b_muTargetPolX = fChain->GetBranch("muTargetPolX");
b_muTargetPolY = fChain->GetBranch("muTargetPolY");
b_muTargetPolZ = fChain->GetBranch("muTargetPolZ");
b_fieldValue = fChain->GetBranch("fieldValue");
b_det_n = fChain->GetBranch("det_n");
b_det_ID = fChain->GetBranch("det_ID");
b_det_edep = fChain->GetBranch("det_edep");
b_det_nsteps = fChain->GetBranch("det_nsteps");
b_det_length = fChain->GetBranch("det_length");
b_det_time_start = fChain->GetBranch("det_time_start");
b_det_time_end = fChain->GetBranch("det_time_end");
b_det_x = fChain->GetBranch("det_x");
b_det_y = fChain->GetBranch("det_y");
b_det_z = fChain->GetBranch("det_z");
return kTRUE;
}
void analysis::Show(Long64_t entry)
{
// Print contents of entry.
// If entry is not specified, print current entry
if (!fChain) return;
fChain->Show(entry);
}
Int_t analysis::Cut(Long64_t entry)
{
// This function may be called from Loop.
// returns 1 if entry is accepted.
// returns -1 otherwise.
return 1;
}
#endif // #ifdef analysis_cxx

View File

@ -0,0 +1,61 @@
% Electric field
% Read file (check the number of header lines) and assing variables.
filename = 'E3D_cap_symm.map';
[x y z Ex Ey Ez] = textread(filename,'%f %f %f %f %f %f','headerlines',8);
z_ind = find (x==0 & y==0); % Determine indices of the points along z axis.
z_p = z(z_ind); % z values along z axis
Ex_p = Ex(z_ind); % Ex field profile along z axis
plot(z_p,-Ex_p*0.8e-8,'m-');
%Calculate effective field length l_eff = Int(Ex dz) / max(Ex).
z_step = max(diff(z));
i_center = find (x==0 & y==0 & z==0); % Field center index
if Ex(i_center) < 0, % If the field at center is negative use min, otherwise max.
l_effE = sum(Ex_p)*z_step/min(Ex_p);
else
l_effE = sum(Ex_p)*z_step/max(Ex_p);
end
ylim = get(gca,'YLim');
hl1 = line([-l_effE/2 -l_effE/2],[ylim(1) ylim(2)]);
hl2 = line([ l_effE/2 l_effE/2],[ylim(1) ylim(2)]);
set(hl1,'Color','m'); set(hl2,'Color','m');
text(0.75,0.8,['Eff. length E = ',num2str(l_effE,3), ' cm'],'Units','Normal','Color','m')
xlabel('z coordinate (cm)','FontSize',14)
ylabel('EM Field (arb. unit)','FontSize',14)
title('Field profile along z axis','FontSize',16)
hold on
% Magnetic field
filename = 'B3D_iron1.map';
[x y z Bx By Bz] = textread(filename,'%f %f %f %f %f %f','headerlines',5);
z_ind = find (x==0 & y==0); % Determine indices of the points along z axis.
z_p = z(z_ind); % z values along z axis
By_p = By(z_ind); % By field profile along z axis
plot(z_p,-By_p,'b-');
%Calculate effective field length l_eff = Int(By dz) / max(By).
z_step = max(diff(z));
i_center = find (x==0 & y==0 & z==0); % Field center index
if By(i_center) < 0, % If the field at center is negative use min, otherwise max.
l_effB = sum(By_p)*z_step/min(By_p);
else
l_effB = sum(By_p)*z_step/max(By_p);
end
ylim = get(gca,'YLim');
hl3 = line([-l_effB/2 -l_effB/2],[ylim(1) ylim(2)]);
hl4 = line([ l_effB/2 l_effB/2],[ylim(1) ylim(2)]);
text(0.75,0.75,['Eff. length B = ',num2str(l_effB,3), ' cm'],'Units','Normal','Color','b')
text(0.03,0.8,{'C-magnet 20x40 cm poles','Symm. capacitor voltages'},'Units','Normal','FontSize',12)
set(gca,'FontSize',14)
axis([-50 50 0 8e-7])

View File

@ -0,0 +1,32 @@
[x y z bx by bz] = textread('sep_0nl.lis','%f %f %f %f %f %f', 'headerlines', 13);
[nx ny nz] = textread('sep_0nl.lis','%*s %*s %*s %*s %u %u %u', 1, 'headerlines', 10);
Bx = reshape(bx,nx,ny,nz);
By = reshape(by,nx,ny,nz);
Bz = reshape(bz,nx,ny,nz);
x1 = reshape(x,nx,ny,nz);
y1 = reshape(y,nx,ny,nz);
z1 = reshape(z,nx,ny,nz);
i=1:32;
test = [x1(i);y1(i);z1(i)]'
%min(x) max(x)
sx = max(diff(x)); % read step size along x
sy = max(diff(y)); % read step size along y
sz = max(diff(z)); % read step size along z
x2 = -x(end:-1:1);
x3 = -x(end:-1:1);
x4 = -x(end:-1:1);
x2 = reshape(x2,nx,ny,nz);
x3 = reshape(x3,nx,ny,nz);
x4 = reshape(x4,nx,ny,nz);
%p = [x(:)'; y(:)'; z(:)'];
%ind2sub
%[i j] = ind2sub(3,x1);

View File

@ -0,0 +1,36 @@
// From www-pnp.physics.ox.ac.uk/~west/root/plotting_with_style.html
// See also http://hep-www.px.tsukuba.ac.jp/~yuji/root/MyRootTips.html
void rootlogon ()
{
printf("rootlogon.C has been loaded.\n");
gStyle->SetCanvasColor(10); // or (kWhite) background is no longer mouse-dropping white
gStyle->SetPalette(1,0); // blue to red false color palette. Use 9 for b/w
gStyle->SetCanvasBorderMode(0); // turn off canvas borders
gStyle->SetPadBorderMode(0);
gStyle->SetPaintTextFormat("5.2f"); // What precision to put numbers if plotted with "TEXT"
gStyle->SetTitleColor(1);
//gStyle->SetFillColor(0);
//gStyle->SetPadColor(10);
gStyle->SetStatColor(10); // sets the filling of statistic box to white
gStyle->SetStatFont(62);
gStyle->SetStatBorderSize(0);
//gStyle->SetPaperSize(20,24);
//gStyle->SetLabelColor(6,"xy");
gStyle->SetTitleBorderSize(0);
gStyle->SetTitleFillColor(10);
//gStyle->SetTitleColor(4,"xy");
//gStyle->SetPadColor(23);
// For publishing:
gStyle->SetLineWidth(2.);
//gStyle->SetTextSize(1.1);
gStyle->SetLabelSize(0.05,"xy");
gStyle->SetTitleSize(0.05,"xy");
gStyle->SetTitleOffset(1.1,"x");
gStyle->SetTitleOffset(1.0,"y");
gStyle->SetPadTopMargin(0.1);
gStyle->SetPadRightMargin(0.1);
gStyle->SetPadBottomMargin(0.16);
gStyle->SetPadLeftMargin(0.12);
}

16
geant4/spin_rot/run/run_jobs Executable file
View File

@ -0,0 +1,16 @@
sr1 1100.mac > data/1100.txt
sr1 1101.mac > data/1101.txt
sr1 1102.mac > data/1102.txt
sr1 1103.mac > data/1103.txt
sr1 1104.mac > data/1104.txt
sr1 1105.mac > data/1105.txt
sr1 1106.mac > data/1106.txt
sr1 1107.mac > data/1107.txt
sr1 1108.mac > data/1108.txt
sr1 1109.mac > data/1109.txt
sr1 1110.mac > data/1110.txt
sr1 1111.mac > data/1111.txt
sr1 1112.mac > data/1112.txt
sr1 1113.mac > data/1113.txt
sr1 1114.mac > data/1114.txt
sr1 1115.mac > data/1115.txt

View File

@ -0,0 +1,161 @@
void spin_rot()
{
TFile* f1 = new TFile("data/sr1_1091.root");
TCanvas* c1 = new TCanvas("c1","Beam Spot",60,40,800,800);
gStyle->SetPalette(1,0);TCanvas* c1 = new TCanvas("c1","Beam Spot",60,40,800,800);
c1->Divide(2,2);
c1->SetGridx();
c1->SetGridy();
//c1->cd(1)->SetGridx();c1->cd(1)->SetGridy();
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("det_y[0]:det_x[0]>>bspot_hist(64,-40,40,64,-40,40)");
//t1->Draw("det_y:det_x>>bspot_hist(64,-80,80,64,-80,80)"); //t1->Draw("mcpHit.positiony*10:mcpHit.positionx*10>>bspot_hist(64,-40,40,64,-40,40)","mcpHit.positionz==2&&mcpHit.ID==0");
bspot_hist->SetTitle("Muon beam cross section [mm]");
//bspot_hist->Draw("cont0");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
// t1->Draw("mcpHit.positionx*10:mcpHit.positiony*10>>bspot_hist(128,-12,12,128,-12,12)","","goff");
bspot_hist->ProjectionX("xsect",31,33);
bspot_hist->ProjectionY("ysect",31,33);
//pl=new TPaveLabel(2, 155, 16, 175, "RA_Up = 11 kV");
pl=new TPaveLabel(-75, 60, -15, 75, "#splitline{RA_Up = 8.6 kV}{L3 = 6.78 kV}");
pl->SetBorderSize(1);
pl->SetTextSize(0.43);
pl->Draw();
c1->cd(3);
ysect->SetTitle("y cross section [mm]");
ysect->SetFillStyle(1001);
//ysect->Draw("F");
ysect->SetFillColor(45);
ysect->Draw();
c1->cd(1);
t1->Draw("save_x","(save_detID==800)"); //xsect->SetTitle("x cross section [mm]");
htemp->SetTitle("Position 1 [mm]");
htemp->Draw("F");
htemp->SetFillStyle(1001);
htemp->SetFillColor(38);
htemp->Draw();
c1->cd(2);
//Energy loss in CFoil
//t1->Draw("10-det_edep_mup*1000","(det_n==1)&&(det_edep<0.1)&&(det_edep_mup*1000>0.1)");
t1->Draw("save_x","(save_detID==900)");
//t1->Draw("det_edep_mup*1000","(det_n==1)&&(det_edep<0.1)");
//t1->Draw("det_edep_mup*1000>>tmp_hist(64,-80,80)","(det_n==1)&&(det_edep<0.1)");
htemp->SetTitle("Position 2 [mm]");
htemp->Draw("F");
htemp->SetFillStyle(1001);
//htemp->SetFillColor(kBlue-5);
htemp->SetFillColor(kGreen-5);
htemp->Draw();
/// Plot of muon beam angular dispersion after lens 1:
TCanvas* c2 = new TCanvas("c2","Angular Dispersion",60,40,800,800);
//gStyle->SetPalette(1,0);
//c1->Divide(2,2);
c2->SetGridx();
c2->SetGridy();
//c2->SetFillColor(kWhite); // Fill background with white for printing
//c2->cd(1);
t1->Draw("save_pz","(save_detID==800)"); //xsect->SetTitle("x cross section [mm]");
htemp->SetTitle("Position 1 [mm]");
htemp->Draw("F");
htemp->SetFillStyle(1001);
htemp->SetFillColor(38);
htemp->Draw();
//t1->Draw("atan2(mcpHit.posyini,mcpHit.posxini)*57.29577951308:atan2(mcpHit.positiony,mcpHit.positionx)*57.29577951308>
// Double_t Pi = TMath::Pi(); // 180/Pi
//c2->cd(2)->SetGridx(); c1->cd(2)->SetGridy();
//t1->Draw("atan2(save_py,save_px)*57.29577951308:atan2(mcpHit.positiony,mcpHit.positionx)*57.29577951308>>bspot_hist3(64,-180,180,64,-180,180)","mcpHit.ID==0");
///t1->Draw("atan2(save_py,save_px)*57.29577951308");
//t1->Draw("det_y[0]:det_x[0]>>bspot_hist(64,-40,40,64,-40,40)");
t1->Draw("90-atan2(save_pz,sqrt(save_px^2 + save_py^2))*57.29577951308>>ang_hist(64,0,5)","(save_detID==900)");
//t1->Draw("90-atan2(save_pz,sqrt(save_px^2 + save_py^2))*57.29577951308");
//htemp->SetTitle("Correlation angle for L1 (8.7 kV, 15 keV); Final angle [deg.]; Inital angle [deg.]");
ang_hist->GetXaxis()->SetRangeUser(-5., 5.);
ang_hist->Draw("F");
ang_hist->SetFillStyle(1001);
ang_hist->SetFillColor(805);
ang_hist->SetTitleSize(0.025);
ang_hist->SetTitle("Muon angular distribution after L1");
ang_hist->GetXaxis()->SetNdivisions(206);
ang_hist->GetXaxis()->SetTickLength(0.018);
ang_hist->GetXaxis()->SetLabelSize(0.04);
ang_hist->GetXaxis()->SetTitleSize(0.045);
ang_hist->GetYaxis()->SetTitleOffset(1.0);
ang_hist->GetXaxis()->SetTitle("Angle (deg.)");
TAxis* yax = ang_hist->GetYaxis();
yax->SetNdivisions(206); // N2*100 + N1 (N2 = secondary, N1 = primary) //yax->SetNdivisions(410);
yax->SetTickLength(0.018);
yax->SetLabelSize(0.04);
yax->SetTitleSize(0.045);
yax->SetTitleOffset(1.3);
yax->SetTitle("Counts (arb. units)");
ang_hist->Draw();
TText *pt = new TText(3.2, 200., "Plane at -300 mm");
pt->SetTextSize(0.03);
pt->Draw();
//TPaveText *pt = new TPaveText(0.6,0.4,0.89, 0.45,"trNDC");
//pt->SetTextSize(0.03);
//pt->SetFillColor(0); //pt->SetFillColor(390);
//pt->SetBorderSize(0);
//pt->SetTextAlign(12);
//pte = pt->AddText("Plane at z = -300 mm");
//pt->Draw();
//c2->SaveAs("data/sr1_ang_dist_300.eps");
}
/// Insert simple text directly from command line
// TText text(3, 300, "Plane at -380 mm");
// text.SetTextSize(0.03);
// text.Draw();
/// Optional changes of statistics
// TPaveStats* sb2=(TPaveStats *)(ang_hist->GetListOfFunctions()->FindObject("stats"));
// sb2->SetTextSize(0.025);
// sb2->SetFillColor(390);
/// An easy way to plot functions
// (new TF1("fun1","(1-x)*sqrt(1-x*x)",-1,1))->Draw();
// (new TF2("fun2","(175/80)**2*(1-y)*(1-x**2)+(1+y)*(1-x)**2",-1,1,-1,1))->Draw("lego")
// From: www-pnp.physics.ox.ac.uk/~west/root/plotting_with_style.html#demo_background
/// For publishing (insert these lines in rootlogon.C):
// gStyle->SetLineWidth(2.);
// gStyle->SetTextSize(1.1);
// gStyle->SetLabelSize(0.06,"xy");
// gStyle->SetTitleSize(0.06,"xy");
// gStyle->SetTitleOffset(1.2,"x");
// gStyle->SetTitleOffset(1.0,"y");
// gStyle->SetPadTopMargin(0.1);
// gStyle->SetPadRightMargin(0.1);
// gStyle->SetPadBottomMargin(0.16);
// gStyle->SetPadLeftMargin(0.12);

View File

@ -0,0 +1,198 @@
// To run root without logon screen and script use: root -l -n
// To run a named macro in root use:
//.L sr_transm.C
// sr_transm("data/sr1_1100.root")
#include <string>
void sr_transm(char* fname)
//void sr_transm() /// Uncomment this and next line to run with a FIXED file name
{
//TFile* f1 = new TFile("data/sr1_1091.root");
TFile* f1 = new TFile(fname);
//TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,800,800);
gStyle->SetPalette(1,0);
TCanvas* c1 = new TCanvas("c1","Polarization and Transmission",60,40,800,800);
c1->Divide(2,2);
c1->SetGridx();
c1->SetGridy();
//c1->cd(1)->SetGridx();c1->cd(1)->SetGridy();
c1->cd(4)->SetGridx();c1->cd(4)->SetGridy(); // c1->cd(3);
t1->Draw("muTargetPosY[0]:muTargetPosX[0]>>bspot_hist(64,-40,40,64,-40,40)","muTargetPolX > -2");
//t1->Draw("det_y:det_x>>bspot_hist(64,-80,80,64,-80,80)"); //t1->Draw("mcpHit.positiony*10:mcpHit.positionx*10>>bspot_hist(64,-40,40,64,-40,40)","mcpHit.positionz==2&&mcpHit.ID==0");
bspot_hist->SetTitle("Muon beam cross section [mm]");
//bspot_hist->Draw("cont0");
bspot_hist->Draw();
bspot_hist->Draw("cont0 same");
// t1->Draw("mcpHit.positionx*10:mcpHit.positiony*10>>bspot_hist(128,-12,12,128,-12,12)","","goff");
bspot_hist->ProjectionX("xsect",31,33);
bspot_hist->ProjectionY("ysect",31,33);
TText *pt = new TText(-36., 33., "L1 = 11.0 kV");
pt->SetTextSize(0.05);
pt->Draw();
//pl=new TPaveLabel(2, 155, 16, 175, "RA_Up = 11 kV");
//pl=new TPaveLabel(-75, 60, -15, 75, "#splitline{RA_Up = 8.6 kV}{L3 = 6.78 kV}");
//pl->SetBorderSize(1);
//pl->SetTextSize(0.43);
//pl->Draw();
c1->cd(3);
xsect->SetTitle("x cross section [mm]");
xsect->SetFillStyle(1001);
//xsect->Draw("F");
xsect->SetFillColor(45);
xsect->Draw();
c1->cd(2);
t1->Draw("muTargetPolZ >> pol_z(128, -1, -0.95)","(muTargetPolZ > -2)"); //-1., 1
pol_z->SetTitle("Polarization Z");
pol_z->GetXaxis()->SetNdivisions(405);
pol_z->GetYaxis()->SetNdivisions(406);
pol_z->GetXaxis()->SetTickLength(0.018);
pol_z->GetYaxis()->SetTickLength(0.018);
pol_z->Draw("F");
pol_z->SetFillStyle(1001);
pol_z->SetFillColor(38);
pol_z->Draw();
c1->cd(1);
//Energy loss in CFoil
//t1->Draw("10-det_edep_mup*1000","(det_n==1)&&(det_edep<0.1)&&(det_edep_mup*1000>0.1)");
t1->Draw("muTargetPolX >> pol_x(128, -0.4, 0)","(muTargetPolX > -2)"); //0, 2
//t1->Draw("det_edep_mup*1000","(det_n==1)&&(det_edep<0.1)");
//t1->Draw("det_edep_mup*1000>>tmp_hist(64,-80,80)","(det_n==1)&&(det_edep<0.1)");
pol_x->SetTitle("Polarization X");
pol_x->GetXaxis()->SetNdivisions(405);
pol_x->GetYaxis()->SetNdivisions(406);
pol_x->GetXaxis()->SetTickLength(0.018);
pol_x->GetYaxis()->SetTickLength(0.018);
pol_x->Draw("F");
pol_x->SetFillStyle(1001);
//pol_x->SetFillColor(kBlue-5);
pol_x->SetFillColor(kGreen-5);
pol_x->Draw();
/*
TAxis* yax = ang_hist->GetYaxis();
yax->SetNdivisions(206); // N2*100 + N1 (N2 = secondary, N1 = primary) //yax->SetNdivisions(410);
yax->SetTickLength(0.018);
yax->SetLabelSize(0.04);
yax->SetTitleSize(0.045);
yax->SetTitleOffset(1.3);
yax->SetTitle("Counts (arb. units)");
*/
c1->SaveAs("data/sr_bspot_0.eps");
}
/*
/// Plot of muon beam angular dispersion after lens 1:
TCanvas* c2 = new TCanvas("c2","Angular Dispersion",60,40,800,800);
//gStyle->SetPalette(1,0);
//c1->Divide(2,2);
c2->SetGridx();
c2->SetGridy();
//c2->SetFillColor(kWhite); // Fill background with white for printing
//c2->cd(1);
t1->Draw("save_pz","(save_detID==800)"); //xsect->SetTitle("x cross section [mm]");
htemp->SetTitle("Position 1 [mm]");
htemp->Draw("F");
htemp->SetFillStyle(1001);
htemp->SetFillColor(38);
htemp->Draw();
//t1->Draw("atan2(mcpHit.posyini,mcpHit.posxini)*57.29577951308:atan2(mcpHit.positiony,mcpHit.positionx)*57.29577951308>
// Double_t Pi = TMath::Pi(); // 180/Pi
//c2->cd(2)->SetGridx(); c1->cd(2)->SetGridy();
//t1->Draw("atan2(save_py,save_px)*57.29577951308:atan2(mcpHit.positiony,mcpHit.positionx)*57.29577951308>>bspot_hist3(64,-180,180,64,-180,180)","mcpHit.ID==0");
///t1->Draw("atan2(save_py,save_px)*57.29577951308");
//t1->Draw("det_y[0]:det_x[0]>>bspot_hist(64,-40,40,64,-40,40)");
t1->Draw("90-atan2(save_pz,sqrt(save_px^2 + save_py^2))*57.29577951308>>ang_hist(64,0,5)","(save_detID==900)");
//t1->Draw("90-atan2(save_pz,sqrt(save_px^2 + save_py^2))*57.29577951308");
//htemp->SetTitle("Correlation angle for L1 (8.7 kV, 15 keV); Final angle [deg.]; Inital angle [deg.]");
ang_hist->GetXaxis()->SetRangeUser(-5., 5.);
ang_hist->Draw("F");
ang_hist->SetFillStyle(1001);
ang_hist->SetFillColor(805);
ang_hist->SetTitleSize(0.025);
ang_hist->SetTitle("Muon angular distribution after L1");
ang_hist->GetXaxis()->SetNdivisions(206);
ang_hist->GetXaxis()->SetTickLength(0.018);
ang_hist->GetXaxis()->SetLabelSize(0.04);
ang_hist->GetXaxis()->SetTitleSize(0.045);
ang_hist->GetYaxis()->SetTitleOffset(1.0);
ang_hist->GetXaxis()->SetTitle("Angle (deg.)");
TAxis* yax = ang_hist->GetYaxis();
yax->SetNdivisions(206); // N2*100 + N1 (N2 = secondary, N1 = primary) //yax->SetNdivisions(410);
yax->SetTickLength(0.018);
yax->SetLabelSize(0.04);
yax->SetTitleSize(0.045);
yax->SetTitleOffset(1.3);
yax->SetTitle("Counts (arb. units)");
ang_hist->Draw();
TText *pt = new TText(3.2, 200., "Plane at -300 mm");
pt->SetTextSize(0.03);
pt->Draw();
//TPaveText *pt = new TPaveText(0.6,0.4,0.89, 0.45,"trNDC");
//pt->SetTextSize(0.03);
//pt->SetFillColor(0); //pt->SetFillColor(390);
//pt->SetBorderSize(0);
//pt->SetTextAlign(12);
//pte = pt->AddText("Plane at z = -300 mm");
//pt->Draw();
//c2->SaveAs("data/sr1_ang_dist_300.eps");
}
*/
/// Insert simple text directly from command line
// TText text(3, 300, "Plane at -380 mm");
// text.SetTextSize(0.03);
// text.Draw();
/// Optional changes of statistics
// TPaveStats* sb2=(TPaveStats *)(ang_hist->GetListOfFunctions()->FindObject("stats"));
// sb2->SetTextSize(0.025);
// sb2->SetFillColor(390);
/// An easy way to plot functions
// (new TF1("fun1","(1-x)*sqrt(1-x*x)",-1,1))->Draw();
// (new TF2("fun2","(175/80)**2*(1-y)*(1-x**2)+(1+y)*(1-x)**2",-1,1,-1,1))->Draw("lego")
// From: www-pnp.physics.ox.ac.uk/~west/root/plotting_with_style.html#demo_background
/// For publishing (insert these lines in rootlogon.C):
// gStyle->SetLineWidth(2.);
// gStyle->SetTextSize(1.1);
// gStyle->SetLabelSize(0.06,"xy");
// gStyle->SetTitleSize(0.06,"xy");
// gStyle->SetTitleOffset(1.2,"x");
// gStyle->SetTitleOffset(1.0,"y");
// gStyle->SetPadTopMargin(0.1);
// gStyle->SetPadRightMargin(0.1);
// gStyle->SetPadBottomMargin(0.16);
// gStyle->SetPadLeftMargin(0.12);

View File

@ -0,0 +1,204 @@
% Test of Spin Rotator transmittance as a function of Lens L1 voltage
% Simulation performed with 1e4 events (9908), B = 360 G, E = 21 kV/m
% To get count number use, e.g.: less sr_bspot_*.eps | grep bspot_histEntries
d0 = [ ... % Parallel beam = > alpha = 0
% L1 (kV) Counts
0.0 8684
1.0 8702
2.0 8755
3.0 8780
4.0 8847
5.0 8927
6.0 8983
7.0 8985
8.0 8982
9.0 8979
9.5 8976
10.0 8964
10.5 8938
11.0 8869
11.3 8758
12.0 8395
13.0 7522
14.0 6345
15.0 5305
];
da = [... % alpha = 1 deg., alpha = 2 deg., third column SR OFF, alpha = 1 deg.
0.0 5777 3688 8794
1.0 5687 3781 8791 % 9005 for alpha = 0
2.0 5780 3790 8804
3.0 5892 3782 8825
4.0 5919 3846 8853
5.0 5986 3881 8898
6.0 6200 3925 8944 % 8991 for alpha = 0
7.0 6453 4192 8959
8.0 6724 4340 8975
9.0 7058 4542 8981
10.0 7337 4889 8981
11.0 7452 5291 8964 % 8972 for alpha = 0
12.0 7322 5663 8881
13.0 6896 5835 8593
14.0 6057 5497 7975
15.0 5082 4827 6911 % 7034 for alpha = 0
];
L1_V0 = d0(:,1); counts0 = d0(:,2); % alpha = 0 deg.
L1_V1 = da(:,1); counts1 = da(:,2); % alpha = 1 deg.
L1_V2 = da(:,1); counts2 = da(:,3); % alpha = 2 deg.
L1_V3 = da(:,1); counts3 = da(:,4); % SR OFF! alpha = 1 deg.
counts0 = counts0/100; counts1 = counts1/100;
counts2 = counts2/100; counts3 = counts3/100;
hp = plot(L1_V0, counts0,'bo-', L1_V1, counts1,'mo-', L1_V2, counts2,'ro-', L1_V3, counts3,'ko-');
set(hp, 'MarkerFaceColor','w','MarkerSize',8)
legend('Ang. spread = 0 deg.','Ang. spread = 1 deg.','Ang. spread = 2 deg.',' SR off, Ang. = 0 deg.',4)
axis([0 15 0 100])
grid on
xlabel('Lens L1 voltage (kV)','FontSize',14)
ylabel('Spin Rotator transmission (%)','FontSize',14)
text(2, 15, [{'E = 21 kV, B = 360 G'},{'10 cm symm. capacitor'},{'E_0 = 20 keV'}],'FontSize',14)
set(gca,'FontSize',14)
% Simulations as a function of L1 voltage for E0 = 15 keV and 20 keV,
% using alpha = 1.4 deg. from previous G3 simulations.
de = [...%E0 = 15, E0 = 20 keV, alpha = 1.4 deg. fixed. Last col. SR OFF at 15 keV.
0.0 4990 5065 7635
1.0 4974 5063 7665
2.0 5056 5048 7722
3.0 5023 5122 7840
4.0 5229 5190 8009
5.0 5457 5269 8253
6.0 5774 5415 8470
7.0 6174 5607 8639
8.0 6716 5784 8714
9.0 6964 6107 8642
10.0 6733 6462 8198
11.0 5628 6790 7043
12.0 4510 6923 5233
13.0 3347 6745 3542
14.0 2431 6077 2175
15.0 1568 5355 1268
];
L1_V = de(:,1); % L1 Voltage values (alpha = 1.4 deg.)
counts1 = de(:,2); counts4 = de(:,4); % E0 = 15 keV, SR ON/OFF
counts2 = de(:,3); % E0 = 20 keV, SR ON
% For SR OFF at 20 keV and alpha = 1 deg. see the variable counts3 above
counts1 = counts1/100; counts2 = counts2/100; counts4 = counts4/100;
figure(2)
hp = plot(L1_V, counts1,'bv-', L1_V, counts2,'mv-', L1_V, counts4,'ro-', L1_V3, counts3,'ko-');
set(hp, 'MarkerFaceColor','w','MarkerSize',8)
legend('E_0 = 15 keV','E_0 = 20 keV','SR off (15 keV)','SR off (20 keV)',4)
axis([0 15 0 100])
grid on
% E_0 = 20 keV, E = 21 kV, B = 360 G
% E_0 = 15 keV, E = 18.1865 kV, B = 311.7691 G
xlabel('Lens L1 voltage (kV)','FontSize',14)
ylabel('Spin Rotator transmission (%)','FontSize',14)
text(1, 15, [{'E_0 = 20 keV, E = 21 kV, B = 360 G'},...
{'E_0 = 15 keV, E = 16.3 kV, B = 312 G'},...
{'10 cm symm. capacitor'},{'Ang. spread = 1.4 deg.'}],'FontSize',14)
set(gca,'FontSize',14)
% Simulations as a function of L1 voltage for E0 = 15 keV and 20 keV,
% using alpha = 1.4 deg. and CONSTANT E and B fields.
dc = [...%E0 = 15: a) EB const; b) only E const. (19.3 kV). c) E0 = 20 keV, EB const,
% d) E0 = 15 keV using E and B MAPS with improved E field map (E3D_cap_symm_edge.map instead of E3D_cap_symm.map)
% e) E0 = 15 keV using E and B MAPS with improved E field map (E3D_edge_rod.map)
% f) E0 = 15 keV using E and B MAPS with improved E field map (E3D_edge_3rods.map) <= BEST E-FIELD MAP!
% a b c d e f
0.0 8065 7192 8082 5772 5669 7250
1.0 8102 7244 8197 5847 5612 7278
2.0 8190 7272 8243 5847 5646 7293
3.0 8257 7365 8190 5915 5825 7374
4.0 8379 7540 8301 6082 5906 7495
5.0 8527 7753 8445 6357 6219 7709
6.0 8659 8016 8541 6673 6515 7965
7.0 8762 8299 8669 7114 6956 8290
8.0 8783 8503 8874 7578 7411 8527
9.0 8747 8497 8792 7738 7653 8577
10.0 8525 8134 8920 7392 7243 8319
11.0 7790 7194 8945 6466 6291 7553
12.0 6301 5714 8889 5215 5057 6225
13.0 4443 4090 8762 4109 3650 4646
14.0 2901 2641 8326 2874 2377 3127
15.0 1674 1607 7605 1834 1485 1982
];
L1_V = dc(:,1); % L1 Voltage values (alpha = 1.4 deg.)
counts1 = dc(:,2); counts2 = dc(:,3); % E0 = 15 keV, EB const, E const only
counts3 = dc(:,4); % E0 = 20 keV, EB const
counts4 = dc(:,5); % E0 = 15 keV, using E and B maps with improved E field map (E3D_cap_symm_edge.map)
counts7 = dc(:,6); % E0 = 15 keV, using E and B maps with improved E field map (E3D_edge_rod.map)
counts8 = dc(:,7); % E0 = 15 keV, using E and B maps with improved E field map (E3D_edge_3rods.map) <- BEST!
% For SR OFF at 20 keV and alpha = 1 deg. see the variable counts3 above
counts1 = counts1/100; counts2 = counts2/100; counts3 = counts3/100;
counts4 = counts4/100; counts5 = de(:,2)/100; counts6 = de(:,4)/100;
counts7 = counts7/100; counts8 = counts8/100;
figure(3)
hp3 = plot(L1_V, counts1,'bv-', L1_V, counts2,'mv-', L1_V, counts4,'yv-', ...
L1_V, counts5,'gv-', L1_V, counts6,'ro-', L1_V3, counts3,'kd-');
set(hp3, 'MarkerFaceColor','w','MarkerSize',8)
legend('Const. E, B (15 keV)','Const. E (15 keV)','Map1 E, B (15 keV)','Map0 E, B (15 keV)','SR off (15 keV)','Const. E, B (20 keV)',3)
axis([0 15 0 100])
grid on
% E_0 = 20 keV, E = 21 kV, B = 360 G
% E_0 = 15 keV, E = 18.1865 kV, B = 311.7691 G
xlabel('Lens L1 voltage (kV)','FontSize',14)
ylabel('Spin Rotator transmission (%)','FontSize',14)
title('Spin Rot. transmission - Constant field vs. Field maps at 15 keV','FontSize',16)
text(1, 15, [{'E_0 = 20 keV, E = 21 kV, B = 360 G'},...
{'E_0 = 15 keV, E = 16.6 kV, B = 312 G'},...
{'12 cm symm. capacitor'},{'Ang. spread = 1.4 deg.'}],'FontSize',14)
set(gca,'FontSize',14)
figure(4)
hp4 = plot(L1_V, counts1,'bv-', L1_V, counts2, 'mv-', L1_V, counts8,'kv-', ...
L1_V, counts7,'rv-', L1_V, counts4,'yv-', L1_V, counts5,'gv-');
set(hp4, 'MarkerFaceColor','w','MarkerSize',8)
legend('Const. E, B fields','Const. E field, B map','Map3 E "3 rods"','Map2 E "1 rod"','Map1 E "12 cm plates"','Map0 E "10 cm plates"',3)
axis([0 15 0 100])
grid on
xlabel('Lens L1 voltage (kV)','FontSize',14)
ylabel('Spin Rotator transmission (%)','FontSize',14)
title('Spin Rot. transmission - E-field const. vs. maps at 15 keV','FontSize',16)
text(6.8, 12, [{'E_0 = 15 keV, E = 15.8 kV, B = 312 G'},...
{'12 cm symm. capacitor'},{'Ang. spread = 1.4 deg.'}],'FontSize',14)
set(gca,'FontSize',14)
return
% SR TOT transm = SR_transm * L1 transm. (SR OFF) => SR_transm. = SR TOT transm / L1 transm. (SR OFF)
hold on
hp1 = plot(L1_V, counts2./counts4*100,'go-');
set(hp1, 'MarkerFaceColor','w','MarkerSize',8)
hp2 = plot(L1_V, da(:,2)./counts4,'yo-');
set(hp2, 'MarkerFaceColor','w','MarkerSize',8)
text(1,82,'L1 transm.','FontSize',12)
text(1,68,'SR transm.','FontSize',12)
text(1,53,'Tot. transm.','FontSize',12)

View File

@ -0,0 +1,44 @@
% Transmission of Spin Rotator as a function of angular dispersion of incoming muons.
% Settings B = 360 G, E = 21 kV.
d = [...
0.0 89.0
0.2 88.8
0.4 86.8
0.7 81.0
1.0 74.1
1.5 62.6
2.0 51.7
2.5 45.6
3.0 42.0
3.5 37.71
4.0 34.03
];
d1 = [...
0.0 92.99
0.2 92.58
0.4 91.80
0.7 93.08
1.0 92.85
1.5 87.82
2.0 75.97
2.5 62.93
3.0 49.73
3.5 40.12
4.0 34.66
];
ang_sigma = d(:,1); sr_trans = d(:,2);
ang_sigm1 = d1(:,1); sr_trans1 = d1(:,2);
hp = plot(ang_sigma, sr_trans,'ro-');
hold on
hp1= plot(ang_sigm1, sr_trans1,'bo-');
set(hp, 'MarkerFaceColor','w','MarkerSize',8)
set(hp1,'MarkerFaceColor','w','MarkerSize',8)
legend('SR ON','SR OFF',3)
axis([0 4 0 100])
grid on
xlabel('Muon beam angular dispersion \sigma (deg.)','FontSize',14)
ylabel('Spin Rotator transmission (%)','FontSize',14)
text(2.5, 85, [{'E = 21 kV, B = 360 G'},{'10 mm symm. capacitor'},{'p rejection: total'}],'FontSize',14)
set(gca,'FontSize',14)

View File

@ -0,0 +1,59 @@
# 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
# 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/lightsThetaPhi 120 60
#/vis/viewer/set/hiddenEdge true
#/vis/viewer/set/style surface
#/vis/viewer/zoom 0.5
# 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
/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 mu+ 1 0 1 1
/vis/modeling/trajectories/drawByParticleID-0/setRGBA e+ 0 0 0.8 0.5
/vis/modeling/trajectories/drawByParticleID-0/setRGBA nu_e 0.7 0.7 0 1
/vis/modeling/trajectories/drawByParticleID-0/setRGBA anti_nu_mu 0.3 1.0 0 0.5
# Verbosity of hits
#/hits/verbose 2
# Output just the detector geometry
/vis/viewer/flush