Functions
def	rotate (px, py, pz, theta, phi)

None	update_file (str filename, final_xsec)

None	inspect_file (str filename)

None	makeMuonDIS ()

def	getMasssq (pid)

def	rotate (ctheta, stheta, cphi, sphi, px, py, pz)

Variables
	level

r	PDG = r.TDatabasePDG.Instance()

argparse	parser = argparse.ArgumentParser(description=__doc__)

	help

	required

	dest

	default

	type

argparse	args = parser.parse_args()

argparse	n_events = args.n_events

argparse	first_mu_event = args.first_mu_event

list	headers

dict	Fixtarget = {1: "p+", 0: "n0"}

int	nJob = 2

int	nMult = 10

str	muonIn = "$SHIPSOFT/data/muConcrete.root"

int	nPerJob = 20000

ROOT	myPythia = ROOT.TPythia6()

ROOT	kc = myPythia.Pycomp(kf)

dict	masssq = {}

int	R = int(time.time() % 900000000)

dict	mutype = {-13: "gamma/mu+", 13: "gamma/mu-"}

ROOT	fout = ROOT.TFile("muonDis_" + str(nJob) + ".root", "recreate")

ROOT	dTree = ROOT.TTree("DIS", "muon DIS")

ROOT	iMuon = ROOT.TClonesArray("TVectorD")

ROOT	iMuonBranch = dTree.Branch("InMuon", iMuon, 32000, -1)

ROOT	dPart = ROOT.TClonesArray("TVectorD")

ROOT	dPartBranch = dTree.Branch("Particles", dPart, 32000, -1)

ROOT	fin = ROOT.TFile(muonIn)

ROOT	sTree = fin.muons

ROOT	nTOT = sTree.GetEntries()

int	nStart = nPerJob * nJob

min	nEnd = min(nTOT, nStart + nPerJob)

int	nMade = 0

ROOT	rc = sTree.GetEvent(k)

	px

	py

	pz

	x

	y

	z

	pid

	w

ROOT	p = ROOT.TMath.Sqrt(px * px + py * py + pz * pz)

ROOT	E = ROOT.TMath.Sqrt(getMasssq(pid) + p * p)

ROOT	theta = ROOT.TMath.ACos(pz / p)

ROOT	phi = ROOT.TMath.ATan2(py, px)

	ctheta

	stheta

	cphi

	sphi

array	mu = array("d", [pid, px, py, pz, E, x, y, z, w])

ROOT	tca_vec = iMuon.ConstructedAt(0)

ROOT	muPart = ROOT.TVectorD(9, mu)

ROOT	did = myPythia.GetK(itrk, 2)

	dpx

	dpy

	dpz

dpx	psq = dpx2 + dpy2 + dpz**2

array	m = array("d", [did, dpx, dpy, dpz, E])

ROOT	part = ROOT.TVectorD(5, m)

len	nPart = len(dPart)

Detailed Description

