pythia8darkphoton_conf.py

Go to the documentation of this file.

# SPDX-License-Identifier: LGPL-3.0-or-later
# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
 
import contextlib
import os
 
import darkphoton
import proton_bremsstrahlung
import readDecayTable
import ROOT
import shipunit as u
from method_logger import MethodLogger
 
# Boundaries for production in meson decays
# mass of the mesons
pi0mass = 0.1349770
etamass = 0.547862
omegamass = 0.78265
eta1mass = 0.95778
 
 
def addDPtoROOT(pid=9900015, m=0.2, g=4.866182e-04):
    pdg = ROOT.TDatabasePDG.Instance()
    pdg.AddParticle("A", "DarkPhoton", m, False, g, 0.0, "A", pid)
 
 
def readFromAscii():
    FairShip = os.environ["FAIRSHIP"]
    with open(FairShip + "/shipgen/branchingratios.dat") as ascii:
        content = ascii.readlines()
    h = {}
    n = 0
    while n < len(content):
        line = content[n]
        if "TH1F" in line:
            keys = line.split("|")
            n += 1
            limits = content[n].split(",")
            hname = keys[1]
            if len(keys) < 5:
                keys.append(",")
            h[hname] = ROOT.TH1F(
                hname, keys[2] + ";" + keys[3] + ";" + keys[4], int(limits[0]), float(limits[1]), float(limits[2])
            )
        else:
            keys = line.split(",")
            h[hname].SetBinContent(int(keys[0]), float(keys[1]))
        n += 1
    return h
 
 
def manipulatePhysics(motherMode, mass, P8gen):
    # changes of the table, now it is deleted and we have each meson mother for each meson production
    # print motherMode
    if motherMode == "pi0" and pi0mass - mass >= 0.0000001:
        # use pi0 -> gamma A'
        selectedMum = 111
        P8gen.SetParameters("111:oneChannel = 1 1 0 22 9900015")
 
    elif motherMode == "eta" and etamass - mass >= 0.000001:
        # use eta -> gamma A'
        # print "eta"
        selectedMum = 221
        P8gen.SetParameters("221:oneChannel = 1 1 0 22 9900015")
 
    elif motherMode == "omega" and omegamass - mass >= 0.00001:
        # use omega -> pi0 A'
        # print "omega"
        selectedMum = 223
        P8gen.SetParameters("223:oneChannel = 1 1 0 111 9900015")
 
    elif motherMode == "eta1" and eta1mass - mass >= 0.00001:
        # use eta' -> gamma A'
        selectedMum = 331
        P8gen.SetParameters("331:oneChannel = 1 1 0 22 9900015")
        # should be considered also for mass < 0.188 GeV....
        # P8gen.SetParameters("331:oneChannel = 1 1 0 223 9900015")29%BR
        # P8gen.SetParameters("331:oneChannel = 1 1 0 113 9900015")2.75%BR
 
    elif motherMode == "eta11" and eta1mass - mass >= 0.00001:
        # use eta' -> gamma A'
        selectedMum = 331
        P8gen.SetParameters("331:oneChannel = 1 1 0 113 9900015")
        # should be considered also for mass < 0.188 GeV....
        # P8gen.SetParameters("331:oneChannel = 1 1 0 223 9900015")29%BR
        # P8gen.SetParameters("331:oneChannel = 1 1 0 113 9900015")2.75%BR
 
    else:
        # print "ERROR: please enter a nicer mass, for meson production it needs to be between %3.3f and %3.3f."%(pi0Start,eta1Stop)
        return -1
    return selectedMum
 
 
