Generate on box surface rather than sphere

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.

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))
+
 
     

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)