Script to generate DIS events for muons in Pythia6, and save them to a ROOT file (along with the original muon's soft interactions).

Function Documentation

◆ getMasssq()

def makeMuonDIS.getMasssq ( pid )

Definition at line 40 of file makeMuonDIS.py.

def getMasssq(pid):
    apid = abs(int(pid))
    if apid not in masssq:
        masssq[apid] = PDG.GetParticle(apid).Mass() ** 2
    return masssq[apid]
 
 

◆ inspect_file()

None makeMuonDIS.inspect_file ( str filename )

Inspect the contents of muonDis file.

Definition at line 106 of file makeMuonDIS.py.

def inspect_file(filename: str) -> None:
    """Inspect the contents of muonDis file."""
    file = r.TFile.Open(filename, "READ")
    tree = file.DIS
 
    table_rows = []
 
    for i, event in enumerate(tree):
        muon = event.InMuon[0]
        isProton = int(muon[9])
        fix_target = Fixtarget.get(isProton, "unknown")
        cross_sec = muon[10]
 
        nParticles = len(event.DISParticles)
        nSoftTracks = len(event.SoftParticles)
        nSBThits = len(event.muon_vetoPoints)
        nUBThits = len(event.muon_UpstreamTaggerPoints)
 
        table_rows.append([i, fix_target, nParticles, nSoftTracks, nSBThits, nUBThits, cross_sec])
 
    file.Close()
    logging.info("\n" + tabulate(table_rows, headers=headers, tablefmt="grid"))
 
 

◆ makeMuonDIS()

None makeMuonDIS.makeMuonDIS ( )

Generate DIS events.

Definition at line 130 of file makeMuonDIS.py.

def makeMuonDIS() -> None:
    """Generate DIS events."""
    final_xsec = {}
 
    logging.info(f"Opening input file: {args.inputFile}")
    muonFile = r.TFile.Open(args.inputFile, "read")
 
    try:
        muon_tree = muonFile.MuonAndSoftInteractions
    except Exception as e:
        logging.error(e)
        muonFile.Close()
        exit(1)
 
    logging.debug(f"Total entries in the tree: {muon_tree.GetEntries()}")
    last_mu_event = min(muon_tree.GetEntries(), first_mu_event + n_events)
 
    logging.info("Creating output file: muonDis.root")
 
    outputFile = r.TFile.Open("muonDis.root", "recreate")
    output_tree = r.TTree("DIS", "muon DIS")
 
    iMuon = r.TClonesArray("TVectorD")
    output_tree.Branch("InMuon", iMuon, 32000, -1)
 
    dPartDIS = r.TClonesArray("TVectorD")
    output_tree.Branch("DISParticles", dPartDIS, 32000, -1)
 
    dPartSoft = r.TClonesArray("TVectorD")
    output_tree.Branch("SoftParticles", dPartSoft, 32000, -1)
 
    muon_vetoPoints = r.TClonesArray("vetoPoint")
    output_tree.Branch("muon_vetoPoints", muon_vetoPoints, 32000, -1)
 
    muon_UpstreamTaggerPoints = r.TClonesArray("UpstreamTaggerPoint")
    output_tree.Branch("muon_UpstreamTaggerPoints", muon_UpstreamTaggerPoints, 32000, -1)
 
    myPythia = r.TPythia6()
    myPythia.SetMSEL(2)
    myPythia.SetPARP(2, 2)
    for kf in [211, 321, 130, 310, 3112, 3122, 3222, 3312, 3322, 3334]:
        kc = myPythia.Pycomp(kf)
        myPythia.SetMDCY(kc, 1, 0)
 
    seed = int(time.time() % 900000000)
    myPythia.SetMRPY(1, seed)
    mutype = {-13: "gamma/mu+", 13: "gamma/mu-"}
 
    myPythia.SetMSTU(11, 11)
    logging.info(f"Processing muon events from {first_mu_event} to {last_mu_event - 1}...")
 
    nMade = 0
 
    for k in range(first_mu_event, last_mu_event):
        DIS_table = []  # debug
        cross_sections = []
 
        muon_tree.GetEvent(k)
 
        imuondata = muon_tree.imuondata
 
        pid = imuondata[0]
        px = imuondata[1]
        py = imuondata[2]
        pz = imuondata[3]
        x = imuondata[4]
        y = imuondata[5]
        z = imuondata[6]
        w = imuondata[7]
        time_muon = imuondata[8]
        nmuons = imuondata[9]  # number of muons in the original MuBack event
 
        p = r.TMath.Sqrt(px**2 + py**2 + pz**2)
        mass = PDG.GetParticle(abs(int(pid))).Mass()
        E = r.TMath.Sqrt(mass**2 + p**2)
 
        theta = r.TMath.ACos(pz / p)
        phi = r.TMath.ATan2(py, px)
 
        isProton = 1
        xsec = 0
 
        mu = array(
            "d",
            [
                pid,
                px,
                py,
                pz,
                E,
                x,
                y,
                z,
                w,
                isProton,
                xsec,
                time_muon,
                args.nDIS,
                nmuons,
            ],
        )
        muPart = r.TVectorD(14, mu)
        myPythia.Initialize("FIXT", mutype[pid], "p+", p)  # target = "p+"
        myPythia.Pylist(1)
 
        for a in range(args.nDIS):
            if a == args.nDIS // 2:
                myPythia.Initialize("FIXT", mutype[pid], "n0", p)  # target = "n0"
                isProton = 0
                # logging.debug("Switching to neutron interaction")
 
            dPartDIS.Clear()
            iMuon.Clear()
            muPart[9] = isProton
            iMuon[0] = muPart
            myPythia.GenerateEvent()
            myPythia.Pyedit(1)
 
            for itrk in range(1, myPythia.GetN() + 1):
                xsec = myPythia.GetPARI(1)
 
                muPart[10] = xsec
                did = myPythia.GetK(itrk, 2)
                dpx, dpy, dpz = rotate(
                    myPythia.GetP(itrk, 1),
                    myPythia.GetP(itrk, 2),
                    myPythia.GetP(itrk, 3),
                    theta,
                    phi,
                )
                psq = dpx**2 + dpy**2 + dpz**2
                masssq = PDG.GetParticle(did).Mass() ** 2
                E = r.TMath.Sqrt(masssq + psq)
                m = array("d", [did, dpx, dpy, dpz, E])
                part = r.TVectorD(5, m)
                nPart = len(dPartDIS)
                if dPartDIS.GetSize() == nPart:
                    dPartDIS.Expand(nPart + 10)
                # dPartDIS.ConstructedAt(nPart).Use(part) #to be adapted later
                dPartDIS[nPart] = part
 
            cross_sections.append(xsec)
 
            dPartSoft.Clear()
 
            for softTrack in muon_tree.tracks:
                did = softTrack.GetPdgCode()
                dpx = softTrack.GetPx()
                dpy = softTrack.GetPy()
                dpz = softTrack.GetPz()
 
                psq = dpx**2 + dpy**2 + dpz**2
                masssq = PDG.GetParticle(did).Mass() ** 2
                E = r.TMath.Sqrt(masssq + psq)
 
                softx = softTrack.GetStartX()
                softy = softTrack.GetStartY()
                softz = softTrack.GetStartZ()
                time_ = softTrack.GetStartT()
 
                m = array("d", [did, dpx, dpy, dpz, E, softx, softy, softz, time_])
 
                part = r.TVectorD(9, m)
                nPart = len(dPartSoft)
                if dPartSoft.GetSize() == nPart:
                    dPartSoft.Expand(nPart + 10)
                # dPartSoft.ConstructedAt(nPart).Use(part) #to be adapted later
                dPartSoft[nPart] = part
 
            muon_vetoPoints.Clear()
 
            index = 0
            for hit in muon_tree.muon_vetoPoints:
                if muon_vetoPoints.GetSize() == index:
                    muon_vetoPoints.Expand(index + 1)
                hit.SetTrackID(0)  # Set TrackID to match for muon ID for new simulation
                muon_vetoPoints[index] = hit
                index += 1
 
            muon_UpstreamTaggerPoints.Clear()
 
            ubt_index = 0
            for hit in muon_tree.muon_UpstreamTaggerPoints:
                if muon_UpstreamTaggerPoints.GetSize() == ubt_index:
                    muon_UpstreamTaggerPoints.Expand(ubt_index + 1)
                hit.SetTrackID(0)  # Set TrackID to match for muon ID for new simulation
                muon_UpstreamTaggerPoints[ubt_index] = hit
                ubt_index += 1
 
            output_tree.Fill()
            DIS_table.append(
                [
                    a,
                    Fixtarget[isProton],
                    myPythia.GetN(),
                    len(dPartSoft),
                    len(muon_vetoPoints),
                    len(muon_UpstreamTaggerPoints),
                    xsec,
                ]
            )
 
        final_xsec[k] = xsec
 
        nMade += 1
        logging.debug(
            f"\nMuon index:{k} \n\tPID = {pid}, weight = {w}, converged_cross_section= {final_xsec[k]}\n\tpx = {px}, py = {py}, pz = {pz}, E = {E},\n\tx = {x}, y = {y}, z = {z}\n\n\tDIS Events Summary\n{tabulate(DIS_table, headers=headers, tablefmt='grid')}"
        )
        if nMade % 10 == 0:
            logging.info(f"Muons processed: {nMade}")
 
    outputFile.cd()
    output_tree.Write()
    myPythia.SetMSTU(11, 6)
    logging.info(
        f"DIS generated for muons (index {first_mu_event} - {last_mu_event - 1}) , output saved in muonDis.root, nDISPerMuon = {args.nDIS}"
    )
    outputFile.Close()
    muonFile.Close()
 
    update_file("muonDis.root", final_xsec)
 
 

◆ rotate() [1/2]

def makeMuonDIS.rotate	(	ctheta,
		stheta,
		cphi,
		sphi,
		px,
		py,
		pz
	)

Definition at line 66 of file makeMuonDIS.py.

def rotate(ctheta, stheta, cphi, sphi, px, py, pz):
    # rotate around y-axis
    px1 = ctheta * px + stheta * pz
    pzr = -stheta * px + ctheta * pz
    # rotate around z-axis
    pxr = cphi * px1 - sphi * py
    pyr = sphi * px1 + cphi * py
    return pxr, pyr, pzr
 
 

◆ rotate() [2/2]

def makeMuonDIS.rotate	(	px,
		py,
		pz,
		theta,
		phi
	)

Rotate the daughter particle momentum to align with respect to the muon's momentum.

Definition at line 55 of file makeMuonDIS.py.

def rotate(px, py, pz, theta, phi):
    """Rotate the daughter particle momentum to align with respect to the muon's momentum."""
    momentum = r.TVector3(px, py, pz)
 
    rotation = r.TRotation()
    rotation.RotateY(theta)  # Rotate around the Y-axis
    rotation.RotateZ(phi)  # Rotate around the Z-axis
 
    # Apply the rotation to the momentum vector
    rotated_momentum = rotation * momentum
 
    return rotated_momentum.X(), rotated_momentum.Y(), rotated_momentum.Z()
 
 

◆ update_file()

None makeMuonDIS.update_file	(	str	filename,
			final_xsec
	)

Update the DIS cross section of the muon to the converged value from Pythia.

Definition at line 69 of file makeMuonDIS.py.

def update_file(filename: str, final_xsec) -> None:
    """Update the DIS cross section of the muon to the converged value from Pythia."""
    file = r.TFile.Open(filename, "read")
 
    original_tree = file.DIS
 
    temp_filename = filename + ".tmp"
    temp_file = r.TFile.Open(temp_filename, "recreate")
 
    updated_tree = original_tree.CloneTree(0)
 
    for i, event in enumerate(original_tree):
        mu = event.InMuon[0]
        mu[10] = final_xsec[int(first_mu_event + i / args.nDIS)]
        updated_tree.Fill()
 
    updated_tree.Write("DIS", r.TObject.kOverwrite)
    temp_file.Close()
    file.Close()
 
    os.remove(filename)
    os.rename(temp_filename, filename)
    logging.info("Muon DIS events successfully updated with converged cross sections.")
 
 

Variable Documentation

◆ args

argparse makeMuonDIS.args = parser.parse_args()

Definition at line 50 of file makeMuonDIS.py.

◆ cphi

makeMuonDIS.cphi

Definition at line 100 of file makeMuonDIS.py.

◆ ctheta

makeMuonDIS.ctheta

Definition at line 99 of file makeMuonDIS.py.

◆ default

makeMuonDIS.default

Definition at line 29 of file makeMuonDIS.py.

◆ dest

makeMuonDIS.dest

Definition at line 26 of file makeMuonDIS.py.

◆ did

ROOT makeMuonDIS.did = myPythia.GetK(itrk, 2)

Definition at line 114 of file makeMuonDIS.py.

◆ dPart

ROOT makeMuonDIS.dPart = ROOT.TClonesArray("TVectorD")

Definition at line 58 of file makeMuonDIS.py.

◆ dPartBranch

ROOT makeMuonDIS.dPartBranch = dTree.Branch("Particles", dPart, 32000, -1)

Definition at line 59 of file makeMuonDIS.py.

◆ dpx

makeMuonDIS.dpx

Definition at line 115 of file makeMuonDIS.py.

◆ dpy

makeMuonDIS.dpy

Definition at line 115 of file makeMuonDIS.py.

◆ dpz

makeMuonDIS.dpz

Definition at line 115 of file makeMuonDIS.py.

◆ dTree

ROOT makeMuonDIS.dTree = ROOT.TTree("DIS", "muon DIS")

Definition at line 55 of file makeMuonDIS.py.

◆ E

ROOT makeMuonDIS.E = ROOT.TMath.Sqrt(getMasssq(pid) + p * p)

Definition at line 95 of file makeMuonDIS.py.

◆ fin

ROOT makeMuonDIS.fin = ROOT.TFile(muonIn)

Definition at line 62 of file makeMuonDIS.py.

◆ first_mu_event

argparse makeMuonDIS.first_mu_event = args.first_mu_event

Definition at line 52 of file makeMuonDIS.py.

◆ Fixtarget

dict makeMuonDIS.Fixtarget = {1: "p+", 0: "n0"}

Definition at line 103 of file makeMuonDIS.py.

◆ fout

ROOT makeMuonDIS.fout = ROOT.TFile("muonDis_" + str(nJob) + ".root", "recreate")

Definition at line 54 of file makeMuonDIS.py.

◆ headers

list makeMuonDIS.headers

Initial value:

=  [
    "DIS_index",
    "Fix_Target",
    "nParticles in event",
    "nSoftTracks_iMuon",
    "nSBThits_iMuon",
    "nUBThits_iMuon",
    "cross_sec",
]

Definition at line 94 of file makeMuonDIS.py.

◆ help

makeMuonDIS.help

Definition at line 22 of file makeMuonDIS.py.

◆ iMuon

ROOT makeMuonDIS.iMuon = ROOT.TClonesArray("TVectorD")

Definition at line 56 of file makeMuonDIS.py.

◆ iMuonBranch

ROOT makeMuonDIS.iMuonBranch = dTree.Branch("InMuon", iMuon, 32000, -1)

Definition at line 57 of file makeMuonDIS.py.

◆ kc

ROOT makeMuonDIS.kc = myPythia.Pycomp(kf)

Definition at line 34 of file makeMuonDIS.py.

◆ level

makeMuonDIS.level

Definition at line 16 of file makeMuonDIS.py.

◆ m

array makeMuonDIS.m = array("d", [did, dpx, dpy, dpz, E])

Definition at line 120 of file makeMuonDIS.py.

◆ masssq

dict makeMuonDIS.masssq = {}

Definition at line 37 of file makeMuonDIS.py.

◆ mu

array makeMuonDIS.mu = array("d", [pid, px, py, pz, E, x, y, z, w])

Definition at line 101 of file makeMuonDIS.py.

◆ muonIn

sys makeMuonDIS.muonIn = "$SHIPSOFT/data/muConcrete.root"

Definition at line 14 of file makeMuonDIS.py.

◆ muPart

ROOT makeMuonDIS.muPart = ROOT.TVectorD(9, mu)

Definition at line 107 of file makeMuonDIS.py.

◆ mutype

dict makeMuonDIS.mutype = {-13: "gamma/mu+", 13: "gamma/mu-"}

Definition at line 49 of file makeMuonDIS.py.

◆ myPythia

ROOT makeMuonDIS.myPythia = ROOT.TPythia6()

Definition at line 30 of file makeMuonDIS.py.

◆ n_events

argparse makeMuonDIS.n_events = args.n_events

Definition at line 51 of file makeMuonDIS.py.

◆ nEnd

ROOT makeMuonDIS.nEnd = min(nTOT, nStart + nPerJob)

Definition at line 79 of file makeMuonDIS.py.

◆ nJob

int makeMuonDIS.nJob = 2

Definition at line 12 of file makeMuonDIS.py.

◆ nMade

int makeMuonDIS.nMade = 0

Definition at line 87 of file makeMuonDIS.py.

◆ nMult

int makeMuonDIS.nMult = 10

Definition at line 13 of file makeMuonDIS.py.

◆ nPart

len makeMuonDIS.nPart = len(dPart)

Definition at line 123 of file makeMuonDIS.py.

◆ nPerJob

int makeMuonDIS.nPerJob = 20000

Definition at line 15 of file makeMuonDIS.py.

◆ nStart

int makeMuonDIS.nStart = nPerJob * nJob

Definition at line 78 of file makeMuonDIS.py.

◆ nTOT

ROOT makeMuonDIS.nTOT = sTree.GetEntries()

Definition at line 76 of file makeMuonDIS.py.

◆ p

ROOT makeMuonDIS.p = ROOT.TMath.Sqrt(px * px + py * py + pz * pz)

Definition at line 94 of file makeMuonDIS.py.

◆ parser

argparse makeMuonDIS.parser = argparse.ArgumentParser(description=__doc__)

Definition at line 21 of file makeMuonDIS.py.

◆ part

ROOT makeMuonDIS.part = ROOT.TVectorD(5, m)

Definition at line 121 of file makeMuonDIS.py.

◆ PDG

ROOT makeMuonDIS.PDG = r.TDatabasePDG.Instance()

Definition at line 17 of file makeMuonDIS.py.

◆ phi

ROOT makeMuonDIS.phi = ROOT.TMath.ATan2(py, px)

Definition at line 98 of file makeMuonDIS.py.

◆ pid

makeMuonDIS.pid

Definition at line 93 of file makeMuonDIS.py.

◆ psq

dpx makeMuonDIS.psq = dpx**2 + dpy**2 + dpz**2

Definition at line 118 of file makeMuonDIS.py.

◆ px

makeMuonDIS.px

Definition at line 91 of file makeMuonDIS.py.

◆ py

makeMuonDIS.py

Definition at line 91 of file makeMuonDIS.py.

◆ pz

makeMuonDIS.pz

Definition at line 91 of file makeMuonDIS.py.

◆ R

int makeMuonDIS.R = int(time.time() % 900000000)

Definition at line 47 of file makeMuonDIS.py.

◆ rc

ROOT makeMuonDIS.rc = sTree.GetEvent(k)

Definition at line 89 of file makeMuonDIS.py.

◆ required

makeMuonDIS.required

Definition at line 22 of file makeMuonDIS.py.

◆ sphi

makeMuonDIS.sphi

Definition at line 100 of file makeMuonDIS.py.

◆ stheta

makeMuonDIS.stheta

Definition at line 99 of file makeMuonDIS.py.

◆ sTree

ROOT makeMuonDIS.sTree = fin.muons

Definition at line 63 of file makeMuonDIS.py.

◆ tca_vec

ROOT makeMuonDIS.tca_vec = iMuon.ConstructedAt(0)

Definition at line 106 of file makeMuonDIS.py.

◆ theta

ROOT makeMuonDIS.theta = ROOT.TMath.ACos(pz / p)

Definition at line 97 of file makeMuonDIS.py.

◆ type

makeMuonDIS.type

Definition at line 30 of file makeMuonDIS.py.

◆ w

makeMuonDIS.w

Definition at line 93 of file makeMuonDIS.py.

◆ x

makeMuonDIS.x

Definition at line 92 of file makeMuonDIS.py.

◆ y

makeMuonDIS.y

Definition at line 92 of file makeMuonDIS.py.

◆ z

makeMuonDIS.z

Definition at line 92 of file makeMuonDIS.py.

Functions

Variables

Detailed Description

Function Documentation

◆ getMasssq()

◆ inspect_file()

◆ makeMuonDIS()

◆ rotate() [1/2]

◆ rotate() [2/2]

◆ update_file()

Variable Documentation

◆ args

◆ cphi

◆ ctheta

◆ default

◆ dest

◆ did

◆ dPart

◆ dPartBranch

◆ dpx

◆ dpy

◆ dpz

◆ dTree

◆ E

◆ fin

◆ first_mu_event

◆ Fixtarget

◆ fout

◆ headers

◆ help

◆ iMuon

◆ iMuonBranch

◆ kc

◆ level

◆ m

◆ masssq

◆ mu

◆ muonIn

◆ muPart

◆ mutype

◆ myPythia

◆ n_events

◆ nEnd

◆ nJob

◆ nMade

◆ nMult

◆ nPart

◆ nPerJob

◆ nStart

◆ nTOT

◆ p

◆ parser

◆ part

◆ PDG

◆ phi

◆ pid

◆ psq

◆ px

◆ py

◆ pz

◆ R

◆ rc

◆ required

◆ sphi

◆ stheta

◆ sTree

◆ tca_vec

◆ theta

◆ type

◆ w

◆ x

◆ y

◆ z