def configure(P8gen, mass, epsilon, inclusive, motherMode, deepCopy=False, debug=True):
    # configure pythia8 for Ship usage
    _exit_stack = contextlib.ExitStack()
    if debug:
        pythia_log = _exit_stack.enter_context(open("pythia8_conf.txt", "w"))  # noqa: SIM115
        P8gen = MethodLogger(P8gen, sink=pythia_log)
    P8gen.UseRandom3()  # TRandom1 or TRandom3 ?
    P8gen.SetMom(400)  # beam momentum in GeV
    if deepCopy:
        P8gen.UseDeepCopy()
    ROOT.TDatabasePDG.Instance()
    if inclusive == "meson":
        # let strange particle decay in Geant4
        p8 = P8gen.getPythiaInstance()
        n = 1
        while n != 0:
            n = p8.particleData.nextId(n)
            p = p8.particleData.particleDataEntryPtr(n)
            if p.tau0() > 1:
                command = str(n) + ":mayDecay = false"
                p8.readString(command)
                print("Pythia8 configuration: Made %s stable for Pythia, should decay in Geant4" % (p.name()))
 
        # Configuring production
        P8gen.SetParameters("SoftQCD:nonDiffractive = on")
 
    elif inclusive == "qcd":
        P8gen.SetDY()
        P8gen.SetMinDPMass(0.7)
 
        if mass < P8gen.MinDPMass():
            print("WARNING! Mass is too small, minimum is set to %3.3f GeV." % P8gen.MinDPMass())
            return 0
 
        # produce a Z' from hidden valleys model
        p8 = P8gen.getPythiaInstance()
        n = 1
        while n != 0:
            n = p8.particleData.nextId(n)
            p = p8.particleData.particleDataEntryPtr(n)
            if p.tau0() > 1:
                command = str(n) + ":mayDecay = false"
                p8.readString(command)
                print("Pythia8 configuration: Made %s stable for Pythia, should decay in Geant4" % (p.name()))
 
        # Configuring production
        P8gen.SetParameters("HiddenValley:ffbar2Zv = on")
        P8gen.SetParameters("HiddenValley:Ngauge = 1")
 
    elif inclusive == "pbrem":
        P8gen.SetParameters("ProcessLevel:all = off")
        # Also allow resonance decays, with showers in them
        # P8gen.SetParameters("Standalone:allowResDec = on")
 
        # Optionally switch off decays.
        # P8gen.SetParameters("HadronLevel:Decay = off")
 
        # Switch off automatic event listing in favour of manual.
        P8gen.SetParameters("Next:numberShowInfo = 0")
        P8gen.SetParameters("Next:numberShowProcess = 0")
        P8gen.SetParameters("Next:numberShowEvent = 0")
        proton_bremsstrahlung.protonEnergy = P8gen.GetMom()
        norm = proton_bremsstrahlung.prodRate(mass, epsilon)
        print("A' production rate per p.o.t: \t %.8g" % norm)
        P8gen.SetPbrem(proton_bremsstrahlung.hProdPDF(mass, epsilon, norm, 350, 1500))
 
    # Define dark photon
    DP_instance = darkphoton.DarkPhoton(mass, epsilon)
    ctau = DP_instance.cTau()
    print("ctau p8dpconf file =%3.6f cm" % ctau)
    print("Initial particle parameters for PDGID %d :" % P8gen.GetDPId())
    P8gen.List(P8gen.GetDPId())
    if inclusive == "qcd":
        dpid = P8gen.GetDPId()
        P8gen.SetParameters(f"{dpid}:m0 = {mass:.12}")
        # P8gen.SetParameters(str(P8gen.GetDPId())+":mWidth = "+str(u.mm/ctau))
        P8gen.SetParameters(f"{dpid}:mWidth = {u.hbarc / ctau:.12}")
        P8gen.SetParameters(f"{dpid}:mMin = 0.001")
        P8gen.SetParameters(f"{dpid}:tau0 = {ctau / u.mm:.12}")
        # P8gen.SetParameters("ParticleData:modeBreitWigner = 0")
        # P8gen.SetParameters(str(P8gen.GetDPId())+":isResonance = false")
        # P8gen.SetParameters(str(P8gen.GetDPId())+":all = A A 3 0 0 "+str(mass)+" 0.0 0.0 0.0 "+str(ctau/u.mm)+"  0   1   0   1   0")
        P8gen.SetParameters(f"{dpid}:onMode = off")
    else:
        P8gen.SetParameters(
            f"{P8gen.GetDPId()}:new = A A 3 0 0 {mass:.12} 0.0 0.0 0.0 {ctau / u.mm:.12}  0   1   0   1   0"
        )
        # if (inclusive=="pbrem"):
        
 
    P8gen.SetParameters("Next:numberCount    =  0")
 
    # Configuring decay modes...
    readDecayTable.addDarkPhotondecayChannels(
        P8gen,
        mass,
        DP_instance,
        conffile=os.path.expandvars("$FAIRSHIP/python/darkphotonDecaySelection.conf"),
        verbose=True,
    )
    # Finish DP setup...
    P8gen.SetParameters(f"{P8gen.GetDPId()}:mayDecay = on")
    # P8gen.SetDPId(P8gen.GetDPId())
    # also add to PDG
    gamma = u.hbarc / float(ctau)  # 197.3269631e-16 / float(ctau) # hbar*c = 197 MeV*fm = 197e-16 GeV*cm
    print("gamma=%e" % gamma)
    addDPtoROOT(pid=P8gen.GetDPId(), m=mass, g=gamma)
 
    if inclusive == "meson":
        # change meson decay to dark photon depending on mass
        selectedMum = manipulatePhysics(motherMode, mass, P8gen)
        print("selected mum is : %d" % selectedMum)
        if selectedMum == -1:
            return 0
 
    # P8gen.SetParameters("Check:particleData = on")
 
    _exit_stack.close()
 
    return 1