Branching Ratio

Asked by NADY

Hi Johan,
I implemented new model with FeynRules1,6
and used it in MG5
 I generated this process p p >Zprime .
How i can get all branching ratio of Zprime from this process where i want to export the matrix elements source code to Pythia8150
Cheers
Nady

Question information

Language:
English Edit question
Status:
Answered
For:
MadGraph5_aMC@NLO Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Arian Abrahantes (arian-abrahantes) said :
#1

Dear Nady:

I am not an expert, I am not sure if you can do both things branching and cross section in a single run but if you split these calculations I am pretty sure it is possible (actually I did it myself). It only needs some extra work frm your side nt on MG team

regarding the partial decay width calculation (branching ratio), you can do the calculation directly in MG just make a process card where you put all possible decay process you have coded like

process zprime > x y @1
add process zprime > a b @2
....

as many as needed to compute the full decay width, then calculate branchings, I am pretty sure pythia has some command to include the new particles and branchings

hope it helps,

Revision history for this message
Johan Alwall (johan-alwall) said :
#2

Hello Nady,

Arian is exactly right, you simply generate all decay processes as he describes. You can then use the attached scripts calculate_decay_widths and collect_decay_widths.py. Simply put the scripts in the bin/ directory of the resulting process directory, make sure that you can run them (chmod +x calculate_decay_widths;chmod +x collect_decay_widths.py), and use generate_decay_widths instead of generate_events. The result is a file with all decay widths and branching ratios, which you can then copy and paste into your param_card. This can then be used directly in Pythia 8.

Cheers,
Johan

Revision history for this message
NADY (nady-bakhet) said :
#3

Hi Johan
Thanks Arian and Johan
I didn't find the scripts calculate_decay_widths , collect_decay_widths.py
and generate_decay_widths in bin directory of the resulting process directory
Where i can find it
Cheers,
Nady

Revision history for this message
Johan Alwall (johan-alwall) said :
#4

calculate_decay_widths:

#!/bin/bash
#
# This runs survey,refine,unweight_events, one after the other
#
# First we need to get into the main directory
#
#
# Usage: generate_events compression events parallel [name]
#
# parallel is 0 for serial, 1 for PBS, 2 for multicore
# name is queue name for PBS or number of cores for multicore
#

function error_exit {
    if [[ -e error ]]; then
 date
 cp -f error status
 rm survey refine refine2 >& /dev/null
 rm RunWeb
 $dirbin/gen_crossxhtml-pl $t
 $dirbin/gen_cardhtml-pl
 exit
    fi
}

if [[ ! -d ./bin ]]; then
    cd ../
    if [[ ! -d ./bin ]]; then
        echo "Error: it must be executed from the main, or bin directory"
        exit
    fi
