madweight stop pair production zero weight

Hi

I am trying to use MadWeight for stop pair production hypothesis. Here is summary of what I noticed.

- bin/internal/madweight/Cards.py had "give_mass_width_dict()" had hardcoded SM mass and width information and didn't contain 1000006 (stop) for now I hard coded and added 1000006 : "MSU3" and 1000006 : "WSU3". That helped push through to prep/compilation steps and started comp_madweight with success

- But I got all 0.0 for weight calculations of all events

- I tracked the problem down to main_code subroutine where the first block calculation yields negative jacobian factor and basically returns null

- The specific block that keeps returning negative jacobian is the block_b where the block_b is applied to t->w+ b -> e+ ve b chain.

- The particle mgid assign seems to be okay... not sure why it's giving me negative jacobian. I haven't further figure out where in block_b routine that seems to be odd so I am posting a help here.

Below I post a whole bunch of additional information for complete information
I hope this is enough!

1. bazaar revision number
2. proc_card
3. param_card
4. one of the lhco event I used (all events gave me zero weight so this should be enough)
5. TF information (i.e. used all_delta, tried with dbl_gauss still gave me zero)
6. MadWeight_card
7. Details on how I changed the "bin/internal/madweight/Cards.py"
8. What the class/block assignment was for my process after generation
9. Information about where in main_code subroutine I *always* get negative jacobian

Thanks!!!
Philip

==========================================================================
I bzr branched checked out ~2 weeks ago. So I think the revno is,

------------------------------------------------------------
revno: 285
committer: olivier Mattelaer <email address hidden>
branch nick: madweight
timestamp: Thu 2014-02-20 15:11:07 +0000
message:
  fix the unittest for the cuts
------------------------------------------------------------

==========================================================================
Here is the proc_card content

set group_subprocesses Auto
set ignore_six_quark_processes False
set gauge unitary
set complex_mass_scheme False
import model mssm
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~
define s1 ul ul~ ur ur~ dl dl~ dr dr~ sl sl~ sr sr~ cl cl~ cr cr~
define s3 b1 b1~ b2 b2~ t2 t2~
define sch x1+ x2+ x1- x2-
generate p p > t1 t1~ / sch s1 s3, t1 > b e+ ve n1 / sch s1 s3, t1~ > b~ e- ve~ n1 / sch s1 s3 @1
output madweight mw_stoppair_v3

==========================================================================
here is the param_card_1.dat

######################################################################
## PARAM_CARD AUTOMATICALY GENERATED BY MG5 ####
######################################################################
###################################
## INFORMATION FOR HMIX
###################################
Block UMIX #
     1 1 1.00000000E+00
     1 2 0.00000000E+00
     2 1 0.00000000E+00
     2 2 1.00000000E+00
Block AE Q= 1.000000e+00 #
     1 1 0.000000e+00
     3 3 -2.517769e+02
     2 2 0.000000e+00
Block AD Q= 1.000000e+00 #
     1 1 0.000000e+00
     3 3 -7.972741e+02
     2 2 0.000000e+00
Block AU Q= 1.000000e+00 #
     1 1 0.000000e+00
     3 3 -4.981299e+02
     2 2 0.000000e+00
Block NMIX #
     1 1 0.98
     1 3 0.00
     1 2 -0.20
     1 4 0.00
     3 1 0.00
     3 3 1.00
     3 2 0.00
     3 4 0.00
     2 1 0.20
     2 3 0.00
     2 2 0.98
     2 4 0.00
     4 1 0.00
     4 3 0.00
     4 2 0.00
     4 4 1.00
Block YE Q= 1.000000e+00 #
     3 3 1.008908e-01
Block SBOTMIX #
     1 1 9.387379e-01
     1 2 3.446319e-01
     2 1 -3.446319e-01
     2 2 9.387379e-01
