From db1402557b27ac686b9c42d05f6c69cc7250d7c6 Mon Sep 17 00:00:00 2001
From: Andy Buckley <a.g.buckley@gmail.com>
Date: Thu, 14 Apr 2016 21:14:39 +0200
Subject: [PATCH] '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.
---
 ...Option.ParticleGun_flatcurvature_flatip.py | 41 +++++++++++++++++++
 Generators/ParticleGun/python/samplers.py     | 25 ++++++++---
 2 files changed, 60 insertions(+), 6 deletions(-)
 create mode 100644 Generators/ParticleGun/examples/jobOption.ParticleGun_flatcurvature_flatip.py

diff --git a/Generators/ParticleGun/examples/jobOption.ParticleGun_flatcurvature_flatip.py b/Generators/ParticleGun/examples/jobOption.ParticleGun_flatcurvature_flatip.py
new file mode 100644
index 00000000000..38233563de9
--- /dev/null
+++ b/Generators/ParticleGun/examples/jobOption.ParticleGun_flatcurvature_flatip.py
@@ -0,0 +1,41 @@
+#! -*- 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")
diff --git a/Generators/ParticleGun/python/samplers.py b/Generators/ParticleGun/python/samplers.py
index b76af6a8889..423f07ac124 100644
--- a/Generators/ParticleGun/python/samplers.py
+++ b/Generators/ParticleGun/python/samplers.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()
-- 
GitLab