Add cosmic generator

Generators/FaserCosmicGenerator/python/cosmicSampler.py

@@ -16,7 +16,7 @@ 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):
+    def __init__(self, maxMuonEnergyGeV = 10000, targetDepthMeters = 85, targetRadiusMm = 3500, thetaMinDegree = 0, thetaMaxDegree = 80, x0Mm = 0, y0Mm = 0, z0Mm = 150, chargeRatio = 1.27):
         if thetaMinDegree > thetaMaxDegree :
             return RuntimeError("thetaMin must be less than or equal to thetaMax")
 
@@ -29,6 +29,7 @@ class CosmicSampler(Sampler):
         self.x0 = x0Mm
         self.y0 = y0Mm
         self.z0 = z0Mm
+        self.chargeRatio = chargeRatio
         # read range data
         # compute rate
 
@@ -38,11 +39,10 @@ class CosmicSampler(Sampler):
         self.genMomentum = ROOT.TLorentzVector()
         self.pdgCode = 13
 
-        CR=CosmicRay(Emax = self.maxEnergy, ddepth = self.depth, pid = self.pdgCode, sampler=self)
-        CR.fill(R = self.radiusMm)
+        CR=CosmicRay(sampler = self)
 
-        x,y,z = CR.pos
-        px,py,pz=CR.mom
+        x, y, z    = CR.pos
+        px, py, pz = CR.mom
 
         self.genPosition.SetXYZT(x + self.x0, y + self.y0, z + self.z0, 0)
         self.genMomentum.SetPxPyPzE(px,py,pz,CR.Efinal)
@@ -53,7 +53,7 @@ class CosmicSampler(Sampler):
         # end placeholders
 
         particles = []
-        p = SampledParticle(self.pdgCode)
+        p = SampledParticle(CR.pid)
         p.mom = self.genMomentum
         p.pos = self.genPosition
         particles.append(p)
@@ -63,9 +63,8 @@ class CosmicSampler(Sampler):
         return self.shoot()
 
 class CosmicRay:
-    def __init__(self, Emax, pid, ddepth, sampler):
-        self.pid=pid
-        self.ddepth=ddepth #detector depth in m
+    def __init__(self, sampler):
+        self.pid=0
         self.sampler=sampler #stores settings from Sampler class
 
         self.mass=105.658 #MeV (muon)
@@ -73,31 +72,35 @@ class CosmicRay:
         self.EiMeV=0 #MeV
         self.Efinal=0 #GeV
         self.EfMeV= 0 #MeV
-        self.Emax = Emax #GeV
         self.pos=(0,0,0)
         self.direction=(0,0,0) #
         self.mom=(0,0,0)
         self.th=0 #theta, radians
         self.phi=0 #angle that particle approaches z (vertical) axis
         self.costh=0 #cos(th)
-
         self.mc_valid=False #monte carlo/keepit bool
         self.deepEnough=False #Does particle reach ddepth?
 
-    def fill(self, R, Edenom=0.09):
+        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.genCosTh()
-            Emin = r.muEnergy(self.ddepth)
-            self.genEnergy(Emin, self.Emax)
+            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:
-                self.Efinal=r.muSlant(self.Einit,self.ddepth,self.costh)
-                self.EfMeV=self.Efinal*1000
-                P=r.getMom(self.EfMeV + self.mass, self.mass)
-                self.genCoords( R )
+            if r.keepit(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.mom=P*array(self.direction)
+                if random() < self.sampler.chargeRatio / (self.sampler.chargeRatio + 1):
+                    self.pid = -13  # mu+
+                else:
+                    self.pid = 13   # mu-
                 self.mc_valid=True
         
             
@@ -127,7 +130,7 @@ class CosmicRay:
     def checkDepth(self):
         mulength=r.muRange(self.Einit)
         mudepth=self.setDepth()
-        return mudepth>=self.ddepth
+        return mudepth>=self.sampler.depth
 
     def genCoords(self, R):
         #ROTATION AROUND X-AXIS TO BE IMPLEMENTED SUCH THAT Z+ -> Y-