diff --git a/Control/CalypsoExample/GenEventExample/CMakeLists.txt b/Control/CalypsoExample/GenEventExample/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..75dc278c2221674e3c42bafeff40470bba512efe
--- /dev/null
+++ b/Control/CalypsoExample/GenEventExample/CMakeLists.txt
@@ -0,0 +1,9 @@
+atlas_subdir( GenEventExample )
+
+atlas_add_component( GenEventExample
+ src/GenEventReadAlg.cxx
+ src/components/GenEventExample_entries.cxx
+ LINK_LIBRARIES AthenaBaseComps GeneratorObjects
+ )
+
+atlas_install_python_modules( python/*.py )
diff --git a/Control/CalypsoExample/GenEventExample/python/GenEventReadExampleConfig.py b/Control/CalypsoExample/GenEventExample/python/GenEventReadExampleConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..3d59aef95f8ab851a0ec8f8985bfb12b397545c4
--- /dev/null
+++ b/Control/CalypsoExample/GenEventExample/python/GenEventReadExampleConfig.py
@@ -0,0 +1,52 @@
+# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+
+#!/usr/bin/env python
+import sys
+from AthenaConfiguration.ComponentFactory import CompFactory
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+
+def GenEventReadExampleCfg(flags, name="GenEventReadExampleAlg", **kwargs):
+
+ a = ComponentAccumulator()
+
+ GenEventReadAlg = CompFactory.GenEventReadAlg
+ a.addEventAlgo(GenEventReadAlg(name, **kwargs))
+
+ thistSvc = CompFactory.THistSvc()
+ thistSvc.Output += ["HIST DATAFILE='genEventHist.root' OPT='RECREATE'"]
+ a.addService(thistSvc)
+ return a
+
+
+if __name__ == "__main__":
+ from AthenaCommon.Configurable import Configurable
+ from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+
+ Configurable.configurableRun3Behavior = True
+
+# Flags for this job
+ ConfigFlags.Input.isMC = True # Needed to bypass autoconfig
+ ConfigFlags.IOVDb.GlobalTag = "OFLCOND-XXXX-XXX-XX" # Needed to bypass autoconfig, only the "OFLCOND" matters at the moment
+ ConfigFlags.GeoModel.FaserVersion = "FASER-01" # Default FASER geometry
+ ConfigFlags.Detector.SimulateFaserSCT = True
+ ConfigFlags.Input.Files = ["my.cosmics.pool.root"]
+ ConfigFlags.lock()
+
+# Configure components
+ from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+ from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+ acc = MainServicesCfg(ConfigFlags)
+ acc.merge(PoolReadCfg(ConfigFlags))
+
+# Set things up to create a conditions DB with neutral Tracker alignment transforms
+ acc.merge(GenEventReadExampleCfg(ConfigFlags))
+
+# Configure verbosity
+ # ConfigFlags.dump()
+ # logging.getLogger('forcomps').setLevel(VERBOSE)
+ # acc.foreach_component("*").OutputLevel = VERBOSE
+ # acc.foreach_component("*ClassID*").OutputLevel = INFO
+ # log.setLevel(VERBOSE)
+
+# Execute and finish
+ sys.exit(int(acc.run(maxEvents=-1).isFailure()))
diff --git a/Control/CalypsoExample/GenEventExample/python/__init__.py b/Control/CalypsoExample/GenEventExample/python/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..bbee1131453aaf11b421ba1b97351dfb8d1ea6fc
--- /dev/null
+++ b/Control/CalypsoExample/GenEventExample/python/__init__.py
@@ -0,0 +1 @@
+# __author__ = 'Ryan Rice-Smith'
diff --git a/Control/CalypsoExample/GenEventExample/src/GenEventReadAlg.cxx b/Control/CalypsoExample/GenEventExample/src/GenEventReadAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2810c5d0e85bca6f566712510c3962e1f7d81bd0
--- /dev/null
+++ b/Control/CalypsoExample/GenEventExample/src/GenEventReadAlg.cxx
@@ -0,0 +1,67 @@
+#include "GenEventReadAlg.h"
+
+
+GenEventReadAlg::GenEventReadAlg(const std::string& name, ISvcLocator* pSvcLocator)
+: AthHistogramAlgorithm(name, pSvcLocator) { m_cosThetaHist = nullptr; m_rVertexHist = nullptr; }
+
+GenEventReadAlg::~GenEventReadAlg() { }
+
+StatusCode GenEventReadAlg::initialize()
+{
+ // initialize a histogram
+ // letter at end of TH1 indicated variable type (D double, F float etc)
+ m_cosThetaHist = new TH1D("CosTheta", "Cosine of ZenithAngle; cos #q; Events/bin", 20, 0., 1.); //first string is root object name, second is histogram title
+ m_rVertexHist = new TH1D("RVertex", "Vertex Radius; r (mm); Events/bin", 40, 0, m_radius * 1.2);
+ m_phiHist = new TH1D("Phi", "Azimuth; phi; Events/bin", 40, -4*atan(1), 4*atan(1));
+ m_rPerpHist = new TH1D("rPerp", "Transverse radius squared; r^2 (mm); Events/bin", 40, 0, pow(m_radius, 2) * 1.2);
+
+ ATH_CHECK(histSvc()->regHist("/HIST/coshist", m_cosThetaHist));
+ ATH_CHECK(histSvc()->regHist("/HIST/rhist", m_rVertexHist));
+ ATH_CHECK(histSvc()->regHist("/HIST/phihist", m_phiHist));
+ ATH_CHECK(histSvc()->regHist("/HIST/rperphist", m_rPerpHist));
+
+ // initialize data handle keys
+ ATH_CHECK( m_mcEventKey.initialize() );
+ ATH_MSG_INFO( "Using GenEvent collection with key " << m_mcEventKey.key());
+ return StatusCode::SUCCESS;
+}
+
+StatusCode GenEventReadAlg::execute()
+{
+ // Handles created from handle keys behave like pointers to the corresponding container
+ SG::ReadHandle<McEventCollection> h_mcEvents(m_mcEventKey);
+ ATH_MSG_INFO("Read McEventContainer with " << h_mcEvents->size() << " events");
+ if (h_mcEvents->size() == 0) return StatusCode::FAILURE;
+ const HepMC::GenEvent* theEvent = h_mcEvents->at(0);
+ if (theEvent->particles_size() != 1 && theEvent->vertices_size() != 1) return StatusCode::FAILURE;
+ // h_mcEvents->at(0)->print( msg().stream() );
+ auto pParticles = theEvent->particles_begin();
+ auto pParticlesEnd = theEvent->particles_end();
+ const HepMC::GenParticle* pHat { nullptr };
+ for (; pParticles != pParticlesEnd ; ++pParticles)
+ {
+ const HepMC::GenParticle* p = *pParticles;
+ if (pHat == nullptr) pHat = p;
+ m_cosThetaHist->Fill(-p->momentum().y()/p->momentum().rho());
+ m_phiHist->Fill(atan2(p->momentum().x(), p->momentum().z()));
+ }
+
+ auto pVertices = theEvent->vertices_begin();
+ auto pVerticesEnd = theEvent->vertices_end();
+ for (; pVertices != pVerticesEnd; ++pVertices)
+ {
+ const HepMC::GenVertex* v = *pVertices;
+ HepMC::FourVector vOffset { v->position().x() , v->position().y(), v->position().z() - m_zOffset, v->position().t() };
+ m_rVertexHist->Fill(vOffset.rho());
+ m_rPerpHist->Fill(pow(vOffset.rho(), 2) - pow((vOffset.x() * pHat->momentum().x() +
+ vOffset.y() * pHat->momentum().y() +
+ vOffset.z() * pHat->momentum().z())/pHat->momentum().rho(), 2));
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode GenEventReadAlg::finalize()
+{
+ return StatusCode::SUCCESS;
+}
\ No newline at end of file
diff --git a/Control/CalypsoExample/GenEventExample/src/GenEventReadAlg.h b/Control/CalypsoExample/GenEventExample/src/GenEventReadAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..b99e07a09e679336b85ab27d464f60cdc18c82d4
--- /dev/null
+++ b/Control/CalypsoExample/GenEventExample/src/GenEventReadAlg.h
@@ -0,0 +1,34 @@
+#include "AthenaBaseComps/AthHistogramAlgorithm.h"
+#include "GeneratorObjects/McEventCollection.h"
+#include <TH1.h>
+#include <math.h>
+// #include <TProfile.h>
+
+/* GenEvent reading example */
+
+class GenEventReadAlg : public AthHistogramAlgorithm
+{
+ public:
+ GenEventReadAlg(const std::string& name, ISvcLocator* pSvcLocator);
+
+ virtual ~GenEventReadAlg();
+
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+ private:
+ TH1* m_cosThetaHist;
+ TH1* m_rVertexHist;
+ TH1* m_phiHist;
+ TH1* m_rPerpHist;
+ // TProfile* m_hprof;
+
+ DoubleProperty m_zOffset { this, "ZOffsetMm", 150.0, "Offset of the generation sphere center" };
+ DoubleProperty m_radius { this, "TargetRadiusMm", 3500.0, "Radius of target sphere in mm" };
+
+ // Read handle keys for data containers
+ // Any other event data can be accessed identically
+ // Note the key names ("GEN_EVENT" or "SCT_Hits") are Gaudi properties and can be configured at run-time
+ SG::ReadHandleKey<McEventCollection> m_mcEventKey { this, "McEventCollection", "GEN_EVENT" };
+};
\ No newline at end of file
diff --git a/Control/CalypsoExample/GenEventExample/src/components/GenEventExample_entries.cxx b/Control/CalypsoExample/GenEventExample/src/components/GenEventExample_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e3a0f06bfd90232e92315168d5b0f1b9f30df3a6
--- /dev/null
+++ b/Control/CalypsoExample/GenEventExample/src/components/GenEventExample_entries.cxx
@@ -0,0 +1,3 @@
+#include "../GenEventReadAlg.h"
+
+DECLARE_COMPONENT( GenEventReadAlg ) //change alg name if update file/class etc
\ No newline at end of file
diff --git a/Generators/FaserCosmicGenerator/README.md b/Generators/FaserCosmicGenerator/README.md
index 4a1e19c73dbd391e77de2b05907fa5778b006165..8ba2dfbc3e2b48000bffb0703aece091684e91de 100644
--- a/Generators/FaserCosmicGenerator/README.md
+++ b/Generators/FaserCosmicGenerator/README.md
@@ -1,3 +1,29 @@
-Faser cosmic-ray generator
+#Faser cosmic-ray generator
-Documentation goes here...
\ No newline at end of file
+## Setup
+
+In Athena 22.0.19, the scipy module is missing from the LCG environment for some reason.
+
+To use the generator, the following command is required after the usual steps to setup the release:
+
+`lsetup "lcgenv -p LCG_98python3_ATLAS_2 x86_64-centos7-gcc8-opt scipy"`
+
+## Usage
+
+The `RunG4Cosmics.py` file shows a default configuration, appropriate for generating cosmic rays in all of FASER.
+
+Options available through the `CosmicSampler` constructor include:
+
+* maxMuonEnergyGeV - the maximum muon energy to sample at the surface (default: 10000 = 10 TeV; setting this higher does not result in appreciably higher rate)
+* 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)
+* 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 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.
diff --git a/Generators/FaserCosmicGenerator/python/Range.py b/Generators/FaserCosmicGenerator/python/Range.py
new file mode 100644
index 0000000000000000000000000000000000000000..1c08928c2be28132f19fea2deddc05ae39e3157d
--- /dev/null
+++ b/Generators/FaserCosmicGenerator/python/Range.py
@@ -0,0 +1,130 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import numpy as np
+
+from math import sqrt,cos,sin
+import scipy.interpolate
+from FaserCosmicGenerator.rockdata import rockdata as rockData
+
+rockDensity = 2.65 # gm/cm^3 for standard rock
+##with open('stdRockData.csv',newline='') as csvFile:
+## reader=csv.DictReader(csvFile)
+## for row in reader:
+## rockData.append(row)
+#print("Read ",len(rockData)," entries")
+# Splines are log(X) vs. log(T), and vice-versa
+# Expect data file has energies in MeV and depth in gm/cm^2
+logE = [np.log(np.float64(d[0])) for d in rockData[1:]]
+logX = [np.log(np.float64(d[3])) for d in rockData[1:]]
+xFromE = scipy.interpolate.CubicSpline(logE, logX, bc_type='natural')
+eFromX = scipy.interpolate.CubicSpline(logX, logE, bc_type='natural')
+
+def muRange(eKin):
+ # returns range in meters of standard rock for specified muon kinetic energy in GeV
+ return np.exp(xFromE(np.log(eKin*1000)))/rockDensity/100
+
+def muEnergy(xRange):
+ # returns muon kinetic energy in GeV to penetrate specified range in meters of standard rock
+ #print('xr:',xRange)
+ return np.exp(eFromX(np.log(xRange*100*rockDensity)))/1000
+
+def muPropagate(eStart, xDepth):
+ # returns the kinetic energy of a muon which starts with kinetic energy eStart after traveling xDepth meters of standard rock
+ fullRange = muRange(eStart)
+ if (fullRange <= xDepth):
+ return 0
+ return muEnergy(fullRange - xDepth)
+
+def muSlant(eStart, xDepth, cosTheta):
+ # returns the kinetic energy of a muon which starts with kinetic energy eStart after penetrating to a depth xDepth in meters of standard rock
+ # cosTheta is the cosine(zenith angle); 1 = vertical downward, 0 = horizontal
+ if (cosTheta <= 0 or cosTheta > 1):
+ print("muSlant: cosTheta (", cosTheta, ") out of range; must be less than or equal to 1 and greater than 0")
+ return 0
+ fullRange = muRange(eStart)
+ slantDepth = xDepth/cosTheta
+ if (fullRange <= slantDepth):
+ return 0
+ return muEnergy(fullRange - slantDepth)
+
+def getCosThetaCorrected(costh,P1=0.102573,P2=-0.068287,P3=0.958633,P4=0.0407253,P5=0.817285):
+ tn=(costh**2)+(P1**2)+(P2*(costh**P3))+(P4*(costh**P5))
+ td=1+(P1**2+P2)+P4
+ return sqrt(tn/td)
+
+def dIdE(costh,Emu):
+ #From DayaBay paper, eqn2
+ #Takes cos(angle of incidence), a muon energy, and returns Flux(?)
+ costhetaS=getCosThetaCorrected(costh)
+
+ B1=(1.1*Emu*costhetaS)/115
+ B2=(1.1*Emu*costhetaS)/850
+ B3=Emu*(costhetaS**1.29)
+
+ A1=Emu*(1+(3.64/B3))
+ A2=(1/(1+B1))+(0.054/(1+B2))
+
+ final=0.14*(A1**-2.7)*A2
+ return final
+
+def keepit(val1,val2,disable=False): #For monte carlo reasons.
+ if disable: return True
+ else: return val1>val2
+
+def getminE(costh): #gfun(x)
+ fdepth=85
+ return muEnergy(fdepth/costh)
+
+def dIdE_swap(Emu,costh): #func(y,x)
+ #From DayaBay paper, eqn2
+ #Takes cos(angle of incidence), a muon energy, and returns Flux(?)
+ costhetaS=getCosThetaCorrected(costh)
+
+ B1=(1.1*Emu*costhetaS)/115
+ B2=(1.1*Emu*costhetaS)/850
+ B3=Emu*(costhetaS**1.29)
+
+ A1=Emu*(1+(3.64/B3))
+ A2=(1/(1+B1))+(0.054/(1+B2))
+
+ final=0.14*(A1**-2.7)*A2
+ return final
+
+def getValforMC(Einit,costh,Edenom=0.09):
+ return (Einit**2.7)*dIdE(costh,Einit)/Edenom
+
+def Rx(th=0):
+ #Rotate on x axis
+ Rx=np.array(((1,0,0),
+ (0,cos(th),-sin(th)),
+ (0,sin(th), cos(th))))
+ return Rx
+
+def Rz(phi=0):
+ #Rotate on x axis
+ Rz = np.array(((cos(phi),-sin(phi),0),
+ (sin(phi), cos(phi),0),
+ (0,0,1)))
+ return Rz
+
+def Ry(psi=0):
+ #Rotate on x axis
+ Ry = np.array(((cos(psi),0,sin(psi)),
+ (0,1,0),
+ (-sin(psi),0,cos(psi))))
+ return Ry
+
+
+def dot(M1,M2):
+ return np.dot(M1,M2)
+
+def getMom(E,m):
+ #assuming c=1
+ return sqrt(E**2-(m**2))
+
+#def getRate(R,cthmin=0,cthmax=1,Emax=100000):
+# a1,a2=dblquad(dIdE_swap,cthmin,cthmax,getminE,Emax)
+# return a1*3.14*(R**2)*2*3.14
+
+
diff --git a/Generators/FaserCosmicGenerator/python/__init__.py b/Generators/FaserCosmicGenerator/python/__init__.py
index 42ea59fd7bebf56ad3497e2c7352f475aaba0950..5258f35af6c9f01055d9fad2188f8a7b6681560c 100644
--- a/Generators/FaserCosmicGenerator/python/__init__.py
+++ b/Generators/FaserCosmicGenerator/python/__init__.py
@@ -1,4 +1,7 @@
# Copyright (C) 2021 CERN for the benefit of the FASER collaboration
from FaserCosmicGenerator.cosmicSampler import CosmicSampler
+import FaserCosmicGenerator.Range
+import FaserCosmicGenerator.rockdata as rockdata
+import numpy,scipy,math, random, ROOT
diff --git a/Generators/FaserCosmicGenerator/python/cosmicSampler.py b/Generators/FaserCosmicGenerator/python/cosmicSampler.py
index 4dca1fed08bdd7321d9bba772677af9af7b5bfa2..9a7ad6bea925b6b55b37c6129b9cf36be9be47eb 100644
--- a/Generators/FaserCosmicGenerator/python/cosmicSampler.py
+++ b/Generators/FaserCosmicGenerator/python/cosmicSampler.py
@@ -1,16 +1,22 @@
# Copyright (C) 2021 CERN for the benefit of the FASER collaboration
import ROOT,math,random
+from random import random
+from math import pi,sin,cos,acos,sqrt
+import FaserCosmicGenerator.Range as r
+from numpy import array,add
+import numpy as np
-PI = math.pi
-TWOPI = 2*math.pi
+
+PI = pi
+TWOPI = 2*pi
from ParticleGun.samplers import Sampler, SampledParticle
class CosmicSampler(Sampler):
"Cosmic ray sampler"
- def __init__(self, maxMuonEnergyGeV = 10000, targetDepthMeters = 85, targetRadiusMm = 3000, thetaMinDegree = 0, thetaMaxDegree = 80):
+ 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")
@@ -20,22 +26,34 @@ class CosmicSampler(Sampler):
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.chargeRatio = chargeRatio
# read range data
# compute rate
- def shoot(self):
+ def shoot(self): #don't, though :(
# generate one event here
self.genPosition = ROOT.TLorentzVector()
self.genMomentum = ROOT.TLorentzVector()
+ self.pdgCode = 13
+
+ CR=CosmicRay(sampler = self)
+
+ 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)
# 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)
- self.pdgCode = 13
+## 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(self.pdgCode)
+ p = SampledParticle(CR.pid)
p.mom = self.genMomentum
p.pos = self.genPosition
particles.append(p)
@@ -44,4 +62,95 @@ class CosmicSampler(Sampler):
def __call__(self):
return self.shoot()
+class CosmicRay:
+ def __init__(self, sampler):
+ self.pid=0
+ self.sampler=sampler #stores settings from Sampler class
+
+ self.mass=105.658 #MeV (muon)
+ self.Einit=0 #GeV
+ self.EiMeV=0 #MeV
+ self.Efinal=0 #GeV
+ self.EfMeV= 0 #MeV
+ 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?
+
+ 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.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.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
+
+
+
+
+ def genCosTh(self):
+ cosmin,cosmax=self.sampler.cosThMin,self.sampler.cosThMax
+ maxdiff=cosmax-cosmin
+ self.costh=cosmin+(maxdiff*random())
+ self.th=acos(self.costh)
+
+ def genEnergy(self, Emin, Emax, gamma=2.7):
+ #Generates initial kinetic energy of muon
+ self.Einit=(Emin**(1-gamma)+random()*(Emax**(1-gamma)-Emin**(1-gamma)))**(1/(1-gamma))
+ self.EiMeV = self.Einit*1000
+
+ def getPtleDepth(self):
+ #Calculates how deep(vertically) a particle with Einit penetrates into earth
+ mulength=r.muRange(self.Einit)
+ return mulength/self.costh
+
+ def setDepth(self):
+ #Calculates how deep(vertically) a particle with Einit penetrates into earth
+ mulength=r.muRange(self.Einit)
+ return mulength*self.costh
+
+ def checkDepth(self):
+ mulength=r.muRange(self.Einit)
+ 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()
+ costheta = self.costh
+ 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))
+
diff --git a/Generators/FaserCosmicGenerator/python/rockdata.py b/Generators/FaserCosmicGenerator/python/rockdata.py
new file mode 100644
index 0000000000000000000000000000000000000000..396116123064a709b7ec42c3385c8603722a5bb4
--- /dev/null
+++ b/Generators/FaserCosmicGenerator/python/rockdata.py
@@ -0,0 +1,146 @@
+rockdata=[['T','a','b','x'],
+ [1,2.643,0.00004771,0.002311],
+ [1.2,35.17,3.98E-05,0.007619],
+[1.4,31.12,3.42E-05,0.01368],
+[1.7,26.67,2.82E-05,0.02413],
+[2,23.42,0.00002408,0.03616],
+[2.5,19.6,0.000019352,0.05961],
+[3,16.94,0.0000162,0.08713],
+[3.5,14.98,1.40E-05,0.1186],
+[4,13.47,1.23E-05,0.1539],
+[4.5,12.27,1.09E-05,0.1928],
+[5,11.29,0.0000099,0.2353],
+[5.5,10.47,0.00000904,0.2814],
+[6,9.784,8.32E-06,0.3308],
+[7,8.681,7.20E-06,0.4395],
+[8,7.834,0.000006355,0.561],
+[9,7.164,5.70E-06,0.6946],
+[10,6.619,0.000005173,0.84],
+[12,5.787,4.39E-06,1.164],
+[14,5.18,3.82E-06,1.53],
+[17,4.524,3.23E-06,2.152],
+[20,4.057,2.81E-06,2.854],
+[25,3.52,2.34E-06,4.183],
+[30,3.157,2.02E-06,5.687],
+[35,2.897,1.80E-06,7.343],
+[40,2.701,1.63E-06,9.133],
+[45,2.55,1.50E-06,11.04],
+[50,2.43,1.40E-06,13.05],
+[55,2.331,1.33E-06,15.15],
+[60,2.249,1.27E-06,17.34],
+[70,2.121,1.18E-06,21.92],
+[80,2.028,1.11E-06,26.75],
+[90,1.958,1.06E-06,31.77],
+[100,1.904,0.00000102,36.95],
+[120,1.828,9.61E-07,47.68],
+[140,1.779,0.00000092,58.78],
+[170,1.734,8.81E-07,75.88],
+[200,1.71,8.56E-07,93.31],
+[250,1.691,8.31E-07,122.7],
+[300,1.688,8.19E-07,152.3],
+[350,1.691,8.13E-07,181.9],
+[400,1.698,8.11E-07,211.4],
+[450,1.707,8.10E-07,240.8],
+[500,1.716,8.12E-07,270],
+[550,1.726,8.14E-07,299.1],
+[600,1.736,8.17E-07,327.9],
+[700,1.756,8.24E-07,385.2],
+[800,1.774,8.32E-07,441.8],
+[900,1.792,8.40E-07,497.9],
+[1000,1.808,8.66E-07,553.4],
+[1200,1.836,9.36E-07,663.1],
+[1400,1.861,9.96E-07,771.2],
+[1700,1.893,1.07E-06,930.9],
+[2000,1.92,0.000001139,1088],
+[2500,1.957,1.24E-06,1346],
+[3000,1.986,1.32E-06,1599],
+[3500,2.011,1.39E-06,1848],
+[4000,2.033,1.45E-06,2095],
+[4500,2.051,1.50E-06,2339],
+[5000,2.067,0.000001549,2581],
+[5500,2.082,1.59E-06,2821],
+[6000,2.095,0.000001631,3059],
+[7000,2.118,1.70E-06,3531],
+[8000,2.138,1.76E-06,3998],
+[9000,2.155,1.82E-06,4461],
+[10000,2.17,0.000001864,4920],
+[12000,2.195,1.95E-06,5827],
+[14000,2.215,2.03E-06,6724],
+[17000,2.241,2.12E-06,8051],
+[20000,2.262,2.20E-06,9360],
+[25000,2.289,0.000002318,11510],
+[30000,2.311,2.41E-06,13620],
+[35000,2.329,2.49E-06,15710],
+[40000,2.344,0.000002565,17760],
+[45000,2.358,2.63E-06,19790],
+[50000,2.369,0.000002682,21800],
+[55000,2.38,2.73E-06,23790],
+[60000,2.389,2.77E-06,25750],
+[70000,2.406,2.85E-06,29630],
+[80000,2.42,2.92E-06,33430],
+[90000,2.433,2.98E-06,37170],
+[100000,2.444,0.000003031,40840],
+[120000,2.463,3.12E-06,48000],
+[140000,2.479,3.19E-06,54950],
+[170000,2.499,3.28E-06,64980],
+[200000,2.515,3.36E-06,74590],
+[250000,2.538,3.45E-06,89770],
+[300000,2.557,3.52E-06,104000],
+[350000,2.573,3.57E-06,117500],
+[400000,2.586,0.000003625,130200],
+[450000,2.598,3.67E-06,142300],
+[500000,2.609,0.000003712,153800],
+[550000,2.619,3.74E-06,164700],
+[600000,2.628,0.00000377,175200],
+[700000,2.644,3.82E-06,194800],
+[800000,2.658,3.86E-06,212900],
+[900000,2.67,0.0000039,229600],
+[1000000,2.681,0.000003934,245300],
+[1200000,2.7,3.98E-06,273700],
+[1400000,2.717,4.02E-06,299000],
+[1700000,2.737,4.06E-06,332400],
+[2000000,2.755,0.000004105,361600],
+[2500000,2.779,0.000004144,403200],
+[3000000,2.798,0.00000418,438400],
+[3500000,2.815,4.21E-06,468800],
+[4000000,2.83,4.23E-06,495700],
+[4500000,2.843,4.25E-06,519700],
+[5000000,2.855,0.000004272,541300],
+[5500000,2.865,4.29E-06,561100],
+[6000000,2.875,4.30E-06,579200],
+[7000000,2.892,4.32E-06,611700],
+[8000000,2.908,0.000004335,640000],
+[9000000,2.921,4.35E-06,665100],
+[10000000,2.933,0.000004365,687700],
+[12000000,2.954,4.38E-06,727000],
+[14000000,2.972,4.40E-06,760300],
+[17000000,2.995,4.42E-06,802500],
+[20000000,3.014,0.000004438,837900],
+[25000000,3.04,0.000004456,886600],
+[30000000,3.062,0.00000447,926400],
+[35000000,3.081,4.48E-06,960100],
+[40000000,3.097,0.000004495,989400],
+[45000000,3.112,4.50E-06,1015000],
+[50000000,3.125,0.000004514,1038000],
+[55000000,3.137,0.00000452,1059000],
+[60000000,3.147,4.53E-06,1078000],
+[70000000,3.167,4.54E-06,1112000],
+[80000000,3.184,4.55E-06,1141000],
+[90000000,3.198,4.56E-06,1166000],
+[100000000,3.212,0.000004563,1189000],
+[120000000,3.235,4.56E-06,1229000],
+[140000000,3.255,4.56E-06,1263000],
+[170000000,3.281,4.56E-06,1305000],
+[200000000,3.302,0.000004563,1340000],
+[250000000,3.332,0.000004564,1389000],
+[300000000,3.356,4.56E-06,1429000],
+[350000000,3.377,4.56E-06,1463000],
+[400000000,3.395,4.56E-06,1492000],
+[450000000,3.411,4.56E-06,1518000],
+[500000000,3.426,0.000004564,1541000],
+[550000000,3.439,4.56E-06,1562000],
+[600000000,3.451,4.56E-06,1581000],
+[700000000,3.473,4.56E-06,1614000],
+[800000000,3.491,4.56E-06,1644000],
+[900000000,3.508,4.56E-06,1669000],
+[1000000000,3.523,0.000004563,1693000]]
diff --git a/Simulation/G4Faser/G4FaserAlg/test/runG4Cosmics.py b/Simulation/G4Faser/G4FaserAlg/test/runG4Cosmics.py
index ad04571192f49a3500e5de4eca60b3a2e12019f5..16153e6b5edc849849bad74fa55b962ae20f57e7 100644
--- a/Simulation/G4Faser/G4FaserAlg/test/runG4Cosmics.py
+++ b/Simulation/G4Faser/G4FaserAlg/test/runG4Cosmics.py
@@ -134,4 +134,4 @@ if __name__ == "__main__":
#
# Execute and finish
#
- sys.exit(int(acc.run(maxEvents=500).isFailure()))
+ sys.exit(int(acc.run(maxEvents=50000).isFailure()))
diff --git a/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py b/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py
index 181a57b6ccc8fe9c16c073852e5c2638ee400b42..2f3f67e6a7bdf8fd0028ffe99b117a01d0a6feb8 100644
--- a/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py
+++ b/Simulation/G4Faser/G4FaserTools/python/G4GeometryToolConfig.py
@@ -98,9 +98,9 @@ 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", 600.0 * mm)
- kwargs.setdefault("dY", 600.0 * mm)
- kwargs.setdefault("dZ", 3000.0 * mm)
+ kwargs.setdefault("dX", 4000.0 * mm)
+ kwargs.setdefault("dY", 4000.0 * mm)
+ kwargs.setdefault("dZ", 4000.0 * mm)
kwargs.setdefault("SubDetectors", SubDetectorList)
diff --git a/faser-common b/faser-common
index 0154f6ecf8286f496911b42665d1a633edd5aca4..a6e556e2b868059d1ef34a5a63075ba1d31c7e6c 160000
--- a/faser-common
+++ b/faser-common
@@ -1 +1 @@
-Subproject commit 0154f6ecf8286f496911b42665d1a633edd5aca4
+Subproject commit a6e556e2b868059d1ef34a5a63075ba1d31c7e6c