fi
#
# Now let shell know where to find important executables
#
main=`pwd`
# Strip off last directory
mainpar=${main%/*}
dirbin=$main/bin
pydir=$mainpar/pythia-pgs/src
pgsdir=$pydir
delphesdir=$mainpar/Delphes
ERAdir=$mainpar/ExRootAnalysis
MAdir=$mainpar/MadAnalysis
webbin=$dirbin
td=$mainpar/td
web=0

echo $$ >> myprocid

# For Linux
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$main/lib
# For Mac OS X
export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:$main/lib

if [[ "$1" == "" ]]; then
    echo 'Enter 2 for multi-core, 1 for parallel, 0 for serial run'
     read mode
else
     mode=$1
fi
n=MadEvent
if [[ $mode -gt 0 ]]; then
   if [[ "$2" == "" ]]; then
     if [[ $mode -eq 1 ]]; then
       echo 'Enter name for jobs on pbs queue'
     elif [[ $mode -eq 2 ]]; then
       echo 'Enter number of cores'
     fi
       read n
   else
       n=$2
   fi
   if [[ "$3" == "" ]]; then
       echo 'Enter run name'
       read t
   else
       t=$3
   fi
else
   if [[ "$2" == "" ]]; then
      echo 'Enter run name'
      read t
   else
      t=$2
   fi
fi
#set t = TeV2
if [[ "$4" -ne "" ]]; then
    web=1
    webbin="$MADGRAPH_BASE/MG_ME/WebBin"
    pydir="$webbin/pythia-pgs"
    pgsdir=$pydir
    delphesdir="$webbin/Delphes"
    ERAdir="$MADGRAPH_BASE/MG_ME/ExRootAnalysis"
    MAdir="$MADGRAPH_BASE/MG_ME/MadAnalysis"
    td="$MADGRAPH_BASE/MG_ME/td"
    touch Online
fi
date
#
# Check if run already exists. If so, store run w/ new name
# and remove old run before starting.
#

if [[ -e status ]]; then
  rm status
fi
if [[ -e error ]]; then
  rm error
fi
touch RunWeb
echo "Cleaning directories" > status
cat status
$dirbin/gen_crossxhtml-pl $t
$dirbin/clean

# Set gridpack in run_card.dat
sed -i.bak "s/^[^#].*=.*gridpack / .true. = gridpack /g" Cards/run_card.dat

# Compile everything in Source
echo "Compiling Source" > status
cat status
$dirbin/gen_crossxhtml-pl $t
$dirbin/compile_Source
# If any compilations failed
error_exit

touch survey
echo "Starting jobs" > status
cat status
$dirbin/gen_crossxhtml-pl $t
$dirbin/survey $mode $n $t
error_exit
#
# Now collect the events
#
echo "Combining Events" >& status
cat status
$dirbin/gen_crossxhtml-pl $t
pushd SubProcesses > /dev/null
$dirbin/run_combine $mode
mv events.lhe ../Events/
mv unweighted_events.lhe ../Events/
popd > /dev/null
#
# do the banner
#
cd ./Events
echo "putting the banner"
$dirbin/put_banner events.lhe
$dirbin/put_banner unweighted_events.lhe

if [[ -e unweighted_events.lhe ]]; then
    $dirbin/extract_banner-pl unweighted_events.lhe banner.txt
fi

cd ..

$dirbin/store $t
$dirbin/gen_crossxhtml-pl $t

python $dirbin/collect_decay_widths.py $t

rm -f RunWeb
echo " " >& status
$dirbin/gen_crossxhtml-pl $t
$dirbin/gen_cardhtml-pl
date

Revision history for this message
Johan Alwall (johan-alwall) said :
#5

collect_decay_widths.py:

#!/usr/bin/python
#
# Collect the decay widths and calculate BRs for all particles, and put in param_card form

import glob
import optparse
import os
import re

root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]

# Main routine
if __name__ == "__main__":

    # Find available runs
    run_names = glob.glob(os.path.join(root_path, "Events", "*_banner.txt"))
    run_names = [os.path.basename(name) for name in run_names]
    run_names = [re.match("(.*)_banner.txt", name).group(1) for name in run_names]
    usage = "usage: %prog run_name"
    usage +="\n where run_name in (%s)" % ",".join(run_names)
    parser = optparse.OptionParser(usage=usage)
    (options, args) = parser.parse_args()
    if len(args) != 1 or args[0] not in run_names:
        if len(run_names) == 1:
            args = run_names
        else:
            parser.error("Please give a run name")
            exit
    run_name = args[0].strip()

    print "Collecting results for run ", run_name

    # Start collecting results
    subproc_dirs = \
           open(os.path.join(root_path, "SubProcesses", "subproc.mg")).read().split('\n')

    particle_dict = {}
    for subdir in subproc_dirs[:-1]:
        subdir = os.path.join(root_path, "SubProcesses", subdir)
        leshouche = open(os.path.join(subdir, 'leshouche.inc')).read().split('\n')[0]
        particles = re.search("/([\d,-]+)/", leshouche)
        if not particles:
            continue
        particles = [int(p) for p in particles.group(1).split(',')]
        results = \
             open(os.path.join(subdir, run_name + '_results.dat')).read().split('\n')[0]
        result = float(results.strip().split(' ')[0])
        try:
            particle_dict[particles[0]].append([particles[1:], result])
        except KeyError:
            particle_dict[particles[0]] = [[particles[1:], result]]

    lines = []
    for key in sorted(particle_dict.keys()):
        width = sum([r for p,r in particle_dict[key]])
        lines.append("#\n# PDG Width")
        lines.append("DECAY %i %e" % (key, width))
        if not width:
            continue
        lines.append("# BR NDA ID1 ID2 ...")
        brs = [[val[1]/width, val[0]] for val in particle_dict[key] if val[1]]
        for val in sorted(brs, reverse=True):
            lines.append(" %e %i %s" % (val[0], len(val[1]),
                                       " ".join([str(v) for v in val[1]])))
    output_name = os.path.join(root_path, run_name + "_decays.dat")
    decay_table = open(output_name, 'w')
    decay_table.write("\n".join(lines))
    print "Results written to ", output_name

Revision history for this message
NADY (nady-bakhet) said :
#6

Thanks Johan

I found " generate_events." only in bin directory of process
and didn't find generate_decay_widths file

Cheers,
Nady

Revision history for this message
Johan Alwall (johan-alwall) said :
#7

I just posted the two files you need here.

Cheers,
Johan

Can you help with this problem?

Provide an answer of your own, or ask NADY for more information if necessary.

To post a message you must log in.