Block MSOFT Q= 1.000000e+00 #
     48 5.231488e+02
     21 3.233749e+04
     22 -1.288001e+05
     49 5.198673e+02
     33 1.944960e+02
     32 1.953348e+02
     31 1.953348e+02
     45 5.295112e+02
     42 5.475735e+02
     36 1.340434e+02
     35 1.364941e+02
     34 1.364941e+02
     1 1.013965e+02
     46 4.232459e+02
     3 5.882630e+02
     2 1.915042e+02
     47 5.231488e+02
     43 4.987639e+02
     41 5.475735e+02
     44 5.295112e+02
Block YD Q= 1.000000e+00 #
     3 3 1.388402e-01
Block VMIX #
     1 1 1.00000000E+00
     1 2 0.00000000E+00
     2 1 0.00000000E+00
     2 2 1.00000000E+00
Block MASS #
     2000013 1.00000e+05
     22 0.000000e+00
     24 7.982901e+01
     25 1.108991e+02
     2000011 1.00000e+05
     1000021 1.00000e+05
     1000022 0.01000e+02
     23 9.118760e+01
     1000006 2.30000e+02
     1000004 1.00000e+05
     1000005 1.00000e+05
     1000002 1.00000e+05
     1000003 1.00000e+05
     1000001 1.00000e+05
     1 0.000000e+00
     3 0.000000e+00
     2 0.000000e+00
     5 4.889917e+00
     4 0.000000e+00
     6 175
     1000037 1.00000e+05
     2000015 1.00000e+05
     2000005 1.00000e+05
     1000024 1.00000e+05
     1000025 1.00000e+05
     1000035 1.00000e+05
     2000001 1.00000e+05
     2000006 1.00000e+05
     2000002 1.00000e+05
     1000015 1.00000e+05
     13 0.000000e+00
     14 0.000000e+00
     1000011 1.00000e+05
     2000004 1.00000e+05
     1000013 1.00000e+05
     1000012 1.00000e+05
     15 1.777000e+00
     1000014 1.00000e+05
     2000003 1.00000e+05
     1000016 1.00000e+05
     1000023 1.00000e+05
     37 4.078790e+02
     36 3.995839e+02
     35 3.999601e+02
     16 0.000000e+00
     21 0.000000e+00
     11 0.000000e+00
     12 0.000000e+00
Block ALPHA #
          -1.138252e-01
Block QNUMBERS 2000005 # b2
     1 -1
     3 3
     2 1
     4 1
Block MODSEL #
     1 1
Block SMINPUTS #
     1 1.279340e+02
     3 1.180000e-01
     2 1.166370e-05
     4 9.118760e+01
     7 1.777000e+00
     6 1.750000e+02
Block YU Q= 1.000000e+00 #
     3 3 8.928445e-01
Block STOPMIX #
     1 1 1.00000000E+00
     1 2 0.00000000E+00
     2 1 0.00000000E+00
     2 2 1.00000000E+00
Block HMIX Q= 1.000000e+00 #
     1 3.576810e+02
     3 2.449353e+02
     2 9.748624e+00
     4 1.664391e+05
Block STAUMIX #
     1 1 2.824872e-01
     1 2 9.592711e-01
     2 1 9.592711e-01
     2 2 -2.824872e-01
