How to run the matrix element program COMPHEP at UC

How to install COMPHEP on one of the CDF machines

Getting the code

If you do not have a copy of the CompHEP package, download one from the official page. Currently a version of 4.1.10 is available locally at /cdf/data23a/COMPHEP/CompHEP_V_41.10.tar.gz. An unofficial but later version can be found at Version 4.2.

Unpacking and compiling the code

Once you have moved the tar file to a convenient location (for example, /cdf/data22d/foo/comphep), unpack and compile the code:

  mkdir /cdf/data22d/foo/comphep
  mv CompHEP_V_41.10.tar.gz /cdf/data22d/foo/comphep
  cd /cdf/data22d/foo/comphep 
  gzip -cd CompHEP_V_41.10.tar.gz | tar xvf -
  cd V_41.10
  ./mk_all

Setting the COMPHEP environment variable

Next, you should define the environment variable COMPHEP to point to your version of the code. In this example, we would use:

  setenv COMPHEP /cdf/data22d/foo/comphep/V_41.10

whereas you should substitute the appropriate location for your CompHEP directory. For added convenience, you should consider adding this command to your login file.

Creating a working area

Each time you wish to generate a new process, you will need to create a working directory and copy over any necessary files from the area pointed to by $COMPHEP. This last part is performed using the install command included with CompHEP. For example, we could create a directory called wgamma and perform an install there:

  mkdir /cdf/data22d/foo/comphep/wgamma
  cd /cdf/data22d/foo/comphep/wgamma
  $COMPHEP/install

You should see the following files and subdirectories:

  batch.pl*  comphep*  comphep.ini  models/  results/  tmp/  usr/  usrlib.a

If you do, you are now ready to start using CompHEP.

Running CompHEP

From the working directory that you created for your process, run the comphep executable:

  ./comphep &

This will launch a GUI interface from which you can perform all of the steps to calculate cross sections, generate events, and produce kinematic distributions. You may also choose to use the GUI for defining the process, cuts, and other properties but then use the batch capabilities for actually carrying out the calculations and event generation.

Selecting the process

The first step involves selecting a Model from menu. We have been using the _SM_ud model, which treats the two lightest generations of quarks as degenerate. After that, select Enter Process, and answer the prompts. For a process such as p pbar -> W + gamma, we would enter:

  _SM_ud
    Enter Process
      Enter process:                p,P -> e1,N1,A
      composite 'p' consists of:    u#,U#,d#,D#,G
      composite 'P' consists of:    u#,U#,d#,D#,G
      Enter CMS Energy in GeV:      1960
      Exclude diagrams with:

Creating the code for performing the calculations

The next several steps are designed to produce the code which will be used later for any numerical calculations. Briefly, you should do:

  View diagrams                   shows all the diagrams for our process

  Squaring technique
    View squared diagrams         
    Symbolic calculations         
    Write results                 
      C code (for num.calc)       outputs the necessary code
    C-compiler                    compiles the code

At this point, a terminal window should appear which shows the compilation of the C code. If compilation is successful, this window will close and a second CompHEP window will pop up.

PDF's, cuts, etc.

In this second window, we set the PDFs for the incoming particles, the event cuts to be applied, and any regularization (to avoid blowing up at the poles in the cross section calculation). For these menus, use the F1 key for help and ESC to return to the previous menu.

For the PDF's, we use CTEQ 5L:

  IN state
    S.F.1:
      PDF
        PDF:CTEQ5L(proton)
    S.F.2:
      PDF
        PDF:CTEQ5L(anti-proton)

The possible kinematic cuts are detailed in the Help menu (F1). These include Tx (transverse momentum, P_t, for particle x), Mxy (invariant mass of x & y), Jxy (jet cone angle between x & y), Nx (pseudorapidity for x). Note that the particle indices x and y correspond to the ordering of the particles in the subprocess, not the process.

