From a72ac14efd29f1630d57c6fa309ad774148eaf6a Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sat, 27 Mar 2021 05:12:36 +0000
Subject: [PATCH] Generate on box surface rather than sphere
---
Generators/FaserCosmicGenerator/README.md | 14 ++--
.../python/cosmicSampler.py | 84 ++++++++++++-------
.../python/G4GeometryToolConfig.py | 4 +-
3 files changed, 64 insertions(+), 38 deletions(-)
diff --git a/Generators/FaserCosmicGenerator/README.md b/Generators/FaserCosmicGenerator/README.md
index 8ba2dfbc..3490a120 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 9a7ad6be..f3dbbaa2 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 2f3f67e6..ba402f7a 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)
--
GitLab