DECAY 2000013 2.161216e-01
DECAY 15 0.000000e+00
DECAY 24 2.002822e+00
DECAY 25 1.986108e-03
DECAY 2000011 2.161216e-01
DECAY 1000021 5.506754e+00
DECAY 22 0.000000e+00
DECAY 23 2.411433e+00
DECAY 1000006 2.021596e+00
DECAY 1000004 5.477195e+00
DECAY 1000005 3.736276e+00
DECAY 1000002 5.477195e+00
DECAY 1000003 5.312788e+00
DECAY 1000001 5.312788e+00
DECAY 1 0.000000e+00
DECAY 3 0.000000e+00
DECAY 2 0.000000e+00
DECAY 5 0.000000e+00
DECAY 4 0.000000e+00
DECAY 6 1.561950e+00
DECAY 1000037 2.486895e+00
DECAY 2000015 2.699061e-01
DECAY 2000005 8.015663e-01
DECAY 1000024 1.704145e-02
DECAY 1000025 1.915985e+00
DECAY 1000035 2.585851e+00
DECAY 2000006 7.373133e+00
DECAY 2000002 1.152973e+00
DECAY 2000001 2.858123e-01
DECAY 13 0.000000e+00
DECAY 14 0.000000e+00
DECAY 1000011 2.136822e-01
DECAY 2000004 1.152973e+00
DECAY 1000013 2.136822e-01
DECAY 1000012 1.498816e-01
DECAY 1000015 1.483273e-01
DECAY 1000014 1.498816e-01
DECAY 2000003 2.858123e-01
DECAY 1000016 1.475190e-01
DECAY 1000022 0.000000e+00
DECAY 1000023 2.077700e-02
DECAY 37 5.469628e-01
DECAY 36 6.321785e-01
DECAY 35 5.748014e-01
DECAY 16 0.000000e+00
DECAY 21 0.000000e+00
DECAY 11 0.000000e+00
DECAY 12 0.000000e+00

==========================================================================
I used all_delta TF

==========================================================================
Here is one LHCO event I used

# typ eta phi pt jmass ntrk btag had/em dummy dummy
0 3 5
1 1 -2.40724711271 4.66537645403 23.2490198061 0.0 1 0 0 0 0
2 4 -0.971694479806 5.19220406648 55.7142783172 4.6999998093 1 2 0 0 0
3 1 0.302029871753 4.30168137228 67.6627564632 0.0 -1 0 0 0 0
4 4 0.101101483177 0.951867778546 73.1067311195 4.6999998093 1 2 0 0 0
5 6 0 2.06025553741 85.1323449643 0.0 1 0 0 0 0

==========================================================================
Here is my MadWeight_card.dat

Block MW_Run
# TAG VALUE UTILITY
     name fermi # name for the run
     nb_exp_events 4 # number of experimental events to consider
     MW_int_points 2000 # number of points (by permutation) in MadWeight integration for survey
     MW_int_refine 10000 # number of points (by permutation) in MadWeight integration for refine
     precision 0.005 # stops computation if precision is reached.
     nb_event_by_node 1 # one job submission compute the weight for N events
     log_level weight # from low level of log to extensive log:
                                  # weight, permutation, channel, full
     use_cut F # use the cut defined in run_card.dat
     bw_cut F # use the BW cut
     nwa 0.1 # width below narrow width approximation is used.
     isr 0 # isr=0 : ignore ISR effect (except if all FS particles are visible)
                                  # isr=1 : correct kinematic based on reconstructed Pt(isr)
     inputfile './Events/input.lhco' # path to the input file (in lhco format)

#*************************************************************************
## define the different param_card's ##
#*************************************************************************
Block MW_parameter
# TAG VALUE UTILITY
    mode 1 # type of input
                           # 0 : inputs are read from the cards: param_card_1.dat, param_card_2.dat,...
                           # 1 : redefines some values from param_card.dat according to the form below
                           # 2 : same but the value for different parameters are modified simultaneously
#
# # first parameter #
     11 mass # Block of the parameter to change
     12 6 # id of the parameter to change
     13 175 # here you can enter the different values:
#
# # second parameter #
# !DELETE ALL THE FOLLOWING TAG IF YOU WANT TO RUN WITH ONLY ONE PARAM!
# 21 MGCKM # Block of the parameter to change
# 22 1 # id of the paramter to change
# 22 2 # id2 of the paramter to change
# 23 1.5E-02 # here you can enter the different values:
# 23 1.8E-02 # add a new line with tag 23 to introduce a new value
#
# use same syntax for parameters 3,4,...
#*************************************************************************
## Permutations ##
#*************************************************************************
Block MW_perm
# TAG VALUE UTILITY
     permutation T # make permutation
     bjet_is_jet T # consider permutation between b-jets and light jets
     montecarlo T # Monte-Carlo over permutation (Huge speed up if many permutation)
     preselect 'default' # How to pre-select the correct permutation set.
                            # put 'None' if no pre-selection to perform.
                            # You can set the path to a fortran file defining the require function
                            # See file SubProcesses/permutation_weight_default.dat for
                            # instructions.
     min_perm_cut 5e-4 # Cut for discarding permutation on the preselected method
#*************************************************************************
## Phase-Space Integration mapping ##
#*************************************************************************
Block MW_gen
    force_nwa 2 # Only consider the change of variable alligning particles
                            # with width smaller than this value. This speed up the code
                            # but can lead to zero weight for background event where the
                            # kinematic doesn't agree with the associated mass.
                            # if "mw_run nwa" parameter is bigger than this value, that
                            # value is used for this parameter automatically.

==========================================================================
The line I changed in the Cards.py looks like this now,

    #3 #########################################################################
    def give_mass_width_dict(self):
        """ return two dict {pid:fortran_mass_code} and {pid:fortran_width_code}
            pid value are always positive
        """
# if not self.charged:
# self.read() #not automaticly read for the moment
# dict_mass={}
# dict_width={}
# for data in self.info['particles']:
# pid=abs(int(data[-1]))
# mass=data[-5]
# width=data[-4]

# dict_mass[pid]=mass
# dict_width[pid]=width

        return {1: 'ZERO', 2: 'ZERO', 3: 'ZERO', 4: 'ZERO', 5: 'MB', 6: 'MT', 11: 'ZERO', 12: 'ZERO', 13: 'ZERO', 14: 'ZERO', 15: 'MTA', 16: 'ZERO', 21: 'ZERO', 22: 'ZERO', 23: 'MZ', 24: 'MW', 25: 'MH', 1000006: 'MSU3'}, {1: 'ZERO', 2: 'ZERO', 3: 'ZERO', 4: 'ZERO', 5: 'ZERO', 6: 'WT', 11: 'ZERO', 12: 'ZERO', 13: 'ZERO', 14: 'ZERO', 15: 'ZERO', 16: 'ZERO', 21: 'ZERO', 22: 'WA', 23: 'WZ', 24: 'WW', 25: 'WH', 1000006: 'WSU3'}

==========================================================================
when I ran madweight.py I got the following block/class assignments

treating P1_gg_t1t1x_t1_bepven1_t1x_bxemvexn1 directory
structure of the configuration 1:
particle: 6 pid : 1000022 level: 1 mother: -6 twin: -5
particle: 10 pid : 1000022 level: 1 mother: -3 twin: -2
particle: 3 pid : 5 level: 2 mother: -5 twin: -4
particle: 7 pid : -5 level: 2 mother: -2 twin: -1
particle: 4 pid : -11 level: 3 mother: -4 twin: 5
particle: 5 pid : 12 level: 3 mother: -4 twin: 4
particle: 8 pid : 11 level: 3 mother: -1 twin: 9
particle: 9 pid : -12 level: 3 mother: -1 twin: 8
particle: -1 pid : -24 level: 3 channel: S des: 9 8 mother: -2 mass/width: 79.82901/2.002822
particle: -2 pid : -6 level: 2 channel: S des: -1 7 mother: -3 mass/width: 175.0/1.56195
particle: -3 pid : -1000006 level: 1 channel: S des: 10 -2 mother: -7 mass/width: 230.0/2.021596
particle: -4 pid : 24 level: 3 channel: S des: 5 4 mother: -5 mass/width: 79.82901/2.002822
particle: -5 pid : 6 level: 2 channel: S des: -4 3 mother: -6 mass/width: 175.0/1.56195
particle: -6 pid : 1000006 level: 1 channel: S des: 6 -5 mother: -7 mass/width: 230.0/2.021596
particle: -7 pid : 21 level: 0 channel: T des: -3 -6 mass/width: 0.0/0.0
3 ECS('s) 7 propagator(s) 4 missing particles(s)
detail :
b(2) b(2) g(2)
4 blob(s) associated
Blob details: main -6
blob generation: [B|5 4 3 :-5 ] [C|6 -5 :-6 ]
Blob details: main -5
blob generation: [B|5 4 3 :-5 ]
Blob details: main -3
blob generation: [B|9 8 7 :-2 ] [C|10 -2 :-3 ]
Blob details: main -2
blob generation: [B|9 8 7 :-2 ]

==========================================================================
Three positions where it returned negative jacobians are pointed by "***HERE***"

C+-----------------------------------------------------------------------+
C| MAIN CODE FOR MADWEIGHT |
C| |
C| Author: Pierre Artoisenet (UCL-CP3) |
C| Olivier Mattelaer (UCL-CP3) |
C+-----------------------------------------------------------------------+
C| This file is generated automaticly by MADWEIGHT-ANALYZER |
C+-----------------------------------------------------------------------+
      subroutine main_code(x,n_var)
C+-----------------------------------------------------------------------+
C| Central Routine for the change of variable choice |
C| - num_sol: number of the solution to charge |
C| - x : random number from vegas |
C+-----------------------------------------------------------------------+
 implicit none

 double precision x(20)
 integer n_var
 double precision S,X1,X2,PSWGT,JAC
        common /PHASESPACE/ S,X1,X2,PSWGT,JAC
c
 double precision multi_channel_weight
 external multi_channel_weight
c
 include 'd_choices.inc'
 include 'phasespace.inc'
 include 'data.inc'
C+-----------------------------------------------------------------------+
C| Scedullar Part |
C+-----------------------------------------------------------------------+

       if (config_pos.eq.1) then

C+-----------------------------------------------------------------------+
C| |
C| ** Enlarged Contraint Sector global information ** |
C| |
C| Class: B |
C| particle in ECS : 6(missing) -5(propagator) |
C| blob linked are generated by :-5 -3 |
C| |
C| |
C+-----------------------------------------------------------------------+

        call block_b(x,5,4,3,-4,-5)***HERE***
        if (jac.le.0d0) return
        call block_b(x,9,8,7,-1,-2)
        if (jac.le.0d0) return
C ++++++++++++
        call block_c(x,10,-2,-3)
        if (jac.le.0d0) return
        call class_b(x,6,-5,-6)
        if (jac.le.0d0) return

        jac=jac*multi_channel_weight(1)

       elseif (config_pos.eq.2) then
C+-----------------------------------------------------------------------+
C| |
C| ** Enlarged Contraint Sector global information ** |
C| |
C| Class: B |
C| particle in ECS : 10(missing) -2(propagator) |
C| blob linked are generated by :-2 -6 |
C| |
C| |
C+-----------------------------------------------------------------------+

        call block_b(x,9,8,7,-1,-2)***HERE***
        if (jac.le.0d0) return
        call block_b(x,5,4,3,-4,-5)
        if (jac.le.0d0) return
C ++++++++++++
        call block_c(x,6,-5,-6)
        if (jac.le.0d0) return
        call class_b(x,10,-2,-3)
        if (jac.le.0d0) return

        jac=jac*multi_channel_weight(2)

       elseif (config_pos.eq.3) then
C+-----------------------------------------------------------------------+
C| |
C| ** Enlarged Contraint Sector global information ** |
C| |
C| Class: G |
C| particle in ECS : 6(missing) 10(missing) -5(propagator) |
C| -2(propagator) |
C| blob linked are generated by :-5 -2 |
C| |
C| |
C+-----------------------------------------------------------------------+

        call block_b(x,5,4,3,-4,-5)***HERE***
        if (jac.le.0d0) return
        call block_b(x,9,8,7,-1,-2)
        if (jac.le.0d0) return
        call class_g(x,6,10,-5,-2,-6,-3)
        if (jac.le.0d0) return

        jac=jac*multi_channel_weight(3)

        endif
        return
        end