diff --git a/Generators/FaserCosmicGenerator/README.md b/Generators/FaserCosmicGenerator/README.md
index 8ba2dfbc3e2b48000bffb0703aece091684e91de..3490a12099ef825fab6fea08fa4fcca8e4b13f90 100644
--- a/Generators/FaserCosmicGenerator/README.md
+++ b/Generators/FaserCosmicGenerator/README.md
@@ -18,12 +18,14 @@ Options available through the `CosmicSampler` constructor include:
* targetDepthMeters - the depth of the detector in meters (default: 85)
* thetaMinDegree - the minimum zenith angle to sample (default: 0 = vertically downward)
* thetaMaxDegree - the maximum zenith angle to sample (default: 80 = almost horizontal)
-* targetRadiusMm - the radius of the sphere around the detector on which to generate muons (default: 3500 = 3.5 m)
-* x0Mm - X offset of the target sphere from the origin in detector coordinates (default: 0)
-* y0Mm - Y offset of the target sphere from the origin in detector coordinates (default: 0)
-* z0Mm - Z offset of the target sphere from the origin in detector coordinates (default: 150 = 15 cm)
+* targetDxMm - the width of the box to generate muons on (default: 7000 = 7m, sufficient to enclose all of FASER including FASERnu and ECAL)
+* targetDyMm - the height of the box to generate muons on (default: 600 = 0.6m, sufficient to enclose all of FASER)
+* targetDzMm - the length of the box to generate muons on (default: 600 = 0.6m, sufficient to enclose all of FASER)
+* x0Mm - X offset of the target box from the origin in detector coordinates (default: 0)
+* y0Mm - Y offset of the target box from the origin in detector coordinates (default: 0)
+* z0Mm - Z offset of the target box from the origin in detector coordinates (default: 150 = 15 cm)
* chargeRatio - the mu+:mu- ratio (default: 1.27; energy and angle dependence not yet implemented)
-The Z offset is chosen to make the radius required to enclose all of FASER (including FASERnu and ECAL) as small as possible. This improves the efficiency for making cosmic rays that hit something. For single-detector studies (tracker, ECAL, FASERnu) it will be more efficient to reduce the target radius and shift the center appropriately.
+The Z offset is chosen to make the box required to enclose all of FASER (including FASERnu and ECAL) as small as possible. This improves the efficiency for making cosmic rays that hit something. For single-detector studies (tracker, ECAL, FASERnu) it will be more efficient to reduce the target length and shift the center appropriately.
-The calculated rate of cosmic rays at 85 meter depth is 3.41 Hz for a 1 meter radius target sphere. This scales with the target radius squared, so for the default radius of 3.5 meters, it is 41.8 Hz. Most of these muons will miss the sensitive parts of FASER, however.
+The calculated rate of cosmic rays at 85 meter depth is 3.41 Hz for a 1 meter radius target sphere, and 5.8 Hz for the default 7000 * 600 * 600 mm^3 box.
diff --git a/Generators/FaserCosmicGenerator/python/cosmicSampler.py b/Generators/FaserCosmicGenerator/python/cosmicSampler.py
index 9a7ad6bea925b6b55b37c6129b9cf36be9be47eb..f3dbbaa2e33b299801922914905868dd641c4180 100644
--- a/Generators/FaserCosmicGenerator/python/cosmicSampler.py
+++ b/Generators/FaserCosmicGenerator/python/cosmicSampler.py
@@ -2,7 +2,7 @@
import ROOT,math,random
from random import random
-from math import pi,sin,cos,acos,sqrt
+from math import pi,sin,cos,acos,asin,sqrt
import FaserCosmicGenerator.Range as r
from numpy import array,add
import numpy as np
@@ -16,19 +16,29 @@ from ParticleGun.samplers import Sampler, SampledParticle
class CosmicSampler(Sampler):
"Cosmic ray sampler"
- def __init__(self, maxMuonEnergyGeV = 10000, targetDepthMeters = 85, targetRadiusMm = 3500, thetaMinDegree = 0, thetaMaxDegree = 80, x0Mm = 0, y0Mm = 0, z0Mm = 150, chargeRatio = 1.27):
+ def __init__(self, chargeRatio = 1.27, maxMuonEnergyGeV = 10000, thetaMinDegree = 0, thetaMaxDegree = 80, targetDepthMeters = 85, targetDxMm = 600, targetDyMm = 600, targetDzMm = 7000, x0Mm = 0, y0Mm = 0, z0Mm = 150):
if thetaMinDegree > thetaMaxDegree :
return RuntimeError("thetaMin must be less than or equal to thetaMax")
self.maxEnergy = maxMuonEnergyGeV
self.depth = targetDepthMeters
- self.radiusMm = targetRadiusMm # used for generating start position
- self.radiusCm = self.radiusMm/10 # used for calculating rates
+ # self.radiusMm = targetRadiusMm # used for generating start position
+ # self.radiusCm = self.radiusMm/10 # used for calculating rates
self.cosThMax = math.cos(thetaMinDegree*PI/180)
self.cosThMin = math.cos(thetaMaxDegree*PI/180)
self.x0 = x0Mm
self.y0 = y0Mm
self.z0 = z0Mm
+ self.dx = targetDxMm
+ self.dy = targetDyMm
+ self.dz = targetDzMm
+ self.areaXZ = self.dx * self.dz # top
+ self.areaXY = 2 * self.dx * self.dy # front/back
+ self.areaYZ = 2 * self.dy * self.dz # sides
+ # Accumulate, including constant solid-angle factors
+ self.weightedXZ = PI * self.areaXZ
+ self.weightedXY = self.weightedXZ + 2 * self.areaXY
+ self.weightedYZ = self.weightedXY + 2 * self.areaYZ
self.chargeRatio = chargeRatio
# read range data
# compute rate
@@ -47,11 +57,6 @@ class CosmicSampler(Sampler):
self.genPosition.SetXYZT(x + self.x0, y + self.y0, z + self.z0, 0)
self.genMomentum.SetPxPyPzE(px,py,pz,CR.Efinal)
- # placeholders to test class structure - should create a vertically down 1 TeV mu- starting above the detector
-## self.genPosition.SetY(self.radiusMm)
-## self.genMomentum.SetPxPyPzE(0,-math.sqrt((1.0e6)**2-(105.71)**2),0,1.0e6)
- # end placeholders
-
particles = []
p = SampledParticle(CR.pid)
p.mom = self.genMomentum
@@ -80,22 +85,31 @@ class CosmicRay:
self.costh=0 #cos(th)
self.mc_valid=False #monte carlo/keepit bool
self.deepEnough=False #Does particle reach ddepth?
+ self.top = False
+ self.frontback = False
+ self.sides = False
self.fill()
def fill(self):
#Generate theta/energy combinations until a combo has enough energy to reach detector
while not self.deepEnough or not self.mc_valid:
+ self.genEntry()
self.genCosTh()
Emin = r.muEnergy(self.sampler.depth)
self.genEnergy(Emin, self.sampler.maxEnergy)
self.deepEnough = self.checkDepth()
- if r.keepit(r.getValforMC(self.Einit, self.costh), random()) and self.deepEnough:
+ if self.top :
+ geoWeight = self.costh
+ else:
+ geoWeight = sqrt(1 - self.costh**2)
+
+ if r.keepit(geoWeight * r.getValforMC(self.Einit, self.costh), random()) and self.deepEnough:
self.Efinal = r.muSlant(self.Einit, self.sampler.depth, self.costh)
self.EfMeV = self.Efinal*1000
P = r.getMom(self.EfMeV + self.mass, self.mass)
- self.genCoords( self.sampler.radiusMm )
+ self.genCoords()
self.mom=P*array(self.direction)
if random() < self.sampler.chargeRatio / (self.sampler.chargeRatio + 1):
self.pid = -13 # mu+
@@ -104,7 +118,15 @@ class CosmicRay:
self.mc_valid=True
-
+ def genEntry(self):
+ # Generate entry surface to try
+ r = random()
+ if r < self.sampler.weightedXZ/self.sampler.weightedYZ :
+ self.top = True
+ elif r < self.sampler.weightedXY/self.sampler.weightedYZ :
+ self.frontback = True
+ else :
+ self.sides = True
def genCosTh(self):
cosmin,cosmax=self.sampler.cosThMin,self.sampler.cosThMax
@@ -132,25 +154,27 @@ class CosmicRay:
mudepth=self.setDepth()
return mudepth>=self.sampler.depth
- def genCoords(self, R):
- #ROTATION AROUND X-AXIS TO BE IMPLEMENTED SUCH THAT Z+ -> Y-
- # Generate random polar coords on circle tangent to point of
- # interception at angle theta
- phi = TWOPI*random()
+ def genCoords(self):
+ face = 1
+ if self.top :
+ phi = TWOPI*random()
+ self.pos = ( (random()-0.5)*self.sampler.dx, self.sampler.dy/2, (random()-0.5)*self.sampler.dz )
+ else:
+ if random() < 0.5 :
+ face = -1
+ if self.frontback :
+ phi = asin(2 * random() - 1)
+ self.pos = ( (random()-0.5)*self.sampler.dx, (random()-0.5)*self.sampler.dy, -face * self.sampler.dz/2 )
+ else :
+ phi = acos(1 - 2 * random())
+ self.pos = ( -face * self.sampler.dx/2, (random()-0.5)*self.sampler.dy, (random()-0.5)*self.sampler.dz )
+
+
+
costheta = self.costh
- sintheta = sqrt(1- costheta**2)
+ sintheta = sqrt(1 - costheta**2)
## Do rotation by hand
- # v = array((cos(phi)*sintheta, sin(phi)*sintheta, -costheta))
- v = array((sin(phi)*sintheta, -costheta, cos(phi)*sintheta))
- r = sqrt(random()*R**2)
- alpha = TWOPI*random()
- ihat = np.cross(v, (1,0,0))
- ihat = ihat/np.linalg.norm(ihat)
- jhat = np.cross(ihat, v)
- self.pos = r*cos(alpha)*ihat + r*sin(alpha)*jhat - sqrt(R**2-r**2)*v
- self.direction=v
-
- #Rotation template
- #finalpos=r.dot(r.Rx(PI/2),finalpos/np.linalg.norm(finalpos))
+ self.direction = array((face*sin(phi)*sintheta, -costheta, face*cos(phi)*sintheta))
+
diff --git a/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py b/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py
index 2f3f67e6a7bdf8fd0028ffe99b117a01d0a6feb8..ba402f7a5ec07fcdd627a4bf250f17b53c1c9d63 100644
--- a/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py
+++ b/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py
@@ -98,8 +98,8 @@ def FASEREnvelopeCfg(ConfigFlags, name="Faser", **kwargs):
# kwargs.setdefault("dX", 16.0 * mm)
# kwargs.setdefault("dY", 16.0 * mm)
# kwargs.setdefault("dZ", 33.0 * mm)
- kwargs.setdefault("dX", 4000.0 * mm)
- kwargs.setdefault("dY", 4000.0 * mm)
+ kwargs.setdefault("dX", 600.0 * mm)
+ kwargs.setdefault("dY", 600.0 * mm)
kwargs.setdefault("dZ", 4000.0 * mm)
kwargs.setdefault("SubDetectors", SubDetectorList)