2-2, decay 1-2, output to pythia8 program, the program is wrong

[tang yilei]

Last time, I reported a bug in the same headline. I just downloaded madgraph v5.1.8, it seems that you didn't fix this problem radically. Now, the function turns to such:

void Sigma_vectorlike_gcx_d5barelm_elm_sux::setIdColAcol()
{
  if(id1 == -4 && id2 == 21)
  {
    // Pick one of the flavor combinations
    int flavors[1][3] = {{-6001, 1000011}};
    vector<double> probs;
    double sum = matrix_element[1];
    probs.push_back(matrix_element[1]/sum);
    int choice = rndmPtr->pick(probs);
    id3 = flavors[choice][0];
    id4 = flavors[choice][1];
    id5 = flavors[choice][2];
  }
  else if(id1 == 21 && id2 == -4)
  {
    // Pick one of the flavor combinations (-6001, 3, -2)
    int flavors[1][3] = {{-6001, 3, -2}};
    vector<double> probs;
    double sum =;
    probs.push_back(/sum);
    int choice = rndmPtr->pick(probs);
    id3 = flavors[choice][0];
    id4 = flavors[choice][1];
    id5 = flavors[choice][2];
  }
  setId(id1, id2, id3, id4, id5);
  // Pick color flow
  int ncolor[1] = {1};
  if(id1 == 21 && id2 == -4 && id3 == -6001 && id4 == 3 && id5 == -2)
  {
    vector<double> probs;
    double sum = jamp2[0][0];
    for(int i = 0; i < ncolor[0]; i++ )
      probs.push_back(jamp2[0][i]/sum);
    int ic = rndmPtr->pick(probs);
    static int colors[1][10] = {{1, 3, 0, 1, 0, 3, 2, 0, 0, 2}};
    setColAcol(colors[ic][0], colors[ic][1], colors[ic][2], colors[ic][3],
        colors[ic][4], colors[ic][5], colors[ic][6], colors[ic][7],
        colors[ic][8], colors[ic][9]);
  }
  else if(id1 == -4 && id2 == 21 && id3 == -6001 && id4 == 1000011)
  {
    vector<double> probs;
    double sum = jamp2[0][0];
    for(int i = 0; i < ncolor[0]; i++ )
      probs.push_back(jamp2[0][i]/sum);
    int ic = rndmPtr->pick(probs);
    static int colors[1][10] = {{0, 1, 1, 3, 0, 3, 2, 0, 0, 2}};
    setColAcol(colors[ic][0], colors[ic][1], colors[ic][2], colors[ic][3],
        colors[ic][4], colors[ic][5], colors[ic][6], colors[ic][7],
        colors[ic][8], colors[ic][9]);
  }
}

The codes under the first "if" is unsuitable since I think it should be the same with the codes under "else if", which should be int flavors[1][3] = {{-6001, 3, -2}} instead of int flavors[1][3] = {{-6001, 1000011}}. However, new bugs are arosed since strange syntax
    double sum =;
    probs.push_back(/sum);
appeared and I think it should be
   double sum = matrix_element[1];
   probs.push_back(matrix_element[1]/sum);

[tang yilei (tangyilei10)]