For a process such as W gamma (where the subprocess looks like u# D# -> e1, N1, A), we might use cuts that look like:

  Cuts
    Parameter  |> Min bound <|> Max bound <|
    T3         |25.0         |             |
    T4         |25.0         |             |
    T5         |25.0         |             |
    N3         |-1.0         |1.0          |
    N4         |-1.0         |1.0          |
    N5         |-6.0         |6.0          |
    J34        |0.4          |             |

The regularization is a slightly more tricky issue. For situations where the cuts are extremely loose or completely undefined, the regularization appears to be necessary to avoid peaks in the cross section where the calculation might blow up. This is true for W's and Z's, but also for photon colinearity.

Entries in this table may look like:

  Regularization
    Momentum    |> Mass  <|> Width <| Power|
    45          |MW       |wW       |2     |
    13          |0        |0        |1     |
    23          |0        |0        |1     |
    34          |0        |0        |1     |
    345         |MW       |wW       |2     |

However, we have found that if a more restrictive set of cuts is used, entries in the regularization table may cause the calculation to crash, with the message ERROR IN REGULARIZATION. If such a message should appear, the job may appear to continue unabated, but surrender all hope of any results being produced. I have found that, in general, it is best to work without anything defined here if any reasonable set of cuts is used.

To batch or not to batch

Once you reach this stage, you may continue with the interactive GUI or exit to use the batch mode. If you continue with the GUI, then do:
(more explanation is forthcoming; ask for details)

  Vegas
    Start integration
    Clear statistics
    nCall = 50000
    Set Distributions
    Start integration
    Generate events
      Start search of maxima
      Number of events=(set value)
      Launch generator

Or, if you decide to go with the batch capability, exit the GUI completely using F9 in both windows. Then, from the working directory for your process (not the results directory):

  ./batch.pl
  ./batch.pl -run vegas     (starts the iterations to calc. crossection)
                            (keep doing the above until you get to 1% or so)
  ./batch.pl -run           (performs  "-run vegas,max,evnt"  all at once)

How to run COMPHEP

To select, use arrows to highlight selection; ENTER to chose. ESC takes one back one menu. If the right choice is already highlighted, just hit ENTER. At the bottom of the main window (first to pop up) are functions you can click on- they are very useful.

cd ~/myareawhateveritis/COMPHEP/V_41.10/USER/
comphep (should pop up a menu)
Now pick a model: either SM (Feynman gauge), or, (and this is the best choice for many things), _SM_ud, e.g. (this is best for Fermilab- uses two generations of light quarks (ud, cs), but doesn't distinguish them so that the number of diagrams is much smaller).
Enter Process
Type in particles for process: e.g., p,P -> e1,E1 for DY production.
In answer to"composit p consists of": u#,U#,d#,D#,G (the # means 2 gens)
In answer to"composit P consists of": u#,U#,d#,D#,G
Root s is 1960 (better than 2 TeV, per David Goldstein)
In answer to "Exclude diagrams with" one could type, in our DY example from above, A to exclude the photon in the s-channel, or Z, similarly.
View diagrams. Look at them, and also make sure the the right hand column has the right number of subprocesses (e.g., in our DY example, should be 8 if have both A and Z, e.g.; only 4 if only A or Z; never should be zero). N.B.- the numbers next to the particles on the diagrams are wrong- ignore them. The numbering is determined by counting in the `subprocess' line (not the `process' line).
Use escape to get back to menu with "squaring technique"
squaring technique (you can look at squared diagrams if you wish, if the process isn't too complicated. Up to 2->4 is ok; but, for example, enu+3 jets is too complicated).
symbolic calculations
write results (don't forget this step!!!)
C-code
C-compiler
IN state
SF 1
PDF then pick your PDF (e.g. CTEQ5L for the proton)
SF2
PDF then pick your PDF (e.g. CTEQ5L for the anti-proton)
Cuts (set cuts so it doesnt blow up: e.g. in our DY example, M34 20 500, or, for a process in which the 5th particle is a radiated gluon, T5 10. One may not need any cuts. See note below on cuts.)
Regularization (see F1- put in poles to be careful- e.g 34 MZ wZ 2 ; If there is a gauge cancellation of massless singularities, the power is 1; otherwise 2. Ask your local theorist. (My simple-minded (and wrong!) formulation is: For out-going photons and gluons power is 1; for other cases (virtual exchanges, e.g., or out-going fermions, unless they're connected by a gauge boson to an incoming line in a way that has a gauge cancellation, power is 2)
Vegas
Start integration (this makes the grid- do a couple of runs and see that convergence looks good)
Clear statistics
ncall = 50000 (e.g.)
set distributions (e.g. M34 20 500)
start integration
could do generate events now (unweighted events)
If you do generate events, check that there is an event_?.txt file with the right number of events before bailing out of this session.
To step through multiple subprocesses for a signature (e.g. pP->eE has 4 subprocesses, uU,Uu,dD,Dd), stay in the 2nd comphep window, and hit escape all the way back up to the top (if it asks you if you want to end the session, you've gone too far- say `no' (!), and you're now in the right place.) Select `subprocess'- move down one to the next subprocesses, and procede as above. The session number will be incremented; the plots, however, just continue with their numbering being bumped (i.e. if you are saving 8 plots per subprocess, the first subprocess will have plots 1-8, 2nd 9-16, etc.).
The log file for each session (e.g. session n) is named prt_n. This records what values you've typed in (VERY good to check).

How to Run CompHEP Numerical Calculations in Batch Mode

If you installed the version of CompHEP with batch mode, the file batch.pl should be in your working directory (e.g. comp_user). You may need to correct the path to perl on your computer in the first line of the script, and/or change the permission of the batch.pl file (chmod +x batch.pl).

The perl script ./batch.pl runs the numerical part of the calculations without the GUI. You can start this script after doing all the symbolic calculations. This script runs the n_comphep file in the subdirectory results. For help type
./batch.pl --help.

First setup the parameters of the calculations; The easiest way is with the GUI:
cd results
./n_comphep

Set the parameters you want (PDF, Cuts, ...) for the first subprocess. Exit the GUI. Type:

cd ../
./batch.pl

The file results/batch.dat will be created automatically with the same parameters for the all of subprocesses.

If you want to change parameters for some specific subprocesses, run the GUI again, choose the subprocess, set the new parameters, and exit to save the session.dat file. Alternatively, another way of saving the session.dat file is by stepping to the last step in the menu, Vegas, and then stepping back. You can add the parameters from the session.dat file for this particular subprocess by:
./batch.pl -add_ses2bat

The default number of generated events is 10000; to change this, edit the line in the results/batch.dat file:
#Events 100 1 0.200000 2.000000 10000
This needs to be changed for each of the subprocesses as appropriate; you can look at the cross-sections per subprocess (see below) to tailor the number of events per subprocess.

After setting the parameters for all of the subprocesses you are ready to run the calculations.
./batch.pl -run vegas
This command (see below) makes the grid and gives you an estimate of the precision- keep doing this step until you get to 1\% or so. Then:
./batch.pl -run
This command starts the steps of the calculations for all subprocesses. If you want to calculate only the cross section, run
./batch.pl -run vegas
which calculates the crossections and the Vegas grid, and stores the crossections in the protocol file results/prt_n and the grid in results/batch.dat. You can run selected subprocesses:
./batch.pl -run vegas -proc 3,1,4-6,5
(the numbers can be in any order)
This calculates the cross section for the subprocesses 1,3-6

If you want to generate events after you get the cross section with good precision, you can run:
./batch.pl -run max,evnt

If you want to know more about the options, run:
./batch.pl -h (or --help for the long help)

Notes and tips (the order is somewhat random- sorry!):

Hitting the key F1 will bring up help on the step you're on.
Hitting the key F2 will bring up general help (the manual).Hitting ESC until the top, and then F9 will get you out cleanly.
Control D deletes a line in one of the parameter boxes
typing n_comphep from one of the renamed results directories restarts the numerical part of comphep in that process.
you can copy the results dir to another name by hand, but you have to make a new results dir afterward. Or, clicking on `F6' in the main menu box will allow one to rename it, and then it automatically makes a new one for you (this is the preferred way).
There are two parts to COMPHEP, each with its own window. The one that comes up first allows you to see the diagrams and squared diagrams; the second, which appears when you `do integration', is the numerical part. Once you have proceded to the second, you can't look at the diagrams any more.
When you step through the diagrams, use the right arrow to go to the next. `Page down' or up will move you also. `O' toggles the DEL or not delete- you may (will) want to delete diagrams which are negligble, such as 2nd order weak diagrams (2 W's) in W+3 jet production.
`D' Deletes all diagrams when you're viewing diagrams. Useful if you want to delete all diagrams except those in one subprocess.
To duplicate a line in `CUTS', CNTRL-D, and then ENTER. TAB moves cursor.
If the subprocess window suddenly vanishes, look at the Xwindow you started from for an explanation. The most common reason I've found is that some divergence has bit you. This can either be from something in the regularization or in the cuts.
An example on the cuts: making W+gamma (really enu+gamma). One subprocess is U# d# -> A, e1,N1. Cuts that seem ok are: M45>20, T3>1, S34>1 (this keeps the photon from being colinear and infrared).
An example on regularization: making W+gamma (really enu+gamma)

 Regularization 
 Momentum    |> Mass  <|> Width <| Power|
45           |MW       |wW       |2      
13           |0        |0        |1      
23           |0        |0        |1      
34           |0        |0        |1      
345          |MW       |wW       |2      
=========================================

Order is important: for example, in the W+gamma example, the subprocess U#d#->A,e1,N1 has a cross-section of 5.2 pb, while d#U#->A,e1,N1 has a cross-section of 38.4 pb (the latter has the ubar in the pbar). Note also the possibility of losing the factor of 2 for W+ and W- (as well as getting the effect of the any EW detector asymmetry such as the chimney wrong.)
You must have a `results' subdirectory in USER. If you don't, when you click `write results' the window vanishes. Just mkdir results.
If after entering one subprocess with the first GUI you copy the session.dat file to a session_save.dat file, you can start the numerical integration over again from it without having to go back to the first gui again by deleting everything in the results file except the very initial files, and renaming session_save.dat, and then rerunning batch.pl. Similarly, when restarting, one can edit the session.dat file. It seems, however, that once you've run the numerical part you have to start back at this juncture (which isn't hard). HF

To start a new session, or resume an old one

To start a new session, in /cdf/data23a/COMPHEP/V_41.10/USR/ (i.e., where COMPHEP was installed) type `comphep'
To resume or change an existing process, in the subdirectory (for example, /cdf/data23a/COMPHEP/V_41.10/USR/W+1j/subprocess_1/) type `../n_comphep'. This will pick up the file `session.dat', and start from there. Remember, once you are in the numerical part of comphep you can't go back and look at diagrams- you need to start at the top.
The file session.dat is an ascii file that tells you what you did.
The recommended way to step through the subprocesses of a signature (e.g., in pP->tT, one has the subprocesses u#U#->tT, d#D#->tT, GG->tT, is to create a directory with the process name (e.g. tT), and then subdirectories of that for each of the subprocesses (e.g., in this case 3 subdirectories uU, dD, and GG). One can run in the process directory, and then use F6 in the first (main) window to rename each subprocess.

To make an input file to Pythia, and then have Pythia write it out in STDHEP format

setenv CPYTH62 /cdf/data23a/COMPHEP/cpyth62/
You will need to copy the pythia library or to compile pythia. In cpythia/pythia62/ do `f77 -c pythia6203.f' Will make pythia6203.o. This version is the first to implement the (famous) Les Houches accord.
The first step will be to mix the subprocess events with the right cross sections. This is done by the program mixPEV: outputs a .PEV file with the right mix. It mixes until one of the input files is exhausted, and then stops. It tells you how many events it has outputted.
In the process directory, make a subdirectory: mkdir EVENTS
From each of the subprocess subdirectories, move (mv) the event files into EVENTS. E.g., In the process directory /cdf/data23a/COMPHEP/V_41.10/USR/tT-eNEnbB_narrowMOdel, there are the subdirectories EVENTS/, GG/, dD/, and uU/. In GG, there is the file events_6.txt, where the 6 refers to the session number. Move this file into EVENTS/events_GG.txt. Similarly, from dD move the events file (in this case already renamed as events_dD.txt) into EVENTS/events_dD.txt; do similarly for uU.
cd EVENTS
run mixevents: type '$CPYTH62/mixPEV events_dD.txt events_dU.txt events_GG.txt'
Makes a new file, Mixed.PEV. For tT dileptons, 118K events gave 44 MB. It also (kindly!) makes a log file in Prt.PEV.
mv Mixed.PEV tT_dilepton.PEV
cp tT_dilepton.PEV $CPYTH62/interf62/tT_dilepton.PEV
cd $CPYTH62/interf62/
edit INPARM.DAT- rename file, change the process string (just text), and the number of events. Warning- the number of events MUST be less than the number generated-- if not, the program will just wrap around without complaining (so you get duplicate events!).
$CPYTH62/makePGEN
Check that there are no files names fort.66, fort.67, etc., leftover from the last time. Pythia won't be able to write if there are. Rename or move them if so.
PGEN.exe (i.e., run it)
When PGEN finishes, it will make files with names like fort.66, fort.67, ... Each of these will have 50K events in STDHEP format, written by Steve's routine (below) in binary (we started in ASCII, but the files were larger).
Move the fort.66 etc. files out of $PYTH62 to an appropriate directory.

Steve Mrenna's Pythia STDHEP Format

This is written by the subroutine pywrite.f, which is called from Pythia's main.f. The subroutine can be found in /cdf/data23a/COMPHEP/cpyth62/interf62. Here it is:

    SUBROUTINE PYWRITE 
      IMPLICIT NONE
      INTEGER I,J,KNT
      INTEGER LUNIT  
      INTEGER NMXHEP 
      INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP 
C...HEPEVT commonblock. 
      PARAMETER (NMXHEP=4000) 
      COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), 
     &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) 
      DOUBLE PRECISION PHEP,VHEP 
      DATA KNT/0/ 
      SAVE KNT 
      SAVE /HEPEVT/ 

      KNT=KNT+1 
      LUNIT=66+INT(KNT/50000) 
       WRITE(LUNIT) NEVHEP,NHEP 
      DO I=1,NHEP 
       WRITE(LUNIT) ISTHEP(I),IDHEP(I),(JMOHEP(J,I),J=1,2), 
     $ (JDAHEP(J,I),J=1,2),(PHEP(J,I),J=1,5),(VHEP(J,I),J=1,4) 
      ENDDO 
      RETURN 
      END

Reading STDHEP Format with AC++, and running CDFSIM

To read STDHEP format into a Hepg bank do the Following, (Stan Thompson mostly responsible for this):

              Easiaest thing is to copy the executable from cdf23 computer;
                  /home/cdf/carron/4.2.0/bin/Linux2-KCC_4_0/cdfGen
              and the .tcl file to run it:
                  /home/cdf/carron/4.2.0/bin/Linux2-KCC_4_0/stdhep.tcl
              In this tcl file you have to change the input file and output file to your
              preference(it is self evident). You also have to adjust the number of events;
              always put one less than the number of events in the binary file.

This will create a .root file that is hepg, this file can
be inputed in ntuples or cdfSim in the usual way.

              To actually build the executable yourself. for this do the following:
                 newrel -t 4.2.0 4.2.0
                 cd 4.2.0
                 addpkg generatorMods
                 copy the following files (these files have not been updated
                 in generatorMods yet) in the corresponding directories:
                 /home/cdf/carron/4.2.0/generatorMods/generatorMods/StdhepInputModule.hh
                 and
                 /home/cdf/carron/4.2.0/generatorMods/src/StdhepInputModule.cc
                                                                                        HadScatGenSequence.cc
                                                                                        hepfil.F
                                                                                        hepred.F
                                                                                        hepcls.F
                Then from base directory of the release do: gmake, this will build cdfGen,
                 Then run cdfGen with the .tcl file mentioned before.

Inputing Hepg bank into the Duke Ntuple.

                First you need to build the Duke Ntuple, for that, from base directory of your relese do:
                 addpkg -h dukehkg
                 gmake clear
                 gmake dukehkg.all
                In the bin directory now you should find an executable called myeleNew. This is theone. To
                run it copy the .tcl filr from cdf23 computer: /home/cdf/carron/4.2.0/TalktoSin.tcl to the
                bin directory and do: myeleNew TalktoSim.tcl >&logfile. You must remember to edit the
                .tcl file to have as input your Hepg bank file, and edit the output file location and name to your
                preference.
                 After this you will end up with a root file containing the filled Duke ntuple, the parton information
                 is in the branch HEPG.

Dumping the Hepg events from ntuple and generating plots with root macro.

              If you have processed the events with a different ntuple from the Duke Ntuple, you will need to
              generate the root skeleton and then copy the relevant parts of the code below. At the en of this
              section there are the instructions.

              Copy the files from cdf23 computer: /cdf/data23a/MC_data/CPythia_HEPG/ttbar/ttbar_hepg.C
              and /cdf/data23a/MC_data/CPythia_HEPG/ttbar/ttbar_hepg.h (This is a sample analysis for the ttbar)
              You MUST edit the .h file above to have the correct input ntuple file to the macro.
              Also edit the .C file at the end to have the name of your output file for your plots (don't overwrite my files!!!)
              open a root session and in the root prompt do:
                Root > gSystem->Load("$ROOTSYS/lib/libPhysics.so");
                Root > .L ttbar_hepg.C
                Root > ttbar_hepg t
                Root > ofstream fout("nameofyourfile.txt", ios_base::app);
                t.GetEntry(12); (Here you select the event you want to Dump)
                Root > t.dump_hepg(fout);
              When you exit root you can see that you have created a text file with the details of the Hepg bank for
              the event selected.

               Now to create the Plots (I am assuming you have edited already the name of the input file and output
               histogram file) open a root session and in root do:
                Root> gSystem->Load("$ROOTSYS/lib/libPhysics.so");
                Root > .L ttbar_hepg.C
                Root > ttbar_hepg t;
                Root > t.Loop();
               The plots will be created (You can alted the .C file to have the plots the way you want them, or select
               the particles you want, this is just an example for the ttbar case).
               To view the plots do: root filename.root
                in root: Root > TBrowser browser
               and clic your way to the plots.

If Ntuple different from Duke: (Comming soon ...)

Inputing the Hepg bank into the Sandard Ntuple (Comming soon...)
Inputing the Hepg bank into cdfSim (Comming soon...)

Links to COMPHEP home page, documentation:

COMPHEP home page
UC CDF Group home page
Pythia home page

Where Monte Carlo Files are Kept

For now, keep 3 kinds of files- the raw output (text) from COMPHEP, the session.dat file from COMPHEP, and the output (STDHEP format a la Mrenna, binary) from Pythia.
These are in /cdf/data23a/ in respective directories
We will want to make a systematic list of which we've done and who and who and when etc.
Intent is to make all this available to anybody who wants it (of course!).
To do this, we need disk space at Fermilab (write your congressman!). Have started with some files on fcdfsg2 /cdf/data04/s5/top/COMPHEP. Have also copied the code there (but not tried it).

Wish List

This is a wish list for features that would be handy. It's made out of ignorance-- some of these features may already exist; others perhaps should never exist. Apologies for the forwardness in making it -- feel free to edit it, adding explanations on how to do these things, or on why one shouldn't want to do them. HJF

Big Wishes

To be able to drive COMPHEP from a script, rather than the GUI. For example, once one has one subprocess done, one could make a script to do all the others. This is important for a hadron collider- it will substantially cut down the chance of mistakes, and speed up generating processes with many subprocesses to the point where it's feasible. It will also allow running COMPHEP on farms where interactive control isn't possible.
The ability to write out ALL the plots as postscript files and/or tables with a single command.

Little Wishes

A little program to read all the tables in a directory (or a list) and make postscript files of them.
Having the program know what p,P, and j are made of (u#,U#,d#,D#,G) when you're in the _SM_ud model.
Not losing the list of plots requested when it blows up during integration.
Deleting or fixing the numbers labelling the lines on the Feynman diagrams.
Once a plot has been declared log, it would be helpful if it stayed that way through all subprocesses (and the minimum would stay what it was previously).

Created by:

Lev Dudko,
Henry J. Frisch (frisch @ hep . uchicago . edu),
Sebastian Carron (carron @ phy . duke . edu), and
Chadd Smith (chadd @ cdf . uchicago . edu)

Last updated May 20, 2003 by Henry Frisch