In This Package:

TestMasses.py

Go to the documentation of this file.

#!/usr/bin/env python
'''
Test that the kinematics are self consistent

nuwa.py -n 10 -m 'GenMuon.Helpers DYB ADE /dd/Structure/AD/db-ade1' -m GenMuon.TestMasses

'''

from DybPython.DybPythonAlg import DybPythonAlg, SUCCESS, FAILURE
from DybPython.Util import irange

from GaudiKernel import SystemOfUnits as units
from GaudiPython import loaddict

loaddict('libHepMCRflx')        # for HepMC Python wrappers

from math import sqrt, fabs

class MassiveFailure(DybPythonAlg):

    def __init__(self, name = 'MassiveFailure'):
        DybPythonAlg.__init__(self, name)
        self.ppsvc = None
        return

    def initialize(self):
        sc = DybPythonAlg.initialize(self)
        if sc.isFailure(): return sc

        self.ppsvc = self.svc('IParticlePropertySvc','ParticlePropertySvc')
        if not self.ppsvc:
            self.error('Failed to get the ParticlePropertySvc')
            return FAILURE

        return SUCCESS

    def execute(self):
        tes = self.evtSvc()
        gh = tes['/Event/Gen/GenHeader']
        if not gh:
            self.error('Failed to get a GenHeader')
            return FAILURE

        ge = gh.event()
        if not ge:
            self.error('Failed to get a GenEvent')
            return FAILURE

        for part in irange(ge.particles_begin(),ge.particles_end()):
            mom = part.momentum() # HepMC::FourVector
            pid = part.pdg_id()   # int

            prop = self.ppsvc.findByStdHepID(pid)
            if not prop:
                self.warning('No particle property for PID %d'%pid)
                continue

            def check_mass(label,mtest,mcalc):
                dM = fabs(mtest-mcalc)
                dM_rat = dM / mtest
                dM_eps = 0.01*units.MeV
                if dM_rat > 0.001 or dM > dM_eps:
                    self.warning('PID %d has %s mass %.2f MeV but calculated as %.2f MeV' % \
                                 (pid, label, mtest / units.MeV, mcalc / units.MeV))
                    return False
                return True

            mcalc = sqrt(fabs( mom.t()**2 - (mom.x()**2 + mom.y()**2 + mom.z()**2) ))
            ok1 = check_mass('generated',part.generated_mass(),mcalc)
            ok2 = check_mass('property',prop.mass(),mcalc)
            if not ok1 and not ok2:
                error("Failed to get reasonable mass")
                return FAILURE
            else:
                self.info('PID %d has mass %.2f MeV, energy %.2f GeV/c' % \
                              (pid, mcalc / units.MeV, mom.e() / units.GeV ))
            continue
        return SUCCESS


def run(app):
    app.ExtSvc += ['ParticlePropertySvc']

    alg = MassiveFailure()
    app.addAlgorithm(alg)

    

---