void Sigma_vectorlike_gc_d4erp_erp_vlep::setIdColAcol()
{
  if(id1 == 4 && id2 == 21)
  {
    // Pick one of the flavor combinations
    int flavors[1][3] = {{7, -2000011}};
    vector<double> probs;
    double sum = matrix_element[12] + matrix_element[13] + matrix_element[14] +
        matrix_element[15] + matrix_element[16] + matrix_element[17];
    probs.push_back(matrix_element[12] + matrix_element[13] +
        matrix_element[14] + matrix_element[15] + matrix_element[16] +
        matrix_element[17]/sum);
    int choice = rndmPtr->pick(probs);
    id3 = flavors[choice][0];
    id4 = flavors[choice][1];
    id5 = flavors[choice][2];
  }
  else if(id1 == -4 && id2 == 21)
  {
    // Pick one of the flavor combinations
    int flavors[1][3] = {{-7, 2000011}};
    vector<double> probs;
    double sum = matrix_element[18] + matrix_element[19] + matrix_element[20] +
        matrix_element[21] + matrix_element[22] + matrix_element[23];
    probs.push_back(matrix_element[18] + matrix_element[19] +
        matrix_element[20] + matrix_element[21] + matrix_element[22] +
        matrix_element[23]/sum);
    int choice = rndmPtr->pick(probs);
    id3 = flavors[choice][0];
    id4 = flavors[choice][1];
    id5 = flavors[choice][2];
  }
  else if(id1 == 21 && id2 == 4)
  {
    // Pick one of the flavor combinations (7, 12, -11), (7, -16, -11), (7,
    // -12, -11), (7, 16, -11), (7, 14, -11), (7, -14, -11)
    int flavors[6][3] = {{7, 12, -11}, {7, -16, -11}, {7, -12, -11}, {7, 16,
        -11}, {7, 14, -11}, {7, -14, -11}};
    vector<double> probs;
    double sum = +++ ++;
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    int choice = rndmPtr->pick(probs);
    id3 = flavors[choice][0];
    id4 = flavors[choice][1];
    id5 = flavors[choice][2];
  }
  else if(id1 == 21 && id2 == -4)
  {
    // Pick one of the flavor combinations (-7, 12, 11), (-7, -16, 11), (-7,
    // -12, 11), (-7, 16, 11), (-7, 14, 11), (-7, -14, 11)
    int flavors[6][3] = {{-7, 12, 11}, {-7, -16, 11}, {-7, -12, 11}, {-7, 16,
        11}, {-7, 14, 11}, {-7, -14, 11}};
    vector<double> probs;
    double sum = +++ ++;
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    probs.push_back(/sum);
    int choice = rndmPtr->pick(probs);
    id3 = flavors[choice][0];
    id4 = flavors[choice][1];
    id5 = flavors[choice][2];
  }
  setId(id1, id2, id3, id4, id5);
  // Pick color flow
  int ncolor[12] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
  if(id1 == 21 && id2 == 4 && id3 == 7 && id4 == 12 && id5 == -11)
  {
    vector<double> probs;
    double sum = jamp2[0][0];
    for(int i = 0; i < ncolor[0]; i++ )
      probs.push_back(jamp2[0][i]/sum);
    int ic = rndmPtr->pick(probs);
    static int colors[1][10] = {{1, 2, 2, 0, 1, 0, 0, 0, 0, 0}};
    setColAcol(colors[ic][0], colors[ic][1], colors[ic][2], colors[ic][3],
        colors[ic][4], colors[ic][5], colors[ic][6], colors[ic][7],
        colors[ic][8], colors[ic][9]);
  }
  else if(id1 == 21 && i...

[tang yilei (tangyilei10)]

I just modify the export_cpp.py code and it gave the correct answer. The code is listed here:
    def getLegsWithDecayChains(self, proc):
        """Generate id lists with decay chains"""
        # decaying id
        decay_id = [d.get('legs')[0].get('id') for d in proc.get('decay_chains')]
        curr_final_id = [l.get('id') for l in proc.get('legs')
                        if l.get('state') and l.get('id') not in decay_id]
        # extend with the decay final state
        curr_final_id += [l.get('id') for dec in \
                        proc.get('decay_chains') for l in \
                        dec.get('legs') if l.get('state')]
        return curr_final_id

    def get_setIdColAcol_lines(self, color_amplitudes):
        """Generate lines to set final-state id and color info for process"""

        res_lines = []

        # Create a set with the pairs of incoming partons
        beams = set([(process.get('legs')[0].get('id'),
                      process.get('legs')[1].get('id')) \
                     for process in self.processes])

        #print beams
        # Now write a selection routine for final state ids
        for ibeam, beam_parts in enumerate(beams):
            #print beam_parts
            if ibeam == 0:
                res_lines.append("if(id1 == %d && id2 == %d){" % beam_parts)
            else:
                res_lines.append("else if(id1 == %d && id2 == %d){" % beam_parts)
            # Pick out all processes with this beam pair
            beam_processes = [(i, me) for (i, me) in \
                              enumerate(self.matrix_elements) if beam_parts in \
                              [(process.get('legs')[0].get('id'),
                                process.get('legs')[1].get('id')) \
                               for process in me.get('processes')]]
            # print len(beam_processes),"!!!"
            # Pick out all mirror processes for this beam pair
            beam_mirror_processes = []
            if beam_parts[0] != beam_parts[1]:
                beam_mirror_processes = [(i, me) for (i, me) in \
                              enumerate(self.matrix_elements) if beam_parts in \
                              [(process.get('legs')[1].get('id'),
                                process.get('legs')[0].get('id')) \
                               for process in me.get('processes')]]

            final_id_list = []
            final_mirror_id_list = []
            for (i, me) in beam_processes:
                valid_proc = [proc for proc in me.get('processes') \
                                      if beam_parts == \
                                      (proc.get('legs')[0].get('id'),
                                       proc.get('legs')[1].get('id'))]
                for proc in valid_proc:
                    # decaying id
                    decay_id = [d.get('legs')[0].get('id') for d in proc.get('decay_chains')]
                    curr_final_id = [l.get('id') for l in proc.get('legs')
                              if l.get('state') and l.get('id') not in decay_id]
                    # extend with the decay final state
                    curr_final_id += [l.get('id') fo...

[Olivier Mattelaer (olivier-mattelaer)]

Thanks,

Johan will look at this and will fix this,

Cheers,

Olivier

[Johan Alwall (johan-alwall)]

Hello Yilei,

Sorry for my poor review of Olivier's fix in v. 1.5.8, and thanks for forcing me to fix this once and for all.

I have fixed the bug in branch lp:~maddevelopers/madgraph5/1.5.9. Feel free to download this branch with
bzr lp:~maddevelopers/madgraph5/1.5.9
and make further tests. Please let us know if there are any further problems.

All the best,
Johan

[tang yilei (tangyilei10)]

Hi Johan
  I mischecked another related bug in the same function. In export_cpp.py, approximately line 1500( Because I modified the program directly so the line mark may be different from yours),
        for ime, me in enumerate(self.matrix_elements):
            if not me.get('has_mirror_process'):
                continue
            res_lines.append("else if(%s){" % \
                                 "||".join(["&&".join(["id%d == %d" % \
                                            (i+1, l.get('id')) for (i, l) in \
                                            enumerate(p.get_legs_with_decays())])\
                                           for p in me.get("legs")]))

should be modified as
        for ime, me in enumerate(self.matrix_elements):
            if not me.get('has_mirror_process'):
                continue
            res_lines.append("else if(%s){" % \
                                 "||".join(["&&".join(["id%d == %d" % \
                                            (i+1, l.get('id')) for (i, l) in \
                                            enumerate(p.get_legs_with_decays())])\
                                           for p in me.get_mirror_processes()])).

And I am puzzled about the way you wrote programs: see, you have both get_mirror_processes, get_legs_with_decays, and get('has_mirror_process') methods in processes class, why you just wrote a lot of lines to exchange the initial state to acquire mirror, combine legs with decay chains etc. in your previous program lines?

[Johan Alwall (johan-alwall)]

Hello Yilei,

I reverted Olivier's change and started from scratch, since it's much more elegant to simply use get_legs_with_decays. I also fixed the resonance treatment, since my previous treatment was not completely correct. Otherwise there are in fact no changes to the file (compared to the original version, before 1.5.8). Please use only the latest version, not your modified one.

All the best,
Johan

[tang yilei (tangyilei10)]

Oh, sorry, me.get_mirror_processes() seems do not exchange the initial legs and the program is still wrong. Sorry to report the wrong modification here. Is me.get_mirror_processes due to get the mirror without exchange the initial legs?

[Johan Alwall (johan-alwall)]

Dear Yilei,

Did you try the fixed version of the program in v. 1.5.9? Are there still problems in that version?

All the best,
Johan

[tang yilei (tangyilei10)]

Dear Johan,
  I'm sorry it's still cannot be run.
  The process I'm trying is such:

generate p p > t t~ , ( t > all all )
output pythia8

Then it returns:
INFO: Organizing processes into subprocess groups
INFO: Organizing processes into subprocess groups
INFO: Generating Helas calls for process: g g > t t~ WEIGHTED=2
INFO: Generating Helas calls for process: u u~ > t t~ WEIGHTED=2
INFO: Combined process c c~ > t t~ WEIGHTED=2 with process u u~ > t t~ WEIGHTED=2
INFO: Combined process d d~ > t t~ WEIGHTED=2 with process u u~ > t t~ WEIGHTED=2
INFO: Combined process s s~ > t t~ WEIGHTED=2 with process u u~ > t t~ WEIGHTED=2
INFO: Generating Helas calls for process: t > g t WEIGHTED=1
INFO: Combine g g > t t~ WEIGHTED=2 with decays t > g t WEIGHTED=1
INFO: Combine u u~ > t t~ WEIGHTED=2 with decays t > g t WEIGHTED=1
INFO: Processing color information for process: g g > t t~ WEIGHTED=2
  Decay: t > g t WEIGHTED=1
INFO: Processing color information for process: u u~ > t t~ WEIGHTED=2
  Decay: t > g t WEIGHTED=1
Command "output pythia8" interrupted with error:
IndexError : list index out of range
Please report this bug on https://bugs.launchpad.net/madgraph5
More information is found in 'MG5_debug'.
Please attach this file to your report.

The debug file shows:
#************************************************************
#* MadGraph 5 *
#* *
#* * * *
#* * * * * *
#* * * * * 5 * * * * *
#* * * * * *
#* * * *
#* *
#* *
#* VERSION 1.5.8 2013-05-05 *
#* *
#* The MadGraph Development Team - Please visit us at *
#* https://server06.fynu.ucl.ac.be/projects/madgraph *
#* *
#************************************************************
#* *
#* Command File for MadGraph 5 *
#* *
#* run as ./bin/mg5 filename *
#* *
#************************************************************
set group_subprocesses Auto
set ignore_six_quark_processes False
set gauge unitary
set complex_mass_scheme False
import model sm
define p = g u c d s u~ c~ d~ s~
define j = g u c d s u~ c~ d~ s~
define l+ = e+ mu+
define l- = e- mu-
define vl = ve vm vt
define vl~ = ve~ vm~ vt~
generate p p > t t~ , ( t > all all )
output pythia8
Traceback (most recent call last):
  File "/home/tangyilei/mgtry/madgraph/interface/extended_cmd.py", line 819, in onecmd
    return self.onecmd_orig(line, **opt)
  File "/home/tangyilei/mgtry/madgraph/interface/e...

[Johan Alwall (johan-alwall)]

Hi Yilei,

Thanks a lot for this, I had missed one line in the modification. Please download the branch again (you can use
bzr pull
from the directory you already downloaded) and try again.

Thanks,
Johan

[tang yilei (tangyilei10)]

Hi Johan
    I downloaded the new version and it runs perfect.
     Anyway, I have an extra suggestion which cannot be regarded as a bug.
    When I tried to calculate the generation of a particle together with decay chains, the output pythia8 function wrote lots of files due to the complicated decay process. e.g.
    A B > C D , ( C > all all )
    This will generate lots of diagrams which even seperate the final state u d into different Sigma_???_AB_CD_C_?? files. The consequence of these too many Sigma files resulted in stuck while pythia is running, because different events are seperated in different sigma classes, and took pythia to do a lot of useless trial. It took me one whole day to run a process with decays while a minute without decay chains. Can you supply some chosen items to combine some of these sigma classes while outputing in order to reduce the number of sigma files?

[Johan Alwall (johan-alwall)]

Hello Yilei,

Could you please give an example, so I can try it out and see what can be done?

Thanks,
Johan

[tang yilei (tangyilei10)]

Johan,
  OK, first I tried to output several 2->2 processes and the generated program look like this:
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4erp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4mulp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4murp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4ta1p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4ta2p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4se7c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4se8c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5elp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5erp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5mulp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5murp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5ta1p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5ta2p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5se7c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d5se8c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4elp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4erp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4mulp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4murp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4ta1p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4ta2p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4se7c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d4se8c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5elp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5erp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5mulp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5murp());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5ta1p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5ta2p());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5se7c());
  pythia.setSigmaPtr(new Sigma_vectorlike_gc_d5se8c());

  These are all processes which p p colide into sleptons and some kind of extra quarks.
  Then I defined a multiparticle called "small" which means "all of the standard model particles" and wanted to generate the sleptons decay chains directly and so I typed in:
  generate p p > slepton dquarks , ( slepton > small small )
  The program generated lots of classes. Then I tried to reduce the number, so I only tried
  generate g u > slepton dquarks , (slepton > small small )
  There're still lots of classes.

    pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_vlep());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_vlmup());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_vltap());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_udx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_usx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_ubx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_tdx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_tsx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_tbx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_cdx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_csx());
  pythia.setSigmaPtr(new Sigma_vectorlike_gu_d4elp_elp_cbx(...

[tang yilei (tangyilei10)]

The reason I tried to consider the 2->3 process including decay chains is actually because the "slepton" in my theory might not be on shell, so I cannot ignore the effects.

[Johan Alwall (johan-alwall)]

Hello Yilei,

I am pretty sure that the processes are already maximally combined. Pythia needs separate process classes for external legs with different mass. So it seems that you have not set light quark and lepton masses to zero. The easiest way to do this is to create a restriction card in your model directory. How this is done is described in some detail in the MadGraph 5 release paper, and you can look at the sm and mssm models for examples (restrict_xxx.dat). I suggest that you start from the param_card written by doing
python write_param_card.py
in your model directory.

All the best,
Johan

[tang yilei (tangyilei10)]

Thanks!

[tang yilei (tangyilei10)]

Oh, sorry. There should be another way to combine the processes and this is very important for your way may cause some problems. Notice that my processes contains lots of processes sharing the similar final states however different intermediate states. For e.g.

gu_d4se7c_se7c_udx
gu_d4se8c_se8c_udx.

Strictly these two process can't be seperated into different classes because they may actually interfere. Can you supply some function to combine processes like this?

[Johan Alwall (johan-alwall)]

Hello Yilei,

Why not simply generate
p p > se7~ | se8~ > d4 j j
(please see "help generate").

Using the > > syntax, you can fine-tune exactly which intermediate s-channel propagators you want to include. However, it still doesn't look like you have set your quark masses to zero. Note that this is crucial to reduce the number of subprocesses. If you have some new particles that have identical masses and widths (typically this is the case for 1st and 2nd generation squarks in unified MSSM models), you should identify those as well in your restriction card, to further reduce the number of processes.

All the best,
Johan

[tang yilei (tangyilei10)]

Hi Johan,
  I tried, and it worked. The reason I ignored this method is due to your mannual describe it as "Alternative required s-channels can be separated by "|":". That misled me into the wrong concept that this "|" symbol can only indicate that the following particles are all decayed by the intermediate state. That is due to the usual "s-channel" definition concept I learnt from general book:
\ /
 -----
/ \
can only be regared as s-channel while

\ /
  |
/ \___
    \
      \
can only be regarded as t-channel together with something decays. Maybe you would change the way you describe in mannual?