Commit db140255 authored by Andy Buckley's avatar Andy Buckley Committed by Graeme Stewart
Browse files

'Add an InvSampler, and an example JO for flat field curvature and beamline IP...

'Add an InvSampler, and an example JO for flat field curvature and beamline IP sampling. Plus some more general code tweaks.' (ParticleGun-02-01-00)

	* Tagging: ParticleGun-02-01-00

	* Add an example JO showing how to use the inverse sampler to
	sample flat in B-field curvature and how to couple together
	momentum and position sampling to give an impact parameter
	distribution.

2016-04-05  Andy Buckley  <andy.buckley@cern.ch>

	* Explicitly force more sampler constructor args to float type.

	* Add InverseSampler, for e.g. B-field curvature ~ 1/pT.
parent 893e0dd1
#! -*- python -*-
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
include("GeneratorUtils/StdEvgenSetup.py")
theApp.EvtMax = 100
import ParticleGun as PG
class MyParticleSampler(PG.ParticleSampler):
"""
A special sampler to generate single particles flat in 1/pT and in
impact parameter to the beam, with flat z0.
"""
def __init__(self):
psamp = PG.PtEtaMPhiSampler(pt=PG.InvSampler(4000, 400000), eta=[0.1,0.3], phi=[0.3, 0.5])
xsamp = PG.PosSampler(0, 0, [-150,150], 0)
PG.ParticleSampler.__init__(self, pid={13,-13}, mom=psamp, pos=xsamp)
self.ip = PG.mksampler([-2,2])
def shoot(self):
"Return a vector of sampled particles"
ps = PG.ParticleSampler.shoot(self)
assert len(ps) == 1
p = ps[0]
from math import sqrt
m = -p.mom.X() / p.mom.Y() #< gradient of azimuthal IP sampling line, perp to mom
x = self.ip() / sqrt(1 + m**2) #< just decomposing sampled IP into x component...
y = m*x #< ... and y-component
p.pos.SetX(x)
p.pos.SetY(m*x)
return [p]
topSeq += PG.ParticleGun()
topSeq.ParticleGun.randomSeed = 123456
topSeq.ParticleGun.sampler = MyParticleSampler()
include("GeneratorUtils/postJO.CopyWeights.py")
include("GeneratorUtils/postJO.PoolOutput.py")
include("GeneratorUtils/postJO.DumpMC.py")
......@@ -148,8 +148,8 @@ class LogSampler(ContinuousSampler):
"Randomly sample from an exponential distribution (i.e. uniformly on a log scale)."
def __init__(self, low, high):
self.low = low
self.high = high
self.low = float(low)
self.high = float(high)
def shoot(self):
rand = random.random()
......@@ -162,13 +162,25 @@ class GaussianSampler(ContinuousSampler):
"Randomly sample from a 1D Gaussian distribution."
def __init__(self, mean, sigma):
self.mean = mean
self.sigma = sigma
self.mean = float(mean)
self.sigma = float(sigma)
def shoot(self):
return random.gauss(self.mean, self.sigma)
class InvSampler(ContinuousSampler):
"Randomly sample from a 1/x distribution."
def __init__(self, low, high):
self.low = float(low)
self.high = float(high)
def shoot(self):
invx = random.uniform(1/self.high, 1/self.low) #< limit inversion not actually necessary
return 1./invx
########################################
......@@ -242,9 +254,9 @@ def mksampler(x):
- if x can be called, i.e. x() is valid, we just return x;
- a Python list (square brackets) will be converted to a continuous
UniformSampler or DisjointUniformSampler;
- TODO: a Python tuple (round brackets/parentheses) will be treated
- a Python tuple (round brackets/parentheses) will be treated
as a discrete CyclicSeqSampler;
- TODO: a Python set (curly brackets/braces) will be treated
- a Python set (curly brackets/braces) will be treated
as a discrete RandomSeqSampler;
- otherwise a ConstSampler will be created from x, so that x is
returned when the sampler is called.
......@@ -889,6 +901,7 @@ class ParticleSampler(Sampler):
m = self.massdict[abs(pid)]
self.mom.mass = m
p.mass = m
# TODO: Should the particle generated_mass be set from the sampler by default?
## Sample momentum and vertex positions into the particle
p.mom = self.mom()
p.pos = self.pos()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment