Skip to content
Snippets Groups Projects
rich_data_gen.py 2.91 KiB
Newer Older
  • Learn to ignore specific revisions
  • import os
    from pyHepMC3 import HepMC3 as hm
    import numpy as np
    import argparse
    
    
    PARTICLES = [
    #    (111, 0.134.9766),      # pi0
        (211, 0.13957018),      # pi+
    #    (-211, 0.13957018),     # pi-
    #    (311, 0.497648),        # K0
        (321, 0.493677),        # K+
    #    (-321, 0.493677),       # K-
        (2212, 0.938272),       # proton
    #    (2112, 0.939565),       # neutron
    #    (11, 0.51099895e-3),    # electron
    #    (-11, 0.51099895e-3),   # positron
    ]
    
    
    # p in GeV, angle in degree, vertex in mm
    def gen_event(prange=(8, 100), arange=(0, 20)):
        evt = hm.GenEvent(momentum_unit=hm.Units.MomentumUnit.GEV, length_unit=hm.Units.LengthUnit.MM)
        pid, mass = PARTICLES[np.random.randint(len(PARTICLES))]
    
        # final state
        state = 1
    
        # momentum, angles, energy
        p = np.random.uniform(*prange)
        theta = np.random.uniform(*arange)*np.pi/180.
        phi = np.random.uniform(0., 2.*np.pi)
        e0 = np.sqrt(p*p + mass*mass)
    
        px = np.cos(phi)*np.sin(theta)
        py = np.sin(phi)*np.sin(theta)
        pz = np.cos(theta)
    
        # beam
        pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
        ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), -11, 4)
    
        hout = hm.GenParticle(hm.FourVector(px*p, py*p, pz*p, e0), pid, state)
        # evt.add_particle(part)
        vert = hm.GenVertex()
        vert.add_particle_in(ebeam)
        vert.add_particle_in(pbeam)
        vert.add_particle_out(hout)
        evt.add_vertex(vert)
        return evt
    
    
    if __name__ == "__main__":
        parser = argparse.ArgumentParser('RICH dataset generator')
    
        parser.add_argument('output', help='path to the output file')
        parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
        parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
        parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
        parser.add_argument('--angmin', type=float, default=0.0, dest='angmin', help='minimum angle in degree')
        parser.add_argument('--angmax', type=float, default=20.0, dest='angmax', help='maximum angle in degree')
        parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
    
        args = parser.parse_args()
    
        # random seed
        if args.seed < 0:
            args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
        np.random.seed(args.seed)
    
        output = hm.WriterAscii(args.output);
        if output.failed():
            print("Cannot open file \"{}\"".format(args.output))
            sys.exit(2)
    
        count = 0
        while count < args.nev:
            if (count % 1000 == 0):
                print("Generated {} events".format(count), end='\r')
            evt = gen_event((args.pmin, args.pmax), (args.angmin, args.angmax))
            output.write_event(evt)
            evt.clear()
            count += 1
    
        print("Generated {} events".format(args.nev))
        output.close()