From b9eeb881c8a2ce1064983b7b70ae90b3c80c063d Mon Sep 17 00:00:00 2001 From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk> Date: Wed, 20 Apr 2022 14:22:25 +0100 Subject: [PATCH] Small fixes to generators + algo to shift LOS --- .../python/FaserParticleGunConfig.py | 16 +- .../python/RadialPosSampler.py | 6 +- .../FlukaReader/python/FlukaReaderAlg.py | 69 ++++---- .../share/generate_forsee_events.py | 26 +-- Generators/GeneratorUtils/CMakeLists.txt | 10 ++ Generators/GeneratorUtils/python/ShiftLOS.py | 165 ++++++++++++++++++ Generators/ParticleGun/python/samplers.py | 9 +- 7 files changed, 247 insertions(+), 54 deletions(-) create mode 100644 Generators/GeneratorUtils/CMakeLists.txt create mode 100644 Generators/GeneratorUtils/python/ShiftLOS.py diff --git a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py index 63571637..866a51be 100644 --- a/Generators/FaserParticleGun/python/FaserParticleGunConfig.py +++ b/Generators/FaserParticleGun/python/FaserParticleGunConfig.py @@ -44,10 +44,18 @@ def FaserParticleGunSingleParticleCfg(ConfigFlags, **kwargs) : theta = kwargs.setdefault("theta", [0, pi/20]), phi = kwargs.setdefault("phi", [0, 2*pi]), mass = kwargs.setdefault("mass", 0.0) ) - pg.sampler.pos = PG.PosSampler(x = kwargs.setdefault("x", [-5, 5]), - y = kwargs.setdefault("y", [-5, 5]), - z = kwargs.setdefault("z", -3750.0), - t = kwargs.setdefault("t", 0.0) ) + + if "radius" in kwargs: + pg.sampler.pos = RadialPosSampler(x = kwargs.setdefault("x", 0.0), + y = kwargs.setdefault("y", 0.0), + z = kwargs.setdefault("z", -3750.0), + r = kwargs.setdefault("radius", 1.0), + t = kwargs.setdefault("t", 0.0) ) + else: + pg.sampler.pos = PG.PosSampler(x = kwargs.setdefault("x", [-5, 5]), + y = kwargs.setdefault("y", [-5, 5]), + z = kwargs.setdefault("z", -3750.0), + t = kwargs.setdefault("t", 0.0) ) return cfg diff --git a/Generators/FaserParticleGun/python/RadialPosSampler.py b/Generators/FaserParticleGun/python/RadialPosSampler.py index 4a14987a..a0ae101c 100644 --- a/Generators/FaserParticleGun/python/RadialPosSampler.py +++ b/Generators/FaserParticleGun/python/RadialPosSampler.py @@ -37,8 +37,10 @@ class RadialPosSampler(Sampler): fwhm = 2*self.radius sig = fwhm/(2 * sqrt(2 * log(2))) - # return random.uniform(0, self.radius) - return random.gauss(0, self.radius) + if self.radius < 0: + return random.uniform(0, abs(self.radius)) + else: + return random.gauss(0, self.radius) # return random.gauss(0, sig) @property diff --git a/Generators/FlukaReader/python/FlukaReaderAlg.py b/Generators/FlukaReader/python/FlukaReaderAlg.py index fb49faba..b54aa07a 100644 --- a/Generators/FlukaReader/python/FlukaReaderAlg.py +++ b/Generators/FlukaReader/python/FlukaReaderAlg.py @@ -14,11 +14,12 @@ import math # TODO: correct angle for beam crossing angle: both in angel itself and position class FlukaReader(EvgenAlg): - def __init__(self, name="FlukaReader", MCEventKey="BeamTruthEvent", file_name = "", dist = 0, randomSeed = None, nsamples = 1, test = True): + def __init__(self, name="FlukaReader", MCEventKey="BeamTruthEvent", file_name = "", dist = 0, z = -1.5, randomSeed = None, nsamples = 1, test = False): super(FlukaReader,self).__init__(name=name) self.McEventKey = MCEventKey self.file_name = file_name - self.dist = dist * 100 + self.dist = dist * 100 # cm + self.z = z * 1000 # mm self.isample = 0 self.nsamples = nsamples self.test = test @@ -84,7 +85,6 @@ class FlukaReader(EvgenAlg): "Convert cos(theta) wrt x or y axis to theta wrt z axis" return np.pi/2 - np.arccos(cosTheta) - def pid(self, ftype): "Convert fluka particle type to PID" if ftype == 10: # mu+ @@ -108,7 +108,7 @@ class FlukaReader(EvgenAlg): def energy_after_loss_exact(self, e, zcorr): "Calculate exact energy after loss in material" - return Range.muPropagate(e, zcorr/100.) + return Range.muPropagate(e, zcorr/100.) # meters def energy_after_loss(self, e, cosThetaX, cosThetaY, zcorr, a = 2e-3, b = 4e-6): "Calculate approximate energy after loss in material" @@ -196,14 +196,15 @@ class FlukaReader(EvgenAlg): xout = self.scattering_postition(x, cosX, cosY, z, theta0, rand1, rand2) yout = self.scattering_postition(y, cosY, cosX, z, theta0, rand1, rand2) - # Update entry info - entry["E"] = eout - entry["x"] = xout - entry["y"] = yout - entry["cosX"] = np.cos(np.pi/2 + thetaXout) - entry["cosY"] = np.cos(np.pi/2 + thetaYout) + # Update entry info using copy for cases when resample so don't change the original + newentry = entry.copy() + newentry["E"] = eout + newentry["x"] = xout + newentry["y"] = yout + newentry["cosX"] = np.cos(np.pi/2 + thetaXout) + newentry["cosY"] = np.cos(np.pi/2 + thetaYout) - return entry + return newentry def process(self, entry, evt): @@ -217,21 +218,25 @@ class FlukaReader(EvgenAlg): if self.dist != 0: # Propoagate to FASER - entry = self.propagate(entry) + newentry = self.propagate(entry) elif self.nsamples != 1: - # Else smear if sampling more than once - entry["E"] *= self.rng.normal(1, 0.05) - entry["x"] *= self.rng.normal(1, 0.05) - entry["y"] *= self.rng.normal(1, 0.05) - entry["cosX"] = np.cos(np.arccos(entry["cosX"]) * self.rng.normal(1, 0.05)) - entry["cosY"] = np.cos(np.arccos(entry["cosY"]) * self.rng.normal(1, 0.05)) + # Else smear if sampling more than once, using copy to avoid changing original + newentry = entry.copy() + newentry["E"] *= self.rng.normal(1, 0.05) + newentry["x"] *= self.rng.normal(1, 0.05) + newentry["y"] *= self.rng.normal(1, 0.05) + newentry["cosX"] = np.cos(np.arccos(entry["cosX"]) * self.rng.normal(1, 0.05)) + newentry["cosY"] = np.cos(np.arccos(entry["cosY"]) * self.rng.normal(1, 0.05)) + else: + # No propagation or smearing + newentry = entry if self.msg.level > DEBUG: - print("Propagated/Smeared Entry", entry) + print("Propagated/Smeared Entry", newentry) if self.test: - for k,v in entry.items(): + for k,v in newentry.items(): self.after[k].append(float(v)) try: @@ -240,14 +245,14 @@ class FlukaReader(EvgenAlg): from AthenaPython.PyAthena import HepMC as HepMC # Add weight, correcting for mutliple sampling - evt.weights().push_back(entry["w"] / self.nsamples) + evt.weights().push_back(newentry["w"] / self.nsamples) # Setup MC event mcEventType = EventType() mcEventType.add_type(EventType.IS_SIMULATION) - mcEventId = EventID(run_number = entry["run"], event_number = entry["event"]) + mcEventId = EventID(run_number = newentry["run"], event_number = newentry["event"]) mcEventInfo = EventInfo(id = mcEventId, type = mcEventType) @@ -258,8 +263,7 @@ class FlukaReader(EvgenAlg): ROOT.SetOwnership(mcEventInfo, False) # Create HepMC Vertex - z = 0 - pos = HepMC.FourVector(entry["x"] * cm, entry["y"] * cm, z, 0) + pos = HepMC.FourVector(newentry["x"] * cm, newentry["y"] * cm, self.z, 0) gv = HepMC.GenVertex(pos) ROOT.SetOwnership(gv, False) @@ -271,17 +275,17 @@ class FlukaReader(EvgenAlg): gp = HepMC.GenParticle() m = 105.66 - e = entry["E"] * 1000. + e = newentry["E"] * 1000. #MeV p = np.sqrt(e**2 - m**2) - thetaX = self.angle(entry["cosX"]) - thetaY = self.angle(entry["cosY"]) + thetaX = self.angle(newentry["cosX"]) + thetaY = self.angle(newentry["cosY"]) # theta: just above z axis as phi deals with negative theta = np.abs(thetaY) # phi: 0 - 2pi - phi = np.arctan2(entry["cosY"], entry["cosX"]) - #phi = np.arctan(entry["cosY"] / entry["cosX"]) + phi = np.arctan2(newentry["cosY"], newentry["cosX"]) + #phi = np.arctan(newentry["cosY"] / newentry["cosX"]) if phi < 0: phi += 2*np.pi if phi == 2*np.pi: phi = 0 @@ -293,7 +297,7 @@ class FlukaReader(EvgenAlg): gp.set_momentum(mom) gp.set_generated_mass(m) - gp.set_pdg_id(self.pid(entry["type"])) + gp.set_pdg_id(self.pid(newentry["type"])) gp.set_status(1) ROOT.SetOwnership(gp, False) @@ -390,7 +394,8 @@ if __name__ == "__main__": import argparse, sys parser = argparse.ArgumentParser(description="Run Fluka reader") parser.add_argument("file", help = "Path to fluka file") - parser.add_argument("--dist", "-z", default = 0, type = float, help = "depth of standard rock to propagate through [m]") + parser.add_argument("--dist", "-d", default = 0, type = float, help = "depth of standard rock to propagate through [m]") + parser.add_argument("--pos", "-z", default = -3.75, type = float, help = "Position in z in FASER coordinate system [m]") parser.add_argument("--output", "-o", default = "evgen.pool.root", help = "Name of output file") parser.add_argument("--mcEventKey", "-k", default = "BeamTruthEvent", help = "Name of MC collection") parser.add_argument("--nevents", "-n", default = -1, type = int, help = "Number of events to process") @@ -421,7 +426,7 @@ if __name__ == "__main__": from AthenaConfiguration.ComponentFactory import CompFactory acc = ComponentAccumulator() - reader = FlukaReader("FlukaReader", MCEventKey=args.mcEventKey, file_name = args.file, dist = args.dist, randomSeed = args.randomSeed, nsamples = args.nsamples, test = args.test) + reader = FlukaReader("FlukaReader", MCEventKey=args.mcEventKey, file_name = args.file, dist = args.dist, z = args.pos, randomSeed = args.randomSeed, nsamples = args.nsamples, test = args.test) reader.OutputLevel = INFO acc.addEventAlgo(reader) cfg.merge(acc) diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py index bd44eb2f..37372b7b 100644 --- a/Generators/ForeseeGenerator/share/generate_forsee_events.py +++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py @@ -9,7 +9,7 @@ class ForeseeGenerator(object): Generate LLP particles within FASER acceptance from FORESEE """ - def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, outdir = None, path = '.'): + def __init__(self, modelname, energy, mass, couplings, daughter1_pid, daughter2_pid, outdir = None, path = '.', randomSeed = 12345): self.modelname = modelname self.energy = energy @@ -20,9 +20,7 @@ class ForeseeGenerator(object): self.outdir = outdir self.path = path self.version = 1 # Forsee "version": 2 is in testing - - #if not os.path.exists("files"): - # os.symlink(os.path.expandvars("$Calypso_DIR/../calypso/Generators/foresee/files"), "files") + self.seed = randomSeed # Set decay mode ... @@ -43,7 +41,9 @@ class ForeseeGenerator(object): self.foresee = Foresee() else: self.foresee = Foresee(path = self.path) - self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1", + + # TODO: relax this a bit as daughters may enter even if mother doesn't + self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1", channels=[self.mode], distance=480, length=1.5 , luminosity=1/1000.) # 1 pb-1 @@ -289,14 +289,13 @@ class ForeseeGenerator(object): def write_hepmc(self, nevents): if self.outdir is None: - self.outdir = f"{self.foresee.dirpath}/Models/{self.modelname}/model/events" - - if not os.path.exists(self.outdir): + self.outdir = "model/events/" + elif not os.path.exists(self.outdir): os.mkdir(self.outdir) - filename = "test.hepmc" #f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.hepmc" + filename = f"{self.outdir}/events_{self.energy}TeV_m{self.mass}GeV_c{self.couplings[0]}to_{self.daughter1_pid}_{self.daughter2_pid}.hepmc" - self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents) + self.foresee.write_events(self.mass, self.couplings[0], self.energy, filename, nevents, zfront = -1500, seed = self.seed) def setup_foresee(path): @@ -361,10 +360,11 @@ if __name__ == "__main__": parser.add_argument("--pid1", required = True, type = int, help = "PID of daughter 1") parser.add_argument("--pid2", default = None, type = int, help = "PID of daughter 2 (if not set then will be -PID1)") parser.add_argument("--Ecom", default = "14", help = "Center of mass energy [TeV]") - parser.add_argument("--outdir", "-o", default = ".", help = "Output path") + parser.add_argument("--outdir", "-o", default = None, help = "Output path") parser.add_argument("--path", default = ".", help = "Path to foresee installation") parser.add_argument("--hepmc", action = "store_true", help = "Write HepMC events") - parser.add_argument("--nevents", "-n", default = 10, type = int, help = "Number of HepMC events ") + parser.add_argument("--nevents", "-n", default = 10, type = int, help = "Number of HepMC events ") + parser.add_argument("--randomSeed", "-s", default = 1234, type = int, help = "Random seed for HepMC generation") args = parser.parse_args() add_to_python_path(f"{args.path}/src") @@ -382,7 +382,7 @@ if __name__ == "__main__": print(f" decay = {args.pid1} {args.pid2}") print(f" couplings = {couplings}") - f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2, outdir = args.outdir, path = args.path) + f = ForeseeGenerator(args.model, args.Ecom, args.mass, couplings, args.pid1, args.pid2, outdir = args.outdir, path = args.path, randomSeed = args.randomSeed) if args.hepmc: f.write_hepmc(args.nevents) diff --git a/Generators/GeneratorUtils/CMakeLists.txt b/Generators/GeneratorUtils/CMakeLists.txt new file mode 100644 index 00000000..d471fbe1 --- /dev/null +++ b/Generators/GeneratorUtils/CMakeLists.txt @@ -0,0 +1,10 @@ +################################################################################ +# Package: GeneratorUtils +################################################################################ + +# Declare the package name: +atlas_subdir( GeneratorUtils ) + +# Install files from the package: +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) + diff --git a/Generators/GeneratorUtils/python/ShiftLOS.py b/Generators/GeneratorUtils/python/ShiftLOS.py new file mode 100644 index 00000000..da6eff7e --- /dev/null +++ b/Generators/GeneratorUtils/python/ShiftLOS.py @@ -0,0 +1,165 @@ + +from AthenaCommon.AppMgr import ServiceMgr as svcMgr +from AthenaPython import PyAthena +from AthenaPython.PyAthena import StatusCode, McEventCollection, CLHEP +from AthenaCommon.SystemOfUnits import m +import ROOT + +try: + from AthenaPython.PyAthena import HepMC3 as HepMC +except ImportError: + from AthenaPython.PyAthena import HepMC as HepMC + +class ShiftLOS(PyAthena.Alg): + def __init__(self, name="ShiftLOS", InputMCEventKey="BeamTruthEvent", OutputMCEventKey="BeamTruthEventShifted", xcross = 0.0, ycross = 0.0, simple = False): + super(ShiftLOS,self).__init__(name=name) + self.InputMCEventKey = InputMCEventKey + self.OutputMCEventKey = OutputMCEventKey + self.xcross = xcross * 1e-6 + self.ycross = ycross * 1e-6 + self.distance = 480*m # Assumes 480m is 0 of FASER coordinate system + return + + def shift_vertices(self, evt): + + # Don't need to shift if at IP + if not self.distance: + return evt + + # Loop over all vertices + for v in evt.vertices: + # Get position + pos = v.position() + x = pos.x() + y = pos.y() + z = pos.z() + dz = self.distance + z + + # Shift x or y by appropriate crossing angle + if self.xcross: + x += dz * self.xcross + self.msg.info(f"Shifting x by {self.xcross} over {dz}: {pos.x()} -> {x} ") + elif self.ycross: + y += dz * self.ycross + self.msg.info(f"Shifting y by {self.ycross} over {dz}: {pos.y()} -> {y} ") + else: + self.msg.error(f"Must suply crossing angle in either x or y (but not both)") + return StatusCode.Failure + + v.set_position(HepMC.FourVector(x, y, z, pos.t())) + + return evt + + + def boost_particles(self, evt): + + pxsum, pysum = 0,0 + pxsum_orig, pysum_orig = 0,0 + + # Loop over all particles + for p in evt.particles: + # Get momentum + mom = p.momentum() + + pxsum_orig += mom.x() + pysum_orig += mom.y() + + # Boost in x or y using CLHEP + boost = CLHEP.Hep3Vector(self.xcross, self.ycross, 0.0) + tmp = CLHEP.HepLorentzVector(mom.px(), mom.py(), mom.pz(), mom.e()) + tmp.boost(boost) + + pxsum += tmp.x() - mom.x() + pysum += tmp.y() - mom.y() + + p.print() + # Convert back to HepMC + p.set_momentum(HepMC.FourVector(tmp.px(), tmp.py(), tmp.pz(), tmp.e())) + p.print() + + self.msg.debug(f"Change in total px = {pxsum:.1f} MeV ({pxsum/pxsum_orig * 100: .3f} %), change in total py = {pysum:.1f} MeV ({pysum/pysum_orig * 100: .3f} %)") + + return evt + + def execute(self): + self.msg.debug(f"Exectuing {self.getName()}") + + self.msg.debug(f"Reading {self.InputMCEventKey}") + inevt = self.evtStore[self.InputMCEventKey][0] + + self.msg.debug("Creating output event and collection") + outcoll = McEventCollection() + ROOT.SetOwnership(outcoll, False) + + # Clone input event + outevt = HepMC.GenEvent(inevt.__follow__()) # go from ElementProxy to element itself + + # Modify + outevt = self.shift_vertices(outevt) + outevt = self.boost_particles(outevt) + + # Write output + outcoll.push_back(outevt) + ROOT.SetOwnership(outevt, False) + + self.msg.debug(f"Recording {self.OutputMCEventKey}") + self.evtStore.record(outcoll, self.OutputMCEventKey, True, False) + + return StatusCode.Success + +if __name__ == "__main__": + import argparse, sys + parser = argparse.ArgumentParser(description="Run ShiftLOS") + parser.add_argument("infile", help = "Path to input EVNT file") + parser.add_argument("outfile", help = "Path to output EVNT file") + parser.add_argument("--InputMCEventKey", "-i", default = "BeamTruthEvent", help = "Name of Input MC collection") + parser.add_argument("--OutputMCEventKey", "-o", default = "BeamTruthEventShifted", help = "Name of Output MC collection") + parser.add_argument("--xcross", "-x", default = 0, type = float, help = "Crossing angle of LHC beam in x [urad]") + parser.add_argument("--ycross", "-y", default = 0, type = float, help = "Crossing angle of LHC beam in y [urad]") + parser.add_argument("--nevents", "-n", default = -1, type = int, help = "Number of events to process") + args = parser.parse_args() + + from AthenaCommon.Logging import log + from AthenaCommon.Constants import DEBUG, INFO + + from AthenaCommon.Configurable import Configurable + Configurable.configurableRun3Behavior = 1 + + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + ConfigFlags.Input.isMC = True + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" # Always needed; must match FaserVersion + ConfigFlags.GeoModel.FaserVersion = "FASER-01" # Default FASER geometry + ConfigFlags.Detector.EnableFaserSCT = True + ConfigFlags.Input.Files = [ args.infile ] + ConfigFlags.Output.EVNTFileName = args.outfile + ConfigFlags.lock() + + # Configure components + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + cfg = MainServicesCfg(ConfigFlags) + + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + + cfg = MainServicesCfg(ConfigFlags) + cfg.merge(PoolReadCfg(ConfigFlags)) + + import McParticleEvent.Pythonizations + + from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator + from AthenaConfiguration.ComponentFactory import CompFactory + + acc = ComponentAccumulator() + alg = ShiftLOS("ShiftLOS", InputMCEventKey=args.InputMCEventKey, OutputMCEventKey=args.OutputMCEventKey, xcross = args.xcross, ycross = args.ycross) + alg.OutputLevel = INFO + acc.addEventAlgo(alg) + cfg.merge(acc) + + itemList = [ "EventInfo#McEventInfo", "McEventCollection#*" ] + from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + cfg.merge(OutputStreamCfg(ConfigFlags, "EVNT", itemList, disableEventTag = True)) + + sc = cfg.run(maxEvents = args.nevents) + sys.exit(not sc.isSuccess()) + + diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py index 4a232f7d..90e9676a 100644 --- a/Generators/ParticleGun/python/samplers.py +++ b/Generators/ParticleGun/python/samplers.py @@ -160,13 +160,16 @@ class LogSampler(ContinuousSampler): class GaussianSampler(ContinuousSampler): "Randomly sample from a 1D Gaussian distribution." - def __init__(self, mean, sigma): + def __init__(self, mean, sigma, oneside = False): self.mean = float(mean) self.sigma = float(sigma) + self.oneside = bool(oneside) def shoot(self): - return random.gauss(self.mean, self.sigma) - + if self.oneside: + return abs(random.gauss(self.mean, self.sigma)) + else: + return random.gauss(self.mean, self.sigma) class InvSampler(ContinuousSampler): "Randomly sample from a 1/x distribution." -- GitLab