Randomly pick one color for each event written to file (jamp2 handling)

[valassi]

Randomly pick one color for each event written to file (jamp2 handling)

Will give details elsewhere

[valassi]

Full details are given in #404.

This is on the critical path to get unweighted event sets to the experiments where each event has one color index (but it should not be needed for cross section calculations)

[valassi]

As I mentioned in oliviermattelaer/mg5amc_test#14 (comment), I would like to go about this issue in two steps:

• First, use a 'dummy' algorithm with a fake random choice of color, where each event gets assigned a (valid) color. In practice, JAMP2 must be filled with something that allows the addmothers algorithm to proceed, otherwise one gets the error 'failed to reduce to color indices' in "Error: failed to reduce to color indices" in gg to ttg (only for event generation, not for cross sections) oliviermattelaer/mg5amc_test#14
• In a second step, ensure a correct random choice of color. One option is to export all relevant JAMP2 through the bridge and use the Fortran addmothers algorithm as it is. Another option is to implement the random choice in cudacpp (this is the option that we said we would target, and @oliviermattelaer I think you are working on it...?).

[valassi]

About the design of the dummy algorithm, I am having a look at addmothers.f

• JAMP2 has dimensions (0:MAXFLOW, NB_PAGE_MAX): essentially, number of colors (plus 1!) and number of events per iteration

• The variable MAXFLOW is essentially NCOLOR? See discussion in JAMP2 dimensions - difference between MAXFLOW and NCOLOR? oliviermattelaer/mg5amc_test#21. Olivier mentioned that MAXFLOW may be bigger than NCOLOR. However, I have the impression that NCOLOR is really the relevant variable for the algorithms, ie the one used in loops, so I will concentrate on that and not on MAXFLOW.

• In particular, the variable nc in addmothers.f is taken from JAMP(0,IVEC) for event IVEC. This is initialized in the following way in code generation:

madgraph/iolibs/template_files/matrix_madevent_v4.inc:        JAMP2(0)=%(ncolor)d
madgraph/iolibs/template_files/matrix_madevent_v4.inc:    PARAMETER (NWAVEFUNCS=%(nwavefuncs)d, NCOLOR=%(ncolor)d) 

In my Fortran patch, I will therefore for the moment assume that JAMP2(0) is equal to parameter NCOLOR. (Suggestion: if this is the case, it would be better to specify this in the Fortran code, rather than relying on code generation to fill the same numerical value).

• The array ICOLAMP seems to be important in defining which colors are "valid" and which ones are not

      nc = int(jamp2(0, ivec))
      is_LC = .true.
      maxcolor=0
      if(nc.gt.0)then
      if(icolamp(1,lconfig,iproc)) then
        targetamp(1)=jamp2(1,ivec)
c        print *,'Color flow 1 allowed for config ',lconfig
      else
        targetamp(1)=0d0
      endif
      do ic =2,nc
        if(icolamp(ic,lconfig,iproc))then
          targetamp(ic) = jamp2(ic,ivec)+targetamp(ic-1)
c          print *,'Color flow ',ic,' allowed for config ',lconfig,targetamp(ic)
        else
          targetamp(ic)=targetamp(ic-1)
        endif
      enddo

• The targetamp here are essentially the cumulative probabilities of the random choice algorithm. Apparently their normalization is not important.

      xtarget=ran1(iseed)*targetamp(nc)

      ic = 1
      do while (targetamp(ic) .lt. xtarget .and. ic .lt. nc)
         ic=ic+1
      enddo
      if(targetamp(nc).eq.0) ic=0 

I should therefore be able to realtively easily design a dummy algorithm...

[valassi]

I have implemented the first step, a 'dummy' choice of color, in #489. This provides the workarund for oliviermattelaer/mg5amc_test#14

Note that in the end addmothers.f has its own way to set a dummy color, as long as JAMP2(0) is correctly set to NCOLORS, which is the only thing I have done.

The next step, doing a real random choice of color, will be more complex...

[oliviermattelaer]

note that oliviermattelaer/mg5amc_test#24
is changing the fortran handling of color/helicity in order to have a clean api
This means:

• the fortran/cpp api is now changed to the what we discussed a while ago (input: 2 random number, output: 2 integer)
• New fortran function has been added to the code
  • select_color(rcol, jamp2, channel, proc_id, icol [out])
    • need to write such function within the CPP plugin (work should have started in that direction before the break)
  • get_helicities(..)
    • no need to convert such function to C++

[valassi]

This is in progress. I have some agorithm but it gives me different results from Fortran and C++ (the last two columns).

Maybe not ideal to use gg_tt for debugging as it has only two colors...

  33  0.66637210E-01  0.66637210E-01          1.00000000000  3  3   1   1
  34  0.36652604E-01  0.36652604E-01          1.00000000000  8  8   1   1
  35  0.41960814E-01  0.41960814E-01          1.00000000000 12 12   1   1
  36  0.64245071E-01  0.64245071E-01          1.00000000000  2  2   1   1
  37  0.67260622E-01  0.67260622E-01          1.00000000000 15 15   1   1
  38  0.72644897E-01  0.72644897E-01          1.00000000000 15 15   1   1
  39  0.48924468E-01  0.48924468E-01          1.00000000000 15 15   1   1
  40  0.71085828E-01  0.71085828E-01          1.00000000000  2  2   1   1
  41  0.71529802E-01  0.71529802E-01          1.00000000000 15 15   1   1
  42  0.70385722E-01  0.70385722E-01          1.00000000000  2  2   1   1
  43  0.31575040E-01  0.31575040E-01          1.00000000000 12 12   1   1
  44  0.69272519E-01  0.69272519E-01          1.00000000000 15 15   2   1
  45  0.71009335E-01  0.71009335E-01          1.00000000000 15 15   1   1
  46  0.68831440E-01  0.68831440E-01          1.00000000000  2  2   1   1
  47  0.69539468E-01  0.69539468E-01          1.00000000000  2  2   1   1
  48  0.70394921E-01  0.70394921E-01          1.00000000000 15 15   1   1
  49  0.73183046E-01  0.73183046E-01          1.00000000000  3  3   1   1
  50  0.60352285E-01  0.60352285E-01          1.00000000000 15 15   1   1
  51  0.18336225E-01  0.18336225E-01          1.00000000000 15 15   1   1
  52  0.47909190E-01  0.47909190E-01          1.00000000000  2  2   1   1
  53  0.12439986E-01  0.12439986E-01          1.00000000000 15 15   1   1
  54  0.66421684E-01  0.66421684E-01          1.00000000000 15 15   1   1
  55  0.68700645E-01  0.68700645E-01          1.00000000000  2  2   1   1
  56  0.66674348E-01  0.66674348E-01          1.00000000000 12 12   1   1
  57  0.34727484E-01  0.34727484E-01          1.00000000000 15 15   1   2
  58  0.71093852E-01  0.71093852E-01          1.00000000000  2  2   1   1
  59  0.71165949E-01  0.71165949E-01          1.00000000000  2  2   1   1
  60  0.70806691E-01  0.70806691E-01          1.00000000000 15 15   1   1
  61  0.55143515E-01  0.55143515E-01          1.00000000000 15 15   1   1
  62  0.59417397E-01  0.59417397E-01          1.00000000000 15 15   1   1
  63  0.42014188E-01  0.42014188E-01          1.00000000000  5  5   1   1
  64  0.70783989E-01  0.70783989E-01          1.00000000000 15 15   2   2

[valassi]

There was a silly bug (jamp2_sv not initialised at every SIMD page) in C++. Actually CUDA was already working.

Now everything looks good. This is again from BothDebug, comparing MEs, hel and col:

   1  0.64067963E-01  0.64067963E-01          1.00000000000 12 12   1   1
   2  0.58379673E-01  0.58379673E-01          1.00000000000  9  9   1   1
   3  0.70810768E-01  0.70810768E-01          1.00000000000 15 15   1   1
   4  0.67192668E-01  0.67192668E-01          1.00000000000 15 15   1   1
   5  0.71590585E-01  0.71590585E-01          1.00000000000 15 15   1   1
   6  0.72862110E-01  0.72862110E-01          1.00000000000  2  2   1   1
   7  0.14271254E-01  0.14271254E-01          1.00000000000  5  5   1   1
   8  0.63986754E-01  0.63986754E-01          1.00000000000 14 14   2   2
   9  0.46316382E-01  0.46316382E-01          1.00000000000 12 12   1   1
  10  0.35372741E-01  0.35372741E-01          1.00000000000  2  2   1   1
  11  0.73958407E-01  0.73958407E-01          1.00000000000  5  5   1   1
  12  0.70691203E-01  0.70691203E-01          1.00000000000 15 15   1   1
  13  0.70805000E-01  0.70805000E-01          1.00000000000  2  2   1   1
  14  0.30801404E-01  0.30801404E-01          1.00000000000  2  2   1   1
  15  0.64111868E-01  0.64111868E-01          1.00000000000 10 10   2   2
  16  0.74312047E-01  0.74312047E-01          1.00000000000  2  2   1   1
  17  0.60961835E-01  0.60961835E-01          1.00000000000  2  2   1   1
  18  0.67698020E-01  0.67698020E-01          1.00000000000  2  2   1   1
  19  0.49748773E-01  0.49748773E-01          1.00000000000  2  2   1   1
  20  0.71951996E-01  0.71951996E-01          1.00000000000  2  2   1   1
  21  0.52116331E-01  0.52116331E-01          1.00000000000  5  5   1   1
  22  0.69245648E-01  0.69245648E-01          1.00000000000 15 15   1   1
  23  0.64808141E-01  0.64808141E-01          1.00000000000 15 15   1   1
  24  0.66861231E-01  0.66861231E-01          1.00000000000  3  3   2   2
  25  0.70041112E-01  0.70041112E-01          1.00000000000  2  2   1   1
  26  0.61135249E-01  0.61135249E-01          1.00000000000 15 15   1   1
  27  0.66574932E-01  0.66574932E-01          1.00000000000 15 15   1   1
  28  0.67312068E-01  0.67312068E-01          1.00000000000  3  3   1   1
  29  0.47056643E-01  0.47056643E-01          1.00000000000  5  5   2   2
  30  0.70509435E-01  0.70509435E-01          1.00000000000  2  2   2   2
  31  0.23138767E-01  0.23138767E-01          1.00000000000 15 15   1   1
  32  0.76096234E-01  0.76096234E-01          1.00000000000 15 15   1   1

This is implemented in MR #573 (currently WIP). When that is merged, this issue will be closed.