1
=== modified file 'MadSpin/decay.py'
2
--- MadSpin/decay.py	2018-04-25 06:47:06 +0000
3
+++ MadSpin/decay.py	2018-12-20 08:18:55 +0000
4
@@ -3963,14 +3963,23 @@
5
3963
                elif mode=='full':
3963
                elif mode=='full':
6
3964
                    stdin_text="5 0 0 0 \n"  # before closing, write down the seed 
3964
                    stdin_text="5 0 0 0 \n"  # before closing, write down the seed 
7
3965
                    external = self.calculator[('full',path)]
3965
                    external = self.calculator[('full',path)]
9
3966
                    external.stdin.write(stdin_text)
3966
                    try:
10
3967
                        external.stdin.write(stdin_text)
11
3968
                    except Exception:
12
3969
                        continue
13
3967
                    ranmar_state=external.stdout.readline()
3970
                    ranmar_state=external.stdout.readline()
14
3968
                    ranmar_file=pjoin(path,'ranmar_state.dat')
3971
                    ranmar_file=pjoin(path,'ranmar_state.dat')
15
3969
                    ranmar=open(ranmar_file, 'w')
3972
                    ranmar=open(ranmar_file, 'w')
16
3970
                    ranmar.write(ranmar_state)
3973
                    ranmar.write(ranmar_state)
17
3971
                    ranmar.close()
3974
                    ranmar.close()
20
3972
                    external.stdin.close()
3975
                    try:
21
3973
                    external.stdout.close()
3976
                        external.stdin.close()
22
3977
                    except Exception:
23
3978
                        continue
24
3979
                    try:
25
3980
                        external.stdout.close()
26
3981
                    except Exception:
27
3982
                        continue
28
3974
                    external.terminate()
3983
                    external.terminate()
29
3975
                    del external
3984
                    del external
30
3976
        else:
3985
        else:
31
3977
3986
32
=== modified file 'MadSpin/interface_madspin.py'
33
--- MadSpin/interface_madspin.py	2018-03-11 12:10:19 +0000
34
+++ MadSpin/interface_madspin.py	2018-12-20 08:18:55 +0000
35
@@ -35,6 +35,7 @@
36
35
import madgraph.interface.madgraph_interface as mg_interface
35
import madgraph.interface.madgraph_interface as mg_interface
37
36
import madgraph.interface.master_interface as master_interface
36
import madgraph.interface.master_interface as master_interface
38
37
import madgraph.interface.madevent_interface as madevent_interface
37
import madgraph.interface.madevent_interface as madevent_interface
39
38
import madgraph.interface.common_run_interface as common_run_interface
40
38
import madgraph.interface.reweight_interface as rwgt_interface
39
import madgraph.interface.reweight_interface as rwgt_interface
41
39
import madgraph.various.misc as misc
40
import madgraph.various.misc as misc
42
40
import madgraph.iolibs.files as files
41
import madgraph.iolibs.files as files
43
@@ -53,6 +54,68 @@
44
53
cmd_logger = logging.getLogger('cmdprint2') # -> print
54
cmd_logger = logging.getLogger('cmdprint2') # -> print
45
54
55
46
55
56
47
57
class MadSpinOptions(banner.ConfigFile):
48
58
    
49
59
    def default_setup(self):
50
60
51
61
        self.add_param("max_weight", -1)
52
62
        self.add_param('curr_dir', os.path.realpath(os.getcwd()))
53
63
        self.add_param('Nevents_for_max_weigth', 0)
54
64
        self.add_param("max_weight_ps_point", 400)
55
65
        self.add_param('BW_cut', -1)
56
66
        self.add_param('nb_sigma', 0.)
57
67
        self.add_param('ms_dir', '')
58
68
        self.add_param('max_running_process', 100)
59
69
        self.add_param('onlyhelicity', False)
60
70
        self.add_param('spinmode', "madspin", allowed=['madspin','none','onshell'])
61
71
        self.add_param('use_old_dir', False, comment='should be use only for faster debugging')
62
72
        self.add_param('run_card', '' , comment='define cut for spinmode==none. Path to run_card to use')
63
73
        self.add_param('fixed_order', False, comment='to activate fixed order handling of counter-event')
64
74
        self.add_param('seed', 0, comment='control the seed of madspin')
65
75
        self.add_param('cross_section', {'__type__':0.}, comment="forcing normalization of cross-section after MS (for none/onshell)" )
66
76
        self.add_param('new_wgt', 'cross-section' ,allowed=['cross-section', 'BR'], comment="if not consistent number of particles, choose what to do for the weight. (BR: means local according to number of part, cross use the force cross-section")
67
77
        self.add_param('input_format', 'auto', allowed=['auto','lhe', 'hepmc', 'lhe_no_banner'])
68
78
        
69
79
        
70
80
    ############################################################################
71
81
    ##  Special post-processing of the options                                ## 
72
82
    ############################################################################
73
83
    def post_set_ms_dir(self, value, change_userdefine, raiseerror):
74
84
        """ special handling for set ms_dir """
75
85
        
76
86
        self.__setitem__('curr_dir', value, change_userdefine=change_userdefine)
77
87
        
78
88
    ############################################################################
79
89
    def post_set_seed(self, value, change_userdefine, raiseerror):
80
90
        """ special handling for set seed """
81
91
        
82
92
        random.seed(value)
83
93
84
94
    ############################################################################        
85
95
    def post_set_run_card(self, value, change_userdefine, raiseerror):
86
96
        """ special handling for set run_card """
87
97
        
88
98
        if value == 'default':
89
99
            self.run_card = None
90
100
        elif os.path.isfile(value):
91
101
            self.run_card = banner.RunCard(value)
92
102
            
93
103
        args = value.split()
94
104
        if  len(args) >2:
95
105
            if not self.options['run_card']:
96
106
                self.run_card =  banner.RunCardLO()
97
107
                self.run_card.remove_all_cut()
98
108
            self.run_card[args[0]] = ' '.join(args[1:])
99
109
            
100
110
        
101
111
    ############################################################################
102
112
    def post_fixed_order(self, value, change_userdefine, raiseerror):
103
113
        """ special handling for set fixed_order """
104
114
        
105
115
        if value:
106
116
            logger.warning('Fix order madspin fails to have the correct scale information. This can bias the results!')
107
117
            logger.warning('Not all functionalities of MadSpin handle this mode correctly (only onshell mode so far).')
108
118
109
56
119
110
57
class MadSpinInterface(extended_cmd.Cmd):
120
class MadSpinInterface(extended_cmd.Cmd):
111
58
    """Basic interface for madspin"""
121
    """Basic interface for madspin"""
112
@@ -77,22 +140,7 @@
113
77
        self.mode = "madspin" # can be flat/bridge change the way the decay is done.
140
        self.mode = "madspin" # can be flat/bridge change the way the decay is done.
114
78
                              # note amc@nlo does not support bridge.
141
                              # note amc@nlo does not support bridge.
115
79
        
142
        
132
80
        self.options = {'max_weight': -1, 
143
        self.options = MadSpinOptions()
117
81
                        'curr_dir': os.path.realpath(os.getcwd()),
118
82
                        'Nevents_for_max_weigth': 0,
119
83
                        'max_weight_ps_point': 400,
120
84
                        'BW_cut':-1,
121
85
                        'nb_sigma':0,
122
86
                        'ms_dir':None,
123
87
                        'max_running_process':100,
124
88
                        'onlyhelicity': False,
125
89
                        'spinmode': "madspin",
126
90
                        'use_old_dir': False, #should be use only for faster debugging
127
91
                        'run_card': None, # define cut for spinmode==none.
128
92
                        'fixed_order': False # to activate fixed order handling of counter-event
129
93
                        }
130
94
        
131
95
133
96
        
144
        
134
97
        self.events_file = None
145
        self.events_file = None
135
98
        self.decay_processes = {}
146
        self.decay_processes = {}
136
@@ -107,7 +155,21 @@
137
107
            logger.info("Extracting the banner ...")
155
            logger.info("Extracting the banner ...")
138
108
            self.do_import(event_path)
156
            self.do_import(event_path)
139
109
            
157
            
141
110
            
158
    
142
159
    def setup_for_pure_decay(self):
143
160
        """this is for spinmode=None -> simple decay
144
161
           We go here if they are no banner.
145
162
           -> this requires that a command import model appears in the card!
146
163
        """
147
164
        
148
165
        logger.info("Setup the code for pure decay mode")
149
166
        self.proc_option = []
150
167
        self.final_state_full = ''
151
168
        self.final_state_compact = ''
152
169
        self.prod_branches = ''
153
170
        self.final_state = set()
154
171
        
155
172
     
156
111
    def do_import(self, inputfile):
173
    def do_import(self, inputfile):
157
112
        """import the event file"""
174
        """import the event file"""
158
113
        
175
        
159
@@ -126,6 +188,7 @@
160
126
        if not os.path.exists(inputfile):
188
        if not os.path.exists(inputfile):
161
127
            if inputfile.endswith('.gz'):
189
            if inputfile.endswith('.gz'):
162
128
                if not os.path.exists(inputfile[:-3]):
190
                if not os.path.exists(inputfile[:-3]):
163
191
                    misc.sprint(os.getcwd(), os.listdir('.'), inputfile, os.path.exists(inputfile), os.path.exists(inputfile[:-3]))
164
129
                    raise self.InvalidCmd('No such file or directory : %s' % inputfile)
192
                    raise self.InvalidCmd('No such file or directory : %s' % inputfile)
165
130
                else: 
193
                else: 
166
131
                    inputfile = inputfile[:-3]
194
                    inputfile = inputfile[:-3]
167
@@ -134,14 +197,22 @@
168
134
            else: 
197
            else: 
169
135
                raise self.InvalidCmd('No such file or directory : %s' % inputfile)
198
                raise self.InvalidCmd('No such file or directory : %s' % inputfile)
170
136
199
171
200
        self.inputfile = inputfile
172
201
        if self.options['spinmode'] == 'none' and \
173
202
           (self.options['input_format'] not in ['lhe','auto'] or 
174
203
             (self.options['input_format'] == 'auto' and '.lhe'  not in inputfile[-7:])):  
175
204
            self.banner = banner.Banner()
176
205
            self.setup_for_pure_decay()
177
206
            return  
178
207
        
179
137
        if inputfile.endswith('.gz'):
208
        if inputfile.endswith('.gz'):
180
138
            misc.gunzip(inputfile)
209
            misc.gunzip(inputfile)
181
139
            inputfile = inputfile[:-3]
210
            inputfile = inputfile[:-3]
182
140
183
141
        # Read the banner of the inputfile
211
        # Read the banner of the inputfile
184
142
        self.events_file = open(os.path.realpath(inputfile))
212
        self.events_file = open(os.path.realpath(inputfile))
185
143
        self.banner = banner.Banner(self.events_file)
213
        self.banner = banner.Banner(self.events_file)
186
144
        
214
        
187
215
188
145
        # Check the validity of the banner:
216
        # Check the validity of the banner:
189
146
        if 'slha' not in self.banner:
217
        if 'slha' not in self.banner:
190
147
            self.events_file = None
218
            self.events_file = None
191
@@ -163,7 +234,6 @@
192
163
                self.options['nb_sigma'] = N_sigma
234
                self.options['nb_sigma'] = N_sigma
193
164
            if self.options['BW_cut'] == -1:
235
            if self.options['BW_cut'] == -1:
194
165
                self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))
236
                self.options['BW_cut'] = float(self.banner.get_detail('run_card', 'bwcutoff'))
195
166
196
167
        else:
237
        else:
197
168
            if not self.options['Nevents_for_max_weigth']:
238
            if not self.options['Nevents_for_max_weigth']:
198
169
                self.options['Nevents_for_max_weigth'] = 75
239
                self.options['Nevents_for_max_weigth'] = 75
199
@@ -273,7 +343,7 @@
200
273
        
343
        
201
274
344
202
275
        #Check the param_card
345
        #Check the param_card
204
276
        if not bypass_check:
346
        if not (bypass_check or self.options['input_format'] in ['hepmc', 'lhe_no_banner']):
205
277
            if not hasattr(self.banner, 'param_card'):
347
            if not hasattr(self.banner, 'param_card'):
206
278
                self.banner.charge_card('slha')
348
                self.banner.charge_card('slha')
207
279
            param_card = check_param_card.ParamCard(card)
349
            param_card = check_param_card.ParamCard(card)
208
@@ -395,42 +465,10 @@
209
395
        
465
        
210
396
        args = self.split_arg(line)
466
        args = self.split_arg(line)
211
397
        self.check_set(args)
467
        self.check_set(args)
248
398
        
468
249
399
        if args[0] not in ['ms_dir', 'run_card']:
469
        self.options[args[0]] = ' '.join(args[1:])
250
400
            args[1] = args[1].lower()
470
        
251
401
        
471
216
402
        if args[0] in  ['max_weight', 'BW_effect','ms_dir', 'spinmode']:
217
403
            self.options[args[0]] = args[1]
218
404
            if args[0] == 'ms_dir':
219
405
                self.options['curr_dir'] = self.options['ms_dir']
220
406
        elif args[0] == 'seed':
221
407
            random.seed(int(args[1]))
222
408
            self.seed = int(args[1])
223
409
        elif args[0] == 'BW_cut':
224
410
            self.options[args[0]] = float(args[1])
225
411
        elif args[0] in ['onlyhelicity', 'use_old_dir']:
226
412
            self.options[args[0]] = banner.ConfigFile.format_variable(args[1], bool, args[0])
227
413
        elif args[0] in ['run_card']:
228
414
            if args[1] == 'default':
229
415
                self.options['run_card'] = None
230
416
            elif os.path.isfile(args[1]):
231
417
                self.options['run_card'] = banner.RunCard(args[1])
232
418
            elif  len(args) >2:
233
419
                if not self.options['run_card']:
234
420
                    self.options['run_card'] =  banner.RunCardLO()
235
421
                    self.options['run_card'].remove_all_cut()
236
422
                self.options['run_card'][args[1]] = args[2]
237
423
        elif args[0] in 'fixed_order':
238
424
            if args[1].lower() in ['t', 'true']:
239
425
                self.options['fixed_order'] = True
240
426
                logger.warning('Fix order madspin fails to have the correct scale information. This can bias the results!')
241
427
                logger.warning('Not all functionality of MadSpin handle this mode correctly (only onshell mode so far).')
242
428
            else:
243
429
                self.options['fixed_order'] = False
244
430
            logger.info('fixed_order options set to %s', self.options['fixed_order'])
245
431
        else:
246
432
            self.options[args[0]] = int(args[1])
247
433
    
252
434
    def complete_set(self,  text, line, begidx, endidx):
472
    def complete_set(self,  text, line, begidx, endidx):
253
435
        
473
        
254
436
474
255
@@ -559,17 +597,17 @@
256
559
597
257
560
        model_line = self.banner.get('proc_card', 'full_model_line')
598
        model_line = self.banner.get('proc_card', 'full_model_line')
258
561
599
264
562
        if not self.seed:
600
        if not self.options['seed']:
265
563
            self.seed = random.randint(0, int(30081*30081))
601
            self.options['seed'] = random.randint(0, int(30081*30081))
266
564
            self.do_set('seed %s' % self.seed)
602
            #self.do_set('seed %s' % self.seed)
267
565
            logger.info('Will use seed %s' % self.seed)
603
            logger.info('Will use seed %s' % self.options['seed'])
268
566
            self.history.insert(0, 'set seed %s' % self.seed)
604
            self.history.insert(0, 'set seed %s' % self.options['seed'])
269
567
605
272
568
        if self.seed > 30081*30081: # can't use too big random number
606
        if self.options['seed'] > 30081*30081: # can't use too big random number
273
569
            msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'
607
            msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
274
570
            raise Exception, msg
608
            raise Exception, msg
275
571
609
277
572
        self.options['seed'] = self.seed
610
        #self.options['seed'] = self.seed
278
573
        text = '%s\n' % '\n'.join([ line for line in self.history if line])
611
        text = '%s\n' % '\n'.join([ line for line in self.history if line])
279
574
        self.banner.add_text('madspin' , text)
612
        self.banner.add_text('madspin' , text)
280
575
        
613
        
281
@@ -678,9 +716,9 @@
282
678
        
716
        
283
679
        # NOW we have all the information available for RUNNING
717
        # NOW we have all the information available for RUNNING
284
680
        
718
        
286
681
        if self.seed:
719
        if self.options['seed']:
287
682
            #seed is specified need to use that one:
720
            #seed is specified need to use that one:
289
683
            open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.seed)
721
            open(pjoin(self.options['ms_dir'],'seeds.dat'),'w').write('%s\n'%self.options['seed'])
290
684
            #remove all ranmar_state
722
            #remove all ranmar_state
291
685
            for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'), 
723
            for name in misc.glob(pjoin('*', 'SubProcesses','*','ranmar_state.dat'), 
292
686
                                                        self.options['ms_dir']):
724
                                                        self.options['ms_dir']):
293
@@ -741,30 +779,61 @@
294
741
            for name in misc.glob("decay_*_*", self.path_me):
779
            for name in misc.glob("decay_*_*", self.path_me):
295
742
                shutil.rmtree(name)
780
                shutil.rmtree(name)
296
743
781
300
744
        self.events_file.close()
782
        if self.events_file:
301
745
        orig_lhe = lhe_parser.EventFile(self.events_file.name)
783
            self.events_file.close()
302
746
        
784
            filename = self.events_file.name
303
785
        else:
304
786
            filename = self.inputfile
305
787
306
788
        if self.options['input_format'] == 'auto':
307
789
            if '.lhe' in filename :
308
790
                self.options['input_format']  = 'lhe'
309
791
            elif '.hepmc' in filename:
310
792
                self.options['input_format']  = 'hepmc'
311
793
            else:
312
794
                raise Exception, "fail to recognized input format automatically"
313
795
                
314
796
        if self.options['input_format'] in ['lhe', 'lhe_no_banner']:
315
797
            orig_lhe = lhe_parser.EventFile(filename)
316
798
            if self.options['input_format'] == 'lhe_no_banner':
317
799
                orig_lhe.allow_empty_event = True
318
800
                
319
801
        elif self.options['input_format'] in ['hepmc']:
320
802
            import madgraph.various.hepmc_parser as hepmc_parser
321
803
            orig_lhe = hepmc_parser.HEPMC_EventFile(filename)
322
804
            logger.info("Parsing input event to know how many decay to generate. This can takes few minuts.")
323
805
        else:
324
806
            raise Exception
325
807
            
326
747
        to_decay = collections.defaultdict(int)
808
        to_decay = collections.defaultdict(int)
327
748
        nb_event = 0
809
        nb_event = 0
328
810
 
329
749
        for event in orig_lhe:
811
        for event in orig_lhe:
330
750
            nb_event +=1
812
            nb_event +=1
331
751
            for particle in event:
813
            for particle in event:
332
752
                if particle.status == 1 and particle.pdg in asked_to_decay:
814
                if particle.status == 1 and particle.pdg in asked_to_decay:
333
753
                    # final state and tag as to decay
815
                    # final state and tag as to decay
334
754
                    to_decay[particle.pdg] += 1
816
                    to_decay[particle.pdg] += 1
335
817
            if self.options['input_format'] == 'hepmc' and nb_event == 250:
336
818
                currpos = orig_lhe.tell()
337
819
                filesize = orig_lhe.getfilesize()
338
820
                for key in to_decay:
339
821
                    to_decay[key] *= 1.05 * filesize/ currpos 
340
822
                    # 1.05 to avoid accidental coincidence with nevents
341
823
                break
342
755
824
343
756
        # Handle the banner of the output file
825
        # Handle the banner of the output file
349
757
        if not self.seed:
826
        if not self.options['seed']:
350
758
            self.seed = random.randint(0, int(30081*30081))
827
            self.options['seed'] = random.randint(0, int(30081*30081))
351
759
            self.do_set('seed %s' % self.seed)
828
            #self.do_set('seed %s' % self.seed)
352
760
            logger.info('Will use seed %s' % self.seed)
829
            logger.info('Will use seed %s' % self.options['seed'])
353
761
            self.history.insert(0, 'set seed %s' % self.seed)
830
            self.history.insert(0, 'set seed %s' % self.options['seed'])
354
762
831
357
763
        if self.seed > 30081*30081: # can't use too big random number
832
        if self.options['seed'] > 30081*30081: # can't use too big random number
358
764
            msg = 'Random seed too large ' + str(self.seed) + ' > 30081*30081'
833
            msg = 'Random seed too large ' + str(self.options['seed']) + ' > 30081*30081'
359
765
            raise Exception, msg
834
            raise Exception, msg
360
766
835
362
767
        self.options['seed'] = self.seed
836
        #self.options['seed'] = self.options['seed']
363
768
        
837
        
364
769
        text = '%s\n' % '\n'.join([ line for line in self.history if line])
838
        text = '%s\n' % '\n'.join([ line for line in self.history if line])
365
770
        self.banner.add_text('madspin' , text)
839
        self.banner.add_text('madspin' , text)
366
@@ -780,7 +849,7 @@
367
780
            for pdg, nb_needed in to_decay.items():
849
            for pdg, nb_needed in to_decay.items():
368
781
                #check if a splitting is needed
850
                #check if a splitting is needed
369
782
                if nb_needed == nb_event:
851
                if nb_needed == nb_event:
371
783
                    evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5)
852
                    evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5)
372
784
                elif nb_needed %  nb_event == 0:
853
                elif nb_needed %  nb_event == 0:
373
785
                    nb_mult = nb_needed // nb_event
854
                    nb_mult = nb_needed // nb_event
374
786
                    part = self.model.get_particle(pdg)
855
                    part = self.model.get_particle(pdg)
375
@@ -788,64 +857,107 @@
376
788
                    if name not in self.list_branches:
857
                    if name not in self.list_branches:
377
789
                        continue
858
                        continue
378
790
                    elif len(self.list_branches[name]) == nb_mult:
859
                    elif len(self.list_branches[name]) == nb_mult:
382
791
                        evt_decayfile[pdg] = self.generate_events(pdg, nb_event, mg5)
860
                        evt_decayfile[pdg] = self.generate_events(pdg, min(nb_event,100000), mg5)
383
792
                    else:
861
                    else:
384
793
                        evt_decayfile[pdg] = self.generate_events(pdg, nb_needed, mg5, cumul=True)
862
                        evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5, cumul=True)
385
863
                elif self.options['cross_section']:
386
864
                    #cross-section hard-coded -> allow 
387
865
                    part = self.model.get_particle(pdg)
388
866
                    name = part.get_name()
389
867
                    
390
868
                    if name not in self.list_branches:
391
869
                        continue
392
870
                    else:
393
871
                        try:
394
872
                            evt_decayfile[pdg] = self.generate_events(pdg, min(nb_needed,100000), mg5, cumul=True)
395
873
                        except common_run_interface.ZeroResult:
396
874
                            logger.warning("Branching ratio is zero for this particle. Not decaying it")
397
875
                            del to_decay[pdg]                    
398
794
                else:
876
                else:
399
795
                    part = self.model.get_particle(pdg)
877
                    part = self.model.get_particle(pdg)
400
796
                    name = part.get_name()
878
                    name = part.get_name()
401
797
                    if name not in self.list_branches or len(self.list_branches[name]) == 0:
879
                    if name not in self.list_branches or len(self.list_branches[name]) == 0:
402
798
                        continue
880
                        continue
404
799
                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
881
                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed. One workaround is to force the final cross-section manually.")
405
800
                    
882
                    
406
801
                     
883
                     
407
802
                
408
803
        
409
804
        # Compute the branching ratio.
884
        # Compute the branching ratio.
435
805
        br = 1
885
        if not self.options['cross_section']:
436
806
        for (pdg, event_files) in evt_decayfile.items():
886
            br = 1
437
807
            if not event_files:
887
            for (pdg, event_files) in evt_decayfile.items():
438
808
                continue
888
                if not event_files:
439
809
            totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
889
                    continue
440
810
            if to_decay[pdg] == nb_event:
890
                totwidth = float(self.banner.get('param', 'decay', abs(pdg)).value)
441
811
                # Exactly one particle of this type to decay by event
891
                if to_decay[pdg] == nb_event:
442
812
                pwidth = sum([event_files[k].cross for k in event_files])
892
                    # Exactly one particle of this type to decay by event
443
813
                if pwidth > 1.01 * totwidth:
893
                    pwidth = sum([event_files[k].cross for k in event_files])
419
814
                    logger.critical("Branching ratio larger than one for %s " % pdg) 
420
815
                br *= pwidth / totwidth
421
816
            elif to_decay[pdg] % nb_event == 0:
422
817
                # More than one particle of this type to decay by event
423
818
                # Need to check the number of event file to check if we have to 
424
819
                # make separate type of decay or not.
425
820
                nb_mult = to_decay[pdg] // nb_event
426
821
                if nb_mult == len(event_files):
427
822
                    for k in event_files:
428
823
                        pwidth = event_files[k].cross
429
824
                        if pwidth > 1.01 * totwidth:
430
825
                            logger.critical("Branching ratio larger than one for %s " % pdg)                       
431
826
                        br *= pwidth / totwidth
432
827
                    br *= math.factorial(nb_mult)
433
828
                else:
434
829
                    pwidth = sum(event_files[k].cross for k in event_files)
444
830
                    if pwidth > 1.01 * totwidth:
894
                    if pwidth > 1.01 * totwidth:
445
831
                        logger.critical("Branching ratio larger than one for %s " % pdg) 
895
                        logger.critical("Branching ratio larger than one for %s " % pdg) 
449
832
                    br *= (pwidth / totwidth)**nb_mult
896
                    br *= pwidth / totwidth
450
833
            else:
897
                elif to_decay[pdg] % nb_event == 0:
451
834
                raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
898
                    # More than one particle of this type to decay by event
452
899
                    # Need to check the number of event file to check if we have to 
453
900
                    # make separate type of decay or not.
454
901
                    nb_mult = to_decay[pdg] // nb_event
455
902
                    if nb_mult == len(event_files):
456
903
                        for k in event_files:
457
904
                            pwidth = event_files[k].cross
458
905
                            if pwidth > 1.01 * totwidth:
459
906
                                logger.critical("Branching ratio larger than one for %s " % pdg)                       
460
907
                            br *= pwidth / totwidth
461
908
                        br *= math.factorial(nb_mult)
462
909
                    else:
463
910
                        pwidth = sum(event_files[k].cross for k in event_files)
464
911
                        if pwidth > 1.01 * totwidth:
465
912
                            logger.critical("Branching ratio larger than one for %s " % pdg) 
466
913
                        br *= (pwidth / totwidth)**nb_mult
467
914
                else:
468
915
                    raise self.InvalidCmd("The bridge mode of MadSpin does not support event files where events do not *all* share the same set of final state particles to be decayed.")
469
916
        else:
470
917
            br = 1
471
835
        self.branching_ratio = br
918
        self.branching_ratio = br
472
836
        self.efficiency = 1
919
        self.efficiency = 1
474
837
        self.cross, self.error = self.banner.get_cross(witherror=True)
920
        try:
475
921
            self.cross, self.error = self.banner.get_cross(witherror=True)
476
922
        except:
477
923
            if self.options['input_format'] != 'lhe':
478
924
                self.cross, self.error = 0, 0
479
838
        self.cross *= br
925
        self.cross *= br
480
839
        self.error *= br
926
        self.error *= br
481
840
        
927
        
482
841
        
483
842
        # modify the cross-section in the init block of the banner
928
        # modify the cross-section in the init block of the banner
487
843
        self.banner.scale_init_cross(self.branching_ratio)
929
        if not self.options['cross_section']:
488
844
                    
930
            self.banner.scale_init_cross(self.branching_ratio)
489
845
        
931
        else:
490
932
            
491
933
            if self.options['input_format'] in ['lhe_no_banner','hepmc'] and 'init' not in self.banner:
492
934
                self.cross = sum(self.options['cross_section'].values())
493
935
                self.error = 0
494
936
                self.branching_ratio = 1
495
937
            else:  
496
938
                self.banner.modify_init_cross(self.options['cross_section'])
497
939
                new_cross, new_error =   self.banner.get_cross(witherror=True)
498
940
                self.branching_ratio = new_cross / self.cross
499
941
                self.cross = new_cross   
500
942
                self.error = new_error
501
943
502
846
        # 3. Merge the various file together.
944
        # 3. Merge the various file together.
505
847
        output_lhe = lhe_parser.EventFile(orig_lhe.name.replace('.lhe', '_decayed.lhe.gz'), 'w')
945
        if self.options['input_format'] == 'hepmc':
506
848
        self.banner.write(output_lhe, close_tag=False)
946
            name = orig_lhe.name.replace('.hepmc', '_decayed.lhe')
507
947
            if not name.endswith('.gz'):
508
948
                name = '%s.gz' % name
509
949
            
510
950
            output_lhe = lhe_parser.EventFile(name, 'w')
511
951
        else:
512
952
            name = orig_lhe.name.replace('.lhe', '_decayed.lhe')
513
953
            if not name.endswith('.gz'):
514
954
                name = '%s.gz' % name
515
955
            output_lhe = lhe_parser.EventFile(name, 'w')
516
956
        try:
517
957
            self.banner.write(output_lhe, close_tag=False)
518
958
        except Exception:
519
959
            if self.options['input_format'] == 'lhe':
520
960
                raise
521
849
        
961
        
522
850
        # initialise object which store not use event due to wrong helicity
962
        # initialise object which store not use event due to wrong helicity
523
851
        bufferedEvents_decay = {}
963
        bufferedEvents_decay = {}
524
@@ -856,8 +968,9 @@
525
856
        start = time.time()
968
        start = time.time()
526
857
        counter = 0
969
        counter = 0
527
858
        orig_lhe.seek(0)
970
        orig_lhe.seek(0)
528
971
529
859
        for event in orig_lhe:
972
        for event in orig_lhe:
531
860
            if counter and counter % 1000 == 0 and float(str(counter)[1:]) ==0:
973
            if counter and counter % 100 == 0 and float(str(counter)[1:]) ==0:
532
861
                print "decaying event number %s [%s s]" % (counter, time.time()-start)
974
                print "decaying event number %s [%s s]" % (counter, time.time()-start)
533
862
            counter +=1
975
            counter +=1
534
863
            
976
            
535
@@ -866,10 +979,17 @@
536
866
            particles = [p for p in event if int(p.status) == 1.0]
979
            particles = [p for p in event if int(p.status) == 1.0]
537
867
            random.shuffle(particles)
980
            random.shuffle(particles)
538
868
            ids = [particle.pid for particle in particles]
981
            ids = [particle.pid for particle in particles]
539
982
            br = 1 #br for that particular events (for special/weighted case)
540
983
            hepmc_output = lhe_parser.Event() #for hepmc case: collect the decay particle
541
869
            for i,particle in enumerate(particles):
984
            for i,particle in enumerate(particles):
542
985
                #misc.sprint(i, particle.pdg, particle.pid)
543
986
                #misc.sprint(self.final_state, evt_decayfile)
544
870
                # check if we need to decay the particle 
987
                # check if we need to decay the particle 
547
871
                if particle.pdg not in self.final_state or particle.pdg not in evt_decayfile:
988
                if self.final_state and particle.pdg not in self.final_state:
548
872
                    continue # nothing to do for this particle
989
                    continue # nothing to do for this particle
549
990
                if particle.pdg not in evt_decayfile:
550
991
                    continue # nothing to do for this particle
551
992
                
552
873
                # check how the decay need to be done
993
                # check how the decay need to be done
553
874
                nb_decay = len(evt_decayfile[particle.pdg])
994
                nb_decay = len(evt_decayfile[particle.pdg])
554
875
                if nb_decay == 0:
995
                if nb_decay == 0:
555
@@ -888,14 +1008,28 @@
556
888
                    cumul = 0
1008
                    cumul = 0
557
889
                    for j,events in evt_decayfile[particle.pdg].items():
1009
                    for j,events in evt_decayfile[particle.pdg].items():
558
890
                        cumul += events.cross
1010
                        cumul += events.cross
560
891
                        if r < cumul:
1011
                        if r <= cumul:
561
892
                            decay_file = events
1012
                            decay_file = events
562
893
                            decay_file_nb = j
1013
                            decay_file_nb = j
563
894
                        else:
564
895
                            break
1014
                            break
566
896
                        
1015
                    else:
567
1016
                        # security for numerical accuracy issue... (unlikely but better safe)
568
1017
                        if (cumul-tot)/tot < 1e-5:
569
1018
                            decay_file = events
570
1019
                            decay_file_nb = j
571
1020
                        else:
572
1021
                            misc.sprint(j,cumul, events.cross, tot, (tot-cumul)/tot)
573
1022
                            raise Exception
574
1023
                
575
1024
                if self.options['new_wgt'] == 'BR':
576
1025
                    tot_width = float(self.banner.get('param', 'decay', abs(pdg)).value)
577
1026
                    if tot_width:
578
1027
                        br = decay_file.cross / tot_width
579
897
                # ok start the procedure
1028
                # ok start the procedure
581
898
                helicity = particle.helicity
1029
                if hasattr(particle,'helicity'):
582
1030
                    helicity = particle.helicity
583
1031
                else:
584
1032
                    helicity = 9
585
899
                bufferedEvents = bufferedEvents_decay[particle.pdg][decay_file_nb]
1033
                bufferedEvents = bufferedEvents_decay[particle.pdg][decay_file_nb]
586
900
                
1034
                
587
901
                # now that we have the file to read. find the associate event
1035
                # now that we have the file to read. find the associate event
588
@@ -910,8 +1044,8 @@
589
910
                        except StopIteration:
1044
                        except StopIteration:
590
911
                            # check how far we are
1045
                            # check how far we are
591
912
                            ratio = counter / nb_event 
1046
                            ratio = counter / nb_event 
594
913
                            needed = 1.05 * to_decay[particle.pdg]/ratio - counter
1047
                            needed = 1.05 * to_decay[particle.pdg] - counter
595
914
                            needed = min(1000, max(needed, 1000))
1048
                            needed = min(100000, max(needed, 6000))
596
915
                            with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
1049
                            with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
597
916
                                new_file = self.generate_events(particle.pdg, needed, mg5, [decay_file_nb])
1050
                                new_file = self.generate_events(particle.pdg, needed, mg5, [decay_file_nb])
598
917
                            evt_decayfile[particle.pdg].update(new_file)
1051
                            evt_decayfile[particle.pdg].update(new_file)
599
@@ -927,14 +1061,36 @@
600
927
                            # only add to the buffering if the buffer is not too large
1061
                            # only add to the buffering if the buffer is not too large
601
928
                            bufferedEvents[helicity].append(decay)
1062
                            bufferedEvents[helicity].append(decay)
602
929
                # now that we have the event make the merge
1063
                # now that we have the event make the merge
604
930
                particle.add_decay(decay)
1064
                if self.options['input_format'] != 'hepmc':
605
1065
                    particle.add_decay(decay)
606
1066
                else:
607
1067
                    if len(hepmc_output) == 0:
608
1068
                        hepmc_output.append(lhe_parser.Particle(event=hepmc_output))
609
1069
                        hepmc_output[0].color2 = 0
610
1070
                    decayed_particle = lhe_parser.Particle(particle, hepmc_output)
611
1071
                    hepmc_output.append(decayed_particle)
612
1072
                    decayed_particle.add_decay(decay)
613
931
            # change the weight associate to the event
1073
            # change the weight associate to the event
620
932
            event.wgt *= self.branching_ratio
1074
            if self.options['new_wgt'] == 'cross-section':
621
933
            wgts = event.parse_reweight()
1075
                event.wgt *= self.branching_ratio
622
934
            for key in wgts:
1076
                br = self.branching_ratio
623
935
                wgts[key] *= self.branching_ratio
1077
            else:
624
936
            # all particle have been decay if needed
1078
                event.wgt *= br
625
937
            output_lhe.write(str(event))
1079
                
626
1080
            if self.options['input_format'] != 'hepmc':
627
1081
                wgts = event.parse_reweight()
628
1082
                for key in wgts:
629
1083
                    wgts[key] *= br
630
1084
                # all particle have been decay if needed
631
1085
                output_lhe.write(str(event))
632
1086
            else:
633
1087
                hepmc_output.wgt = event.wgt
634
1088
                hepmc_output.nexternal = len(hepmc_output) # the append does not update nexternal
635
1089
                hepmc_output.assign_mother()
636
1090
                output_lhe.write(str(hepmc_output))
637
1091
        else:
638
1092
            if counter==0:
639
1093
                raise Exception
640
938
        output_lhe.write('</LesHouchesEvents>\n')        
1094
        output_lhe.write('</LesHouchesEvents>\n')        
641
939
                    
1095
                    
642
940
    
1096
    
643
@@ -971,7 +1127,12 @@
644
971
           cumul allow to merge all the definition in one run (add process)
1127
           cumul allow to merge all the definition in one run (add process)
645
972
                 to generate events according to cross-section
1128
                 to generate events according to cross-section
646
973
        """
1129
        """
648
974
        
1130
        if not hasattr(self, 'me_int'):
649
1131
            self.me_int = {}
650
1132
            
651
1133
        
652
1134
        
653
1135
        nb_event = int(nb_event) # in case of hepmc request the nb_event is not an integer
654
975
        if cumul:
1136
        if cumul:
655
976
            width = 0.
1137
            width = 0.
656
977
        else:   
1138
        else:   
657
@@ -1003,31 +1164,36 @@
658
1003
                options = dict(mg5.options)
1164
                options = dict(mg5.options)
659
1004
                if self.options['ms_dir']:
1165
                if self.options['ms_dir']:
660
1005
                    # we are in gridpack mode -> create it
1166
                    # we are in gridpack mode -> create it
673
1006
                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
1167
                    if decay_dir in self.me_int:
674
1007
                                            decay_dir), options=options)
1168
                        me5_cmd = self.me_int[decay_dir]
675
1008
                    me5_cmd.options["automatic_html_opening"] = False
1169
                    else:
676
1009
                    me5_cmd.options["madanalysis5_path"] = None
1170
                        me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
677
1010
                    me5_cmd.options["madanalysis_path"] = None
1171
                                                decay_dir), options=options)
678
1011
                    me5_cmd.allow_notification_center = False
1172
                        me5_cmd.options["automatic_html_opening"] = False
679
1012
                    try:
1173
                        me5_cmd.options["madanalysis5_path"] = None
680
1013
                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
1174
                        me5_cmd.options["madanalysis_path"] = None
681
1014
                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
1175
                        me5_cmd.allow_notification_center = False
682
1015
                    except Exception,error:
1176
                        try:
683
1016
                        logger.debug(error)
1177
                            os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
684
1017
                        pass            
1178
                            os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
685
1179
                        except Exception,error:
686
1180
                            logger.debug(error)
687
1181
                            pass 
688
1182
                        self.me_int[decay_dir] = me5_cmd
689
1183
690
1018
                    if self.options["run_card"]:
1184
                    if self.options["run_card"]:
692
1019
                        run_card = self.options["run_card"]
1185
                        run_card = self.run_card
693
1020
                    else:
1186
                    else:
694
1021
                        run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
1187
                        run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))                        
695
1022
                    
1188
                    
697
1023
                    run_card["iseed"] = self.seed
1189
                    run_card["iseed"] = self.options['seed']
698
1024
                    run_card['gridpack'] = True
1190
                    run_card['gridpack'] = True
699
1025
                    run_card['systematics_program'] = 'False'
1191
                    run_card['systematics_program'] = 'False'
700
1026
                    run_card['use_syst'] = False
1192
                    run_card['use_syst'] = False
701
1027
                    run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
1193
                    run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
702
1028
                    param_card = self.banner['slha']
1194
                    param_card = self.banner['slha']
703
1029
                    open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
1195
                    open(pjoin(decay_dir, "Cards", "param_card.dat"),"w").write(param_card)
705
1030
                    self.seed += 1
1196
                    self.options['seed'] += 1
706
1031
                    # actually creation
1197
                    # actually creation
707
1032
                    me5_cmd.exec_cmd("generate_events run_01 -f")
1198
                    me5_cmd.exec_cmd("generate_events run_01 -f")
708
1033
                    if output_width:
1199
                    if output_width:
709
@@ -1041,26 +1207,35 @@
710
1041
                    misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
1207
                    misc.call(['tar', '-xzpvf', 'run_01_gridpack.tar.gz'], cwd=decay_dir)
711
1042
            
1208
            
712
1043
            # Now generate the events
1209
            # Now generate the events
713
1044
714
1045
            if not self.options['ms_dir']:
1210
            if not self.options['ms_dir']:
728
1046
                me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
1211
                if decay_dir in self.me_int:
729
1047
                                                decay_dir), options=mg5.options)
1212
                        me5_cmd = self.me_int[decay_dir]
730
1048
                me5_cmd.options["automatic_html_opening"] = False
1213
                else:
731
1049
                me5_cmd.options["automatic_html_opening"] = False
1214
                    me5_cmd = madevent_interface.MadEventCmdShell(me_dir=os.path.realpath(\
732
1050
                me5_cmd.options["madanalysis5_path"] = None
1215
                                                    decay_dir), options=mg5.options)
733
1051
                me5_cmd.options["madanalysis_path"] = None
1216
                    me5_cmd.options["automatic_html_opening"] = False
734
1052
                me5_cmd.allow_notification_center = False
1217
                    me5_cmd.options["automatic_html_opening"] = False
735
1053
                try:
1218
                    me5_cmd.options["madanalysis5_path"] = None
736
1054
                    os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
1219
                    me5_cmd.options["madanalysis_path"] = None
737
1055
                    os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
1220
                    me5_cmd.allow_notification_center = False
738
1056
                except Exception,error:
1221
                    try:
739
1057
                    logger.debug(error)
1222
                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card_default.dat'))
740
1058
                    pass                    
1223
                        os.remove(pjoin(decay_dir, 'Cards', 'madanalysis5_parton_card.dat'))
741
1224
                    except Exception,error:
742
1225
                        logger.debug(error)
743
1226
                        pass                 
744
1227
                    self.me_int[decay_dir] = me5_cmd
745
1059
                if self.options["run_card"]:
1228
                if self.options["run_card"]:
747
1060
                    run_card = self.options["run_card"]
1229
                    run_card = self.run_card
748
1061
                else:
1230
                else:
749
1062
                    run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
1231
                    run_card = banner.RunCard(pjoin(decay_dir, "Cards", "run_card.dat"))
750
1063
                run_card["nevents"] = int(1.2*nb_event)
1232
                run_card["nevents"] = int(1.2*nb_event)
751
1233
                # Handle the banner of the output file
752
1234
                if not self.seed:
753
1235
                    self.seed = random.randint(0, int(30081*30081))
754
1236
                    self.do_set('seed %s' % self.seed)
755
1237
                    logger.info('Will use seed %s' % self.seed)
756
1238
                    self.history.insert(0, 'set seed %s' % self.seed)
757
1064
                run_card["iseed"] = self.seed
1239
                run_card["iseed"] = self.seed
758
1065
                run_card["systematics_program"] = 'None'
1240
                run_card["systematics_program"] = 'None'
759
1066
                run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
1241
                run_card.write(pjoin(decay_dir, "Cards", "run_card.dat"))
760
@@ -1140,8 +1315,7 @@
761
1140
                if particle.status == 1 and particle.pdg in asked_to_decay:
1315
                if particle.status == 1 and particle.pdg in asked_to_decay:
762
1141
                    # final state and tag as to decay
1316
                    # final state and tag as to decay
763
1142
                    to_decay[particle.pdg] += 1
1317
                    to_decay[particle.pdg] += 1
766
1143
        #misc.sprint(to_decay)
1318
765
1144
        #misc.sprint("import the mode -> temporary with logging")   
767
1145
        with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
1319
        with misc.MuteLogger(["madgraph", "madevent", "ALOHA", "cmdprint"], [50,50,50,50]):
768
1146
            mg5 = self.mg5cmd
1320
            mg5 = self.mg5cmd
769
1147
            if not self.model:
1321
            if not self.model:
770
@@ -1436,10 +1610,6 @@
771
1436
        full_event = lhe_parser.Event(str(production))
1610
        full_event = lhe_parser.Event(str(production))
772
1437
        full_event = full_event.add_decays(decays)
1611
        full_event = full_event.add_decays(decays)
773
1438
        full_me = self.calculate_matrix_element(full_event)
1612
        full_me = self.calculate_matrix_element(full_event)
774
1439
        #misc.sprint(full_event)
775
1440
        #misc.sprint([p.pdg for p in production])
776
1441
        #misc.sprint([p.pdg for p in full_event])
777
1442
        #misc.sprint(full_me, production_me, decay_me)
778
1443
        return full_event, full_me/(production_me*decay_me)
1613
        return full_event, full_me/(production_me*decay_me)
779
1444
        
1614
        
780
1445
        
1615
        
781
@@ -1457,7 +1627,6 @@
782
1457
            init = (-init[0],)
1627
            init = (-init[0],)
783
1458
            final = tuple(-i for i in final)
1628
            final = tuple(-i for i in final)
784
1459
            tag = (init, final)
1629
            tag = (init, final)
785
1460
            misc.sprint([k for k in  self.all_me.keys() if len(k[0])==1])
786
1461
            orig_order = self.all_me[tag]['order']
1630
            orig_order = self.all_me[tag]['order']
787
1462
            
1631
            
788
1463
        pdir = self.all_me[tag]['pdir']
1632
        pdir = self.all_me[tag]['pdir']
789
@@ -1520,7 +1689,6 @@
790
1520
                
1689
                
791
1521
            rwgt_interface.ReweightInterface.get_LO_definition_from_NLO()        
1690
            rwgt_interface.ReweightInterface.get_LO_definition_from_NLO()        
792
1522
        
1691
        
793
1523
        
794
1524
1692
795
1525
1693
796
1526
if __name__ == '__main__':
1694
if __name__ == '__main__':
797
1527
1695
798
=== modified file 'MadSpin/madspin'
799
--- MadSpin/madspin	2013-11-12 23:52:05 +0000
800
+++ MadSpin/madspin	2018-12-20 08:18:55 +0000
801
@@ -49,6 +49,7 @@
802
49
if len(args) == 0:
49
if len(args) == 0:
803
50
    args = ''
50
    args = ''
804
51
51
805
52
806
52
import subprocess
53
import subprocess
807
53
54
808
54
# Check if optimize mode is (and should be) activated
55
# Check if optimize mode is (and should be) activated
809
@@ -114,23 +115,6 @@
810
114
    pass
115
    pass
811
115
import MadSpin.interface_madspin as interface
116
import MadSpin.interface_madspin as interface
812
116
117
813
117
try:
814
118
    cmd_line = interface.MadSpinInterface()
815
119
    cmd_line.use_rawinput = True
816
120
    cmd_line.cmdloop()
817
121
except:
818
122
    pass
819
123
try:
820
124
    cmd_line.exec_cmd('quit all', printcmd=False)
821
125
    readline.set_history_length(100)
822
126
    if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
823
127
        os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
824
128
    readline.write_history_file(history_file)
825
129
except Exception, error:
826
130
    pass
827
131
828
132
sys.exit(0)
829
133
830
134
# Call the cmd interface main loop
118
# Call the cmd interface main loop
831
135
try:
119
try:
832
136
    if options.file or args:
120
    if options.file or args:
833
@@ -143,12 +127,14 @@
834
143
            cmd_line = interface.MasterCmdWeb()
127
            cmd_line = interface.MasterCmdWeb()
835
144
            cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log')
128
            cmd_line.debug_output = os.path.join(os.path.dirname(input_file),'generation.log')
836
145
            cmd_line.use_rawinput = False
129
            cmd_line.use_rawinput = False
838
146
            cmd_line.run_cmd('import ' + input_file)
130
            cmd_line.import_command_file(input_file)
839
147
            cmd_line.run_cmd('quit')
131
            cmd_line.run_cmd('quit')
840
148
        else:
132
        else:
841
133
            print "using input+file", input_file
842
149
            cmd_line = interface.MadSpinInterface()
134
            cmd_line = interface.MadSpinInterface()
843
150
            cmd_line.use_rawinput = False
135
            cmd_line.use_rawinput = False
845
151
            cmd_line.run_cmd('import ' + input_file)
136
            cmd_line.haspiping = False
846
137
            cmd_line.import_command_file(input_file)
847
152
            cmd_line.run_cmd('quit')
138
            cmd_line.run_cmd('quit')
848
153
    else:
139
    else:
849
154
        # Interactive mode
140
        # Interactive mode
850
@@ -160,8 +146,20 @@
851
160
            cmd_line = interface.MasterCmdWeb()
146
            cmd_line = interface.MasterCmdWeb()
852
161
            cmd_line.cmdloop()
147
            cmd_line.cmdloop()
853
162
        else:
148
        else:
856
163
            cmd_line = interface.MasterCmd(mgme_dir = options.mgme_dir)
149
            try:
857
164
            cmd_line.cmdloop()
150
                cmd_line = interface.MadSpinInterface()
858
151
                cmd_line.use_rawinput = True
859
152
                cmd_line.cmdloop()
860
153
            except:
861
154
                pass
862
155
        try:
863
156
            cmd_line.exec_cmd('quit all', printcmd=False)
864
157
            readline.set_history_length(100)
865
158
            if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
866
159
                os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
867
160
            readline.write_history_file(history_file)
868
161
        except Exception, error:
869
162
            pass
870
165
except KeyboardInterrupt:
163
except KeyboardInterrupt:
871
166
    print 'writting history and quit on KeyboardInterrupt' 
164
    print 'writting history and quit on KeyboardInterrupt' 
872
167
    pass
165
    pass
873
168
166
874
=== modified file 'MadSpin/src/driver.f'
875
--- MadSpin/src/driver.f	2017-08-03 18:49:25 +0000
876
+++ MadSpin/src/driver.f	2018-12-20 08:18:55 +0000
877
@@ -340,10 +340,10 @@
878
340
           if (jac.lt.0d0) then
340
           if (jac.lt.0d0) then
879
341
             counter2=counter2+1 
341
             counter2=counter2+1 
880
342
             counter3=counter3+1 
342
             counter3=counter3+1 
885
343
c             if (counter3.gt.500) then
343
             if (counter3.gt.500*8*100) then
886
344
c               write(*,*) "500_pts_failed_stop_executation"
344
               write(*,*) "500_pts_failed_stop_executation"
887
345
c               stop
345
               stop
888
346
c             endif
346
             endif
889
347
             if (counter2.ge.8) then ! use another topology to generate PS points
347
             if (counter2.ge.8) then ! use another topology to generate PS points
890
348
               do k=1,n_max_cg
348
               do k=1,n_max_cg
891
349
                 amp2(k)=0d0
349
                 amp2(k)=0d0
892
@@ -960,7 +960,7 @@
893
960
      include 'nexternal_prod.inc'
960
      include 'nexternal_prod.inc'
894
961
      !include 'run.inc'
961
      !include 'run.inc'
895
962
c arguments
962
c arguments
897
963
      double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal)
963
      double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal_prod)
898
964
      integer itree(2,-nexternal:-1)
964
      integer itree(2,-nexternal:-1)
899
965
      double precision qmass(-nexternal:0),qwidth(-nexternal:0)
965
      double precision qmass(-nexternal:0),qwidth(-nexternal:0)
900
966
c common
966
c common
901
967
967
902
=== modified file 'Template/LO/Cards/run_card.dat'
903
--- Template/LO/Cards/run_card.dat	2018-04-25 21:20:23 +0000
904
+++ Template/LO/Cards/run_card.dat	2018-12-20 08:18:55 +0000
905
@@ -12,7 +12,7 @@
906
12
#   mind the format:   value    = variable     ! comment             *
12
#   mind the format:   value    = variable     ! comment             *
907
13
#                                                                    *
13
#                                                                    *
908
14
#   To display more options, you can type the command:               *
14
#   To display more options, you can type the command:               *
910
15
*      update full_run_card                                          *
15
#      update full_run_card                                          *
911
16
#*********************************************************************
16
#*********************************************************************
912
17
#
17
#
913
18
#*******************                                                 
18
#*******************                                                 
914
19
19
915
=== modified file 'Template/LO/SubProcesses/genps.f'
916
--- Template/LO/SubProcesses/genps.f	2018-04-20 08:38:20 +0000
917
+++ Template/LO/SubProcesses/genps.f	2018-12-20 08:18:55 +0000
918
@@ -41,7 +41,7 @@
919
41
      integer mincfig,maxcfig                  !Range of configurations
41
      integer mincfig,maxcfig                  !Range of configurations
920
42
      integer invar
42
      integer invar
921
43
      double precision wgt                     !(input and output)
43
      double precision wgt                     !(input and output)
923
44
      double precision x(maxdim),p(maxdim)     !x,p (output) [p(0:3,nexternal)]
44
      double precision x(*),p(*)     !x,p (output) [p(0:3,nexternal)]
924
45
c
45
c
925
46
c     Local
46
c     Local
926
47
c
47
c
927
@@ -104,7 +104,7 @@
928
104
c
104
c
929
105
      integer iconfig,mincfig,maxcfig,invar
105
      integer iconfig,mincfig,maxcfig,invar
930
106
      double precision p1(0:3,nexternal+1)
106
      double precision p1(0:3,nexternal+1)
932
107
      double precision x(maxinvar)
107
      double precision x(*)
933
108
      double precision wgt
108
      double precision wgt
934
109
c
109
c
935
110
c     Local
110
c     Local
936
@@ -191,8 +191,6 @@
937
191
c
191
c
938
192
      include 'configs.inc'
192
      include 'configs.inc'
939
193
      data firsttime/.true./
193
      data firsttime/.true./
940
194
      integer isym(0:100)
941
195
c      data isym /2,1,5,27,42,47,0,0,0,0,0/
942
196
      data jfig/1/
194
      data jfig/1/
943
197
c-----
195
c-----
944
198
c  Begin Code
196
c  Begin Code
945
@@ -215,11 +213,6 @@
946
215
         write(*,'(a,12e10.3)') ' Masses:',(m(i),i=1,nparticles)
213
         write(*,'(a,12e10.3)') ' Masses:',(m(i),i=1,nparticles)
947
216
      endif                          !First_time
214
      endif                          !First_time
948
217
215
949
218
      if (.false.) then
950
219
         iconfig = isym(jfig)
951
220
         jfig = jfig+1
952
221
         if (jfig .gt. isym(0)) jfig=1      
953
222
      endif
954
223
      this_config = iconfig             !Pass iconfig to amplitude routine
216
      this_config = iconfig             !Pass iconfig to amplitude routine
955
224
C
217
C
956
225
C     Get fraction of beam energy if pdf's are used
218
C     Get fraction of beam energy if pdf's are used
957
@@ -279,19 +272,18 @@
958
279
         cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
272
         cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
959
280
         set_cm_rap=.true.
273
         set_cm_rap=.true.
960
281
c        Set shat
274
c        Set shat
963
282
         s(-nbranch) = m2**2+2*xbk(1)*ebeam(1) *
275
         s(-nbranch) = x(ndim)*stot         
962
283
     $                 (ebeam(2)+sqrt(ebeam(2)**2-m2**2))
964
284
      elseif (abs(lpp(2)) .ge. 1) then
276
      elseif (abs(lpp(2)) .ge. 1) then
965
285
         call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
277
         call sample_get_x(sjac,x(ndim),ndim,mincfig,0d0,1d0)
966
286
         xbk(2) = x(ndim)
278
         xbk(2) = x(ndim)
967
279
968
287
c        Set CM rapidity for use in the rap() function
280
c        Set CM rapidity for use in the rap() function
969
288
         p0=ebeam(1)+xbk(2)*ebeam(2)
281
         p0=ebeam(1)+xbk(2)*ebeam(2)
970
289
         p3=sqrt(ebeam(1)**2-m1**2)-xbk(2)*ebeam(2)
282
         p3=sqrt(ebeam(1)**2-m1**2)-xbk(2)*ebeam(2)
971
290
         cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
283
         cm_rap=.5d0*dlog((p0+p3)/(p0-p3))
972
291
         set_cm_rap=.true.
284
         set_cm_rap=.true.
973
292
c        Set shat
285
c        Set shat
976
293
         s(-nbranch) = m1**2+2*(ebeam(1)+sqrt(ebeam(1)**2-m1**2))
286
         s(-nbranch) =  x(ndim)*stot         
975
294
     $                 * xbk(2)*ebeam(2)
977
295
      else
287
      else
978
296
c        Set CM rapidity for use in the rap() function
288
c        Set CM rapidity for use in the rap() function
979
297
         p0=ebeam(1) + ebeam(2)
289
         p0=ebeam(1) + ebeam(2)
980
@@ -576,7 +568,7 @@
981
576
            if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
568
            if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
982
577
            if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
569
            if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
983
578
            if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
570
            if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
985
579
            if (mass_ion(2).ge.0d0) m1 = mass_ion(2)
571
            if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
986
580
            if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
572
            if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
987
581
            if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
573
            if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
988
582
            pi1(0)=ebeam(1)
574
            pi1(0)=ebeam(1)
989
@@ -598,7 +590,6 @@
990
598
         maxcfig=iconfig
590
         maxcfig=iconfig
991
599
         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
591
         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
992
600
         maxwgt=0d0
592
         maxwgt=0d0
993
601
c         write(*,'(a,12i4)') 'Summing configs',(isym(i),i=1,isym(0))
994
602
         nparticles   = nexternal
593
         nparticles   = nexternal
995
603
         nfinal       = nparticles-nincoming
594
         nfinal       = nparticles-nincoming
996
604
         nbranch      = nparticles-2
595
         nbranch      = nparticles-2
997
@@ -617,8 +608,6 @@
998
617
c
608
c
999
618
c         do i=1,mapconfig(0)
609
c         do i=1,mapconfig(0)
1000
619
         do i=mincfig,maxcfig
610
         do i=mincfig,maxcfig
1001
620
c         do k=1,isym(0)
1002
621
c            i = isym(k)
1003
622
            write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)
611
            write(*,'(15i4)') i,(minvar(j,i),j=1,ndim)
1004
623
            do j=1,ndim
612
            do j=1,ndim
1005
624
               ipole = minvar(j,i)
613
               ipole = minvar(j,i)
1006
@@ -663,7 +652,7 @@
1007
663
c      double precision spole(-max_branch:0),swidth(-max_branch:0)
652
c      double precision spole(-max_branch:0),swidth(-max_branch:0)
1008
664
      double precision jac,pswgt
653
      double precision jac,pswgt
1009
665
      integer nbranch
654
      integer nbranch
1011
666
      double precision x(40) ! ja 3/2/11 21->40 after strange segfault
655
      double precision x(*) ! ja 3/2/11 21->40 after strange segfault
1012
667
c
656
c
1013
668
c     Local
657
c     Local
1014
669
c
658
c
1015
670
659
1016
=== modified file 'Template/LO/SubProcesses/reweight.f'
1017
--- Template/LO/SubProcesses/reweight.f	2018-03-11 12:29:51 +0000
1018
+++ Template/LO/SubProcesses/reweight.f	2018-12-20 08:18:55 +0000
1019
@@ -1,3 +1,21 @@
1020
1
c     for cross-checking change in this file. 
1021
2
c     here is a minimal list of process that we have to 
1022
3
c     test
1023
4
c    
1024
5
c     SM
1025
6
c     ----
1026
7
c     p p > t t~ (up to 2jet)
1027
8
c     p p > w+ (up to 3 jet)
1028
9
c     p p > j j w+ (ordering seems important) 
1029
10
c     p p > z t t~ j j  (no MLM needed)
1030
11
c
1031
12
c
1032
13
c     HEFT
1033
14
c     ----
1034
15
c     p p > h j b b~
1035
16
c     q q > a a g q q
1036
17
c     g g > h g q q
1037
18
1038
1
      double precision function gamma(q0)
19
      double precision function gamma(q0)
1039
2
c**************************************************
20
c**************************************************
1040
3
c   calculates the branching probability
21
c   calculates the branching probability
1041
@@ -226,13 +244,15 @@
1042
226
      integer iddgluon, iddother, idgluon, idother
244
      integer iddgluon, iddother, idgluon, idother
1043
227
      logical isqcd
245
      logical isqcd
1044
228
      external isqcd
246
      external isqcd
1045
247
      integer get_color
1046
248
      external get_color 
1047
229
249
1048
230
      idmo=ipdg(imo)
250
      idmo=ipdg(imo)
1049
231
      idda1=ipdg(ida1)
251
      idda1=ipdg(ida1)
1050
232
      idda2=ipdg(ida2)
252
      idda2=ipdg(ida2)
1051
233
253
1052
234
      if (btest(mlevel,4)) then
254
      if (btest(mlevel,4)) then
1054
235
        write(*,*) ' updating ipart for: ',ida1,ida2,' -> ',imo
255
        write(*,*) 'updating ipart for: ',ida1,ida2,' -> ',imo
1055
236
      endif
256
      endif
1056
237
257
1057
238
      if (btest(mlevel,4)) then
258
      if (btest(mlevel,4)) then
1058
@@ -274,7 +294,6 @@
1059
274
     $        ' (',ipdg(imo),')'
294
     $        ' (',ipdg(imo),')'
1060
275
         return
295
         return
1061
276
      endif        
296
      endif        
1062
277
1063
278
c     FS clustering
297
c     FS clustering
1064
279
c     Transmit parton PDG code for parton vertex
298
c     Transmit parton PDG code for parton vertex
1065
280
      if(isjet(idmo)) then
299
      if(isjet(idmo)) then
1066
@@ -353,14 +372,26 @@
1067
353
c     quark -> gluon-quark or Z-quark or h-quark or W-quark
372
c     quark -> gluon-quark or Z-quark or h-quark or W-quark
1068
354
        ipart(1,imo)=ipart(1,ida2)
373
        ipart(1,imo)=ipart(1,ida2)
1069
355
        ipart(2,imo)=0
374
        ipart(2,imo)=0
1071
356
      else
375
      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda1)).eq.3.and.get_color(idda2).eq.1) then
1072
376
c       exotic q > q' Scalar         
1073
377
        ipart(1,imo)=ipart(1,ida1)
1074
378
        ipart(2,imo)=0
1075
379
      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda2)).eq.3.and.get_color(idda1).eq.1) then
1076
380
c       exotic q > Scalar q'
1077
381
        ipart(1,imo)=ipart(1,ida2)
1078
382
        ipart(2,imo)=0
1079
383
      else if (get_color(idmo).eq.1) then
1080
357
c     Color singlet
384
c     Color singlet
1081
358
         ipart(1,imo)=ipart(1,ida1)
385
         ipart(1,imo)=ipart(1,ida1)
1082
359
         ipart(2,imo)=ipart(1,ida2)
386
         ipart(2,imo)=ipart(1,ida2)
1083
387
      else
1084
388
         write(*,*) idmo,'>', idda1, idda2, 'color', get_color(idmo),'>', get_color(idda1), get_color(idda2)
1085
389
         write(*,*) "failed for ipartupdate. Please retry without MLM/default dynamical scale"
1086
390
         stop 3
1087
360
      endif
391
      endif
1088
361
      
392
      
1089
362
      if (btest(mlevel,4)) then
393
      if (btest(mlevel,4)) then
1091
363
        write(*,*) ' -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
394
        write(*,*) 'XY -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
1092
364
      endif
395
      endif
1093
365
396
1094
366
      return
397
      return
1095
@@ -522,7 +553,8 @@
1096
522
      common/to_stot/stot,m1,m2
553
      common/to_stot/stot,m1,m2
1097
523
554
1098
524
C   local variables
555
C   local variables
1100
525
      integer i, j, idi, idj, k
556
      integer i, j, idi, idj, k,m
1101
557
      integer get_color
1102
526
      real*8 PI
558
      real*8 PI
1103
527
      parameter( PI = 3.14159265358979323846d0 )
559
      parameter( PI = 3.14159265358979323846d0 )
1104
528
      integer iforest(2,-max_branch:-1,lmaxconfigs)
560
      integer iforest(2,-max_branch:-1,lmaxconfigs)
1105
@@ -545,6 +577,7 @@
1106
545
      logical chclusold,fail,increasecode
577
      logical chclusold,fail,increasecode
1107
546
      save chclusold
578
      save chclusold
1108
547
      integer tmpindex
579
      integer tmpindex
1109
580
      integer pdgm, pdgid1, pdgid2
1110
548
581
1111
549
      logical isqcd,isjet,isparton,cluster,isjetvx,is_octet
582
      logical isqcd,isjet,isparton,cluster,isjetvx,is_octet
1112
550
      integer ifsno
583
      integer ifsno
1113
@@ -676,6 +709,7 @@
1114
676
c     increasecode gives whether we should increase jcode at next vertex
709
c     increasecode gives whether we should increase jcode at next vertex
1115
677
      increasecode=.false.
710
      increasecode=.false.
1116
678
      do n=1,nexternal-2
711
      do n=1,nexternal-2
1117
712
c         write(*,*) 'QCD jet status (before n= ',n,'):',(iqjets(i),i=3,nexternal)
1118
679
        do i=1,2 ! index of the child in the interaction
713
        do i=1,2 ! index of the child in the interaction
1119
680
          do j=1,2 ! j index of the beam
714
          do j=1,2 ! j index of the beam
1120
681
            if(idacl(n,i).eq.ibeam(j))then
715
            if(idacl(n,i).eq.ibeam(j))then
1121
@@ -746,8 +780,16 @@
1122
746
c          The ishft gives the FS particle corresponding to imocl
780
c          The ishft gives the FS particle corresponding to imocl
1123
747
              if(.not.is_octet(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc)))then
781
              if(.not.is_octet(ipdgcl(ishft(1,ipart(1,imocl(n))-1),igraphs(1),iproc)))then
1124
748
                 ! split case for q a > q and for g > g h (with the gluon splitting into quark)
782
                 ! split case for q a > q and for g > g h (with the gluon splitting into quark)
1127
749
                 if (ipart(2,imocl(n)).eq.0) then ! q a > q case
783
                 ! also check for case of three scalar interaction (then do nothing)
1128
750
                    iqjets(ipart(1,imocl(n)))=0
784
                 pdgm   = ipdgcl(imocl(n),igraphs(1),iproc)
1129
785
                 pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
1130
786
                 pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
1131
787
                 if (.not.isqcd(pdgm).and..not.isqcd(pdgid1).and..not.isqcd(pdgid2)) then
1132
788
                    ! this is to avoid to do weird stuff for w+ w- z (or h h h) 
1133
789
                    ! this fix an issue for qq_zttxqq G1594.08
1134
790
                     continue
1135
791
                 elseif (ipart(2,imocl(n)).eq.0) then ! q a > q case
1136
792
                       iqjets(ipart(1,imocl(n)))=0
1137
751
                 else ! octet. want to be sure that both are tagged as jet before removing one
793
                 else ! octet. want to be sure that both are tagged as jet before removing one
1138
752
                    ! this prevent that both are removed in case of g > g h , g > q1 q2, q1 > a q1.
794
                    ! this prevent that both are removed in case of g > g h , g > q1 q2, q1 > a q1.
1139
753
                    ! at least one of the two should be kept as jet
795
                    ! at least one of the two should be kept as jet
1140
@@ -791,7 +833,57 @@
1141
791
           goodjet(imocl(n))=
833
           goodjet(imocl(n))=
1142
792
     $          (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.
834
     $          (isjet(ipdgcl(imocl(n),igraphs(1),iproc)).and.
1143
793
     $          goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
835
     $          goodjet(idacl(n,1)).and.goodjet(idacl(n,2)))
1145
794
        endif
836
            
1146
837
c       check case with g > g g
1147
838
c       where the hardest gluon is not goodjet but the other is.
1148
839
c       in that case change ipart(1,) of the mother gluon
1149
840
c             pure QCD jet
1150
841
c             need to take care of the following case:
1151
842
c                                      tttttttttt
1152
843
c                      gggggggggggggggg  
1153
844
c             ggggggggg                tttttttttt 
1154
845
c                      gggggggg
1155
846
c                        
1156
847
c             in that case the up gluon can be tag as the hardest one
1157
848
c             but this one is also lead to no QCD one.
1158
849
c             so in that case we have to change ipart(1) to the sofest gluon
1159
850
           pdgm   = ipdgcl(imocl(n),igraphs(1),iproc)
1160
851
           pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
1161
852
           pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
1162
853
           if (is_octet(pdgm).and.is_octet(pdgid1).and.is_octet(pdgid2))then
1163
854
c              write(*,*) 'pure QCD vertex (2)'
1164
855
c              write(*,*) pdgm , '>', pdgid1,' ', pdgid2
1165
856
c              write(*,*) 'ipart', ipart(1,imocl(n)), ipart(2,imocl(n))
1166
857
c              write(*,*) 'id', imocl(n), idacl(n,1),idacl(n,2)
1167
858
c              write(*,*) 'ipart',ipart(1,imocl(n)),'/', ipart(2,imocl(n)), ipart(1,idacl(n,1)),'/', ipart(2,idacl(n,1)),
1168
859
c     $             ipart(1,idacl(n,2)),'/', ipart(2,idacl(n,2))
1169
860
c              write(*,*) 'googjet', goodjet(idacl(n,1)),goodjet(idacl(n,2))
1170
861
              if (ipart(1,imocl(n)).eq.ipart(1, idacl(n,1))) then
1171
862
                  if (.not.goodjet(idacl(n,1)).and.goodjet(idacl(n,2))) then
1172
863
c                     write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,2)), ipart(2,idacl(n,2))
1173
864
                     do m =n_max_cl,n,-1
1174
865
                        if(ipart(1,m).eq.ipart(1,imocl(n)).and.ipart(2,m).eq.ipart(2,imocl(n)))then
1175
866
                           ipart(1,m) = ipart(1,idacl(n,2))
1176
867
                           ipart(2,m) = ipart(2,idacl(n,2))
1177
868
                        endif
1178
869
                     enddo
1179
870
                  endif
1180
871
               else
1181
872
                  if (.not.goodjet(idacl(n,2)).and.goodjet(idacl(n,1))) then
1182
873
c                     write(*,*) 'ggg with hard jet set a QED the second jet lead to', ipart(1,idacl(n,1)), ipart(2,idacl(n,1))
1183
874
                     do m =n_max_cl,n,-1
1184
875
                        if(ipart(1,m).eq.ipart(1,imocl(n)).and.ipart(2,m).eq.ipart(2,imocl(n)))then
1185
876
                           ipart(1,m) = ipart(1,idacl(n,1))
1186
877
                           ipart(2,m) = ipart(2,idacl(n,1))
1187
878
                        endif
1188
879
                     enddo
1189
880
                  endif
1190
881
               endif
1191
882
            endif
1192
883
1193
884
1194
885
1195
886
       endif
1196
795
      enddo
887
      enddo
1197
796
888
1198
797
      if (btest(mlevel,4))then
889
      if (btest(mlevel,4))then
1199
@@ -836,8 +928,8 @@
1200
836
            endif
928
            endif
1201
837
         enddo
929
         enddo
1202
838
         njetstore(iconfig)=njets
930
         njetstore(iconfig)=njets
1205
839
         if (btest(mlevel,4))
931
        if (btest(mlevel,4))
1206
840
     $        write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
932
     $       write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
1207
841
c     Recluster without requiring chcluster
933
c     Recluster without requiring chcluster
1208
842
         goto 100
934
         goto 100
1209
843
      else
935
      else
1210
844
936
1211
=== modified file 'Template/LO/SubProcesses/unwgt.f'
1212
--- Template/LO/SubProcesses/unwgt.f	2017-06-22 08:03:29 +0000
1213
+++ Template/LO/SubProcesses/unwgt.f	2018-12-20 08:18:55 +0000
1214
@@ -513,6 +513,10 @@
1215
513
513
1216
514
      double precision stot,m1,m2
514
      double precision stot,m1,m2
1217
515
      common/to_stot/stot,m1,m2
515
      common/to_stot/stot,m1,m2
1218
516
1219
517
      INTEGER IMIRROR
1220
518
      INTEGER IPROC
1221
519
      COMMON/TO_MIRROR/IMIRROR, IPROC
1222
516
c
520
c
1223
517
c     Data
521
c     Data
1224
518
c
522
c
1225
@@ -641,6 +645,13 @@
1226
641
            pb(4,isym(j,jsym))=pmass(j)
645
            pb(4,isym(j,jsym))=pmass(j)
1227
642
         enddo
646
         enddo
1228
643
      endif
647
      endif
1229
648
1230
649
      if (IMIRROR.eq.2.and.pmass(1).ne.pmass(2)) then
1231
650
c        Note that in this context isym(1,jsym) should never be "2" since the mass differ 
1232
651
         pb(4,isym(1,jsym))=pmass(2)
1233
652
         pb(4,isym(2,jsym))=pmass(1)
1234
653
      endif
1235
654
1236
644
c
655
c
1237
645
c     Add info on resonant mothers
656
c     Add info on resonant mothers
1238
646
c
657
c
1239
647
658
1240
=== added file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f'
1241
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	1970-01-01 00:00:00 +0000
1242
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	2018-12-20 08:18:55 +0000
1243
@@ -0,0 +1,880 @@
1244
1
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1245
2
      subroutine analysis_begin(nwgt,weights_info)
1246
3
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1247
4
      implicit none
1248
5
      integer nwgt
1249
6
      character*(*) weights_info(*)
1250
7
      integer i,l
1251
8
      character*6 cc(2)
1252
9
      data cc/'|T@NLO','|T@LO '/
1253
10
1254
11
c     Also specific perturbative orders can be directly plotted, adding for examples
1255
12
c     the following further entries in the variable data
1256
13
c     $        ' |T@QCD4QED0 ',
1257
14
c     $        ,' |T@QCD2QED2 ',' |T@QCD0QED4 ',' |T@QCD6QED0 '
1258
15
c     $        ,' |T@QCD4QED2 ',' |T@QCD2QED4 ',' |T@QCD0QED6 '
1259
16
c     $        ,' |T@QCD8QED0 ',' |T@QCD6QED2 ',' |T@QCD4QED4 '
1260
17
c     $        ,' |T@QCD2QED6 ',' |T@QCD0QED8 '
1261
18
c     
1262
19
c     See also line 376 in this file 
1263
20
      call HwU_inithist(nwgt,weights_info)
1264
21
      do i=1,2
1265
22
      l=(i-1)*55
1266
23
c transverse momenta
1267
24
      call HwU_book(l+ 1,'total rate                 '//cc(i),1,0.5d0,1.5d0)
1268
25
      call HwU_book(l+ 2,'1st charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
1269
26
      call HwU_book(l+ 3,'2nd charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
1270
27
      call HwU_book(l+ 4,'3rd charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
1271
28
      call HwU_book(l+ 5,'4th charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
1272
29
      call HwU_book(l+ 6,'electron log pT            '//cc(i),25,-0.2d0,3.8d0)
1273
30
      call HwU_book(l+ 7,'positron log pT            '//cc(i),25,-0.2d0,3.8d0)
1274
31
      call HwU_book(l+ 8,'mu-plus log pT             '//cc(i),25,-0.2d0,3.8d0)
1275
32
      call HwU_book(l+ 9,'mu-minus log pT            '//cc(i),25,-0.2d0,3.8d0)
1276
33
      call HwU_book(l+14,'epem log pT                '//cc(i),25,-0.2d0,3.8d0)
1277
34
      call HwU_book(l+15,'mupmum log pT              '//cc(i),25,-0.2d0,3.8d0)
1278
35
      call HwU_book(l+16,'1st OSSF lep pair log pT   '//cc(i),25,-0.2d0,3.8d0)
1279
36
      call HwU_book(l+17,'2nd OSSF lep pair log pT   '//cc(i),25,-0.2d0,3.8d0)
1280
37
      call HwU_book(l+18,'epve log pT                '//cc(i),25,-0.2d0,3.8d0)
1281
38
      call HwU_book(l+19,'mumvm log pT               '//cc(i),25,-0.2d0,3.8d0)
1282
39
      call HwU_book(l+20,'1st SF lep-neu pair log pT '//cc(i),25,-0.2d0,3.8d0)
1283
40
      call HwU_book(l+21,'2nd SF lep-neu pair log pT '//cc(i),25,-0.2d0,3.8d0)
1284
41
      call HwU_book(l+22,'epemmupmum log pT          '//cc(i),25,-0.2d0,3.8d0)
1285
42
      call HwU_book(l+23,'epvemumvm log pT           '//cc(i),25,-0.2d0,3.8d0)
1286
43
      call HwU_book(l+24,'top quark log pT           '//cc(i),25,-0.2d0,3.8d0)
1287
44
      call HwU_book(l+25,'anti-top quark log pT      '//cc(i),25,-0.2d0,3.8d0)
1288
45
      call HwU_book(l+26,'ttbar log pT               '//cc(i),25,-0.2d0,3.8d0)
1289
46
      call HwU_book(l+27,'log missing Et (neutrinos) '//cc(i),25,-0.2d0,3.8d0)
1290
47
      call HwU_book(l+28,'1st higgs log pT           '//cc(i),25,-0.2d0,3.8d0)
1291
48
      call HwU_book(l+29,'2nd higgs log pT           '//cc(i),25,-0.2d0,3.8d0)
1292
49
      call HwU_book(l+30,'higgs-pair log pT          '//cc(i),25,-0.2d0,3.8d0)
1293
50
      call HwU_book(l+31,'1st V-boson log pT         '//cc(i),25,-0.2d0,3.8d0)
1294
51
      call HwU_book(l+32,'2nd V-boson log pT         '//cc(i),25,-0.2d0,3.8d0)
1295
52
      call HwU_book(l+33,'3rd V-boson log pT         '//cc(i),25,-0.2d0,3.8d0)
1296
53
      call HwU_book(l+34,'1st jet log pT             '//cc(i),25,-0.2d0,3.8d0)
1297
54
      call HwU_book(l+35,'2nd jet log pT             '//cc(i),25,-0.2d0,3.8d0)
1298
55
      call HwU_book(l+36,'3rd jet log pT             '//cc(i),25,-0.2d0,3.8d0)
1299
56
c
1300
57
c invariant masses
1301
58
      call HwU_book(l+37,'epem log M                 '//cc(i),25,-0.2d0,3.8d0)
1302
59
      call HwU_book(l+38,'mupmum log M               '//cc(i),25,-0.2d0,3.8d0)
1303
60
      call HwU_book(l+39,'1st OSSF lep pair log M    '//cc(i),25,-0.2d0,3.8d0)
1304
61
      call HwU_book(l+40,'2nd OSSF lep pair log M    '//cc(i),25,-0.2d0,3.8d0)
1305
62
      call HwU_book(l+41,'epve log M                 '//cc(i),25,-0.2d0,3.8d0)
1306
63
      call HwU_book(l+42,'mumvm log M                '//cc(i),25,-0.2d0,3.8d0)
1307
64
      call HwU_book(l+43,'1st SF lep-neu pair log M  '//cc(i),25,-0.2d0,3.8d0)
1308
65
      call HwU_book(l+44,'2nd SF lep-neu pair log M  '//cc(i),25,-0.2d0,3.8d0)
1309
66
      call HwU_book(l+45,'epemmupmum log M           '//cc(i),25,-0.2d0,3.8d0)
1310
67
      call HwU_book(l+46,'epvemumvm log M            '//cc(i),25,-0.2d0,3.8d0)
1311
68
      call HwU_book(l+47,'ttbar log M                '//cc(i),25,-0.2d0,3.8d0)
1312
69
      call HwU_book(l+48,'higgs-pair log M           '//cc(i),25,-0.2d0,3.8d0)
1313
70
      call HwU_book(l+49,'j1j2 log M                 '//cc(i),25,-0.2d0,3.8d0)
1314
71
      call HwU_book(l+50,'j2j3 log M                 '//cc(i),25,-0.2d0,3.8d0)
1315
72
      call HwU_book(l+51,'j1j3 log M                 '//cc(i),25,-0.2d0,3.8d0)
1316
73
      call HwU_book(l+52,'j1j2j3 log M               '//cc(i),25,-0.2d0,3.8d0)
1317
74
      call HwU_book(l+53,'3w log M                   '//cc(i),25,-0.2d0,3.8d0)
1318
75
c
1319
76
c HT
1320
77
      call HwU_book(l+54,'log HT (partons)           '//cc(i),25,-0.2d0,3.8d0)
1321
78
      call HwU_book(l+55,'log HT (reconstructed)     '//cc(i),25,-0.2d0,3.8d0)
1322
79
      enddo
1323
80
      return
1324
81
      end
1325
82
1326
83
1327
84
1328
85
1329
86
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1330
87
      subroutine analysis_end(dummy)
1331
88
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1332
89
      implicit none
1333
90
      double precision dummy
1334
91
      call HwU_write_file
1335
92
      return                
1336
93
      end
1337
94
1338
95
1339
96
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1340
97
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
1341
98
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1342
99
      implicit none
1343
100
      include 'nexternal.inc'
1344
101
      include 'cuts.inc'
1345
102
      integer istatus(nexternal)
1346
103
      integer iPDG(nexternal)
1347
104
      double precision p(0:4,nexternal)
1348
105
      double precision wgts(*)
1349
106
      integer ibody
1350
107
      integer i,j,k,kk,l
1351
108
      double precision www,pQCD(0:3,nexternal),palg,rfj,sycut,yjmax
1352
109
     $     ,pjet(0:3,nexternal),tmp,ptlep(4),ptem,ptep,ptmm,ptmp,ptepem
1353
110
     $     ,ptmpmm,ptepve,ptmmvm,ptepemmpmm,ptepvemmvm,ptt,ptat,pttt
1354
111
     $     ,etmiss,pth(2),pthh,ptv(3),ptjet(3),Mepem,Mmpmm,Mepve,Mmmvm
1355
112
     $     ,Mepemmpmm,Mepvemmvm,Mtt,Mhh,Mj1j2,Mj1j3,Mj2j3,Mj1j2j3,Mvvv
1356
113
     $     ,HTparton,HTreco,p_reco(0:4,nexternal)
1357
114
      integer nQCD,jet(nexternal),njet,itop,iatop,iem,iep,imp,imm,ive
1358
115
     $     ,ivm,iv1,iv2,iv3,ih1,ih2,il,ipdg_reco(nexternal)
1359
116
      double precision getptv4,getptv4_2,getptv4_4,getinvm4_2
1360
117
     $     ,getinvm4_3,getinvm4_4,l10
1361
118
      external getptv4,getptv4_2,getptv4_4,getinvm4_2,getinvm4_3
1362
119
     $     ,getinvm4_4,l10
1363
120
      integer orders_tag_plot
1364
121
      common /corderstagplot/ orders_tag_plot
1365
122
c First, try to recombine photons with leptons      
1366
123
      if (.not.quarkphreco) then
1367
124
         write (*,*) 'quark-photon recombination is turned off. '/
1368
125
     $        /'Do need it'
1369
126
         stop
1370
127
      endif
1371
128
      if (.not. lepphreco) then
1372
129
         write (*,*) 'lepton-photon recombination is turned off. '
1373
130
     $        //'Do need it.'
1374
131
         stop
1375
132
      endif
1376
133
      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
1377
134
     $                       p, iPDG, p_reco, iPDG_reco)
1378
135
1379
136
c Put all (light) QCD partons(+photon) in momentum array for jet clustering.
1380
137
      nQCD=0
1381
138
      do j=nincoming+1,nexternal
1382
139
         if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or.
1383
140
     $        ipdg_reco(j).eq.22)then
1384
141
            nQCD=nQCD+1
1385
142
            do i=0,3
1386
143
               pQCD(i,nQCD)=p_reco(i,j)
1387
144
            enddo
1388
145
         endif
1389
146
      enddo
1390
147
      
1391
148
C---CLUSTER THE EVENT
1392
149
      palg  = jetalgo
1393
150
      rfj   = jetradius
1394
151
      sycut = ptj
1395
152
      yjmax = etaj
1396
153
      
1397
154
c******************************************************************************
1398
155
c     call FASTJET to get all the jets
1399
156
c
1400
157
c     INPUT:
1401
158
c     input momenta:               pQCD(0:3,nexternal), energy is 0th component
1402
159
c     number of input momenta:     nQCD
1403
160
c     radius parameter:            rfj
1404
161
c     minumum jet pt:              sycut
1405
162
c     jet algorithm:               palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
1406
163
c
1407
164
c     OUTPUT:
1408
165
c     jet momenta:                             pjet(0:3,nexternal), E is 0th cmpnt
1409
166
c     the number of jets (with pt > SYCUT):    njet
1410
167
c     the jet for a given particle 'i':        jet(i),   note that this is
1411
168
c     the particle in pQCD, which doesn't necessarily correspond to the particle
1412
169
c     label in the process
1413
170
c
1414
171
      call amcatnlo_fastjetppgenkt_etamax(pQCD,nQCD,rfj,sycut,yjmax,palg
1415
172
     $     ,pjet,njet,jet)
1416
173
c
1417
174
c******************************************************************************
1418
175
1419
176
c Look for the other physics objects
1420
177
      itop=0
1421
178
      iatop=0
1422
179
      iem=0
1423
180
      iep=0
1424
181
      iem=0
1425
182
      iep=0
1426
183
      imp=0
1427
184
      imm=0
1428
185
      ive=0
1429
186
      ivm=0
1430
187
      iv1=0
1431
188
      iv2=0
1432
189
      iv3=0
1433
190
      ih1=0
1434
191
      ih2=0
1435
192
      do i=1,nexternal
1436
193
         if (ipdg_reco(i).eq.6) then
1437
194
            itop=i
1438
195
         elseif(ipdg_reco(i).eq.-6) then
1439
196
            iatop=i
1440
197
         elseif(ipdg_reco(i).eq.11) then
1441
198
            iem=i
1442
199
         elseif(ipdg_reco(i).eq.13) then
1443
200
            imm=i
1444
201
         elseif(ipdg_reco(i).eq.-11) then
1445
202
            iep=i
1446
203
         elseif(ipdg_reco(i).eq.-13) then
1447
204
            imp=i
1448
205
         elseif(abs(ipdg_reco(i)).eq.12) then
1449
206
            ive=i
1450
207
         elseif(abs(ipdg_reco(i)).eq.14) then
1451
208
            ivm=i
1452
209
         elseif(abs(ipdg_reco(i)).eq.24 .or. ipdg_reco(i).eq.23) then
1453
210
            if (iv1.eq.0) then
1454
211
               iv1=i
1455
212
            else
1456
213
               if (iv2.eq.0) then
1457
214
                  iv2=i
1458
215
               else
1459
216
                  if(iv3.eq.0) then
1460
217
                     iv3=i
1461
218
                  else
1462
219
                     write (*,*) 'too many vector bosons'
1463
220
                     stop
1464
221
                  endif
1465
222
               endif
1466
223
            endif
1467
224
         elseif(abs(ipdg_reco(i)).eq.25) then
1468
225
            if (ih1.eq.0) then
1469
226
               ih1=i
1470
227
            else
1471
228
               if (ih2.eq.0) then
1472
229
                  ih2=i
1473
230
               else
1474
231
                  write (*,*) 'too many higgs bosons'
1475
232
                  stop
1476
233
               endif
1477
234
            endif
1478
235
         endif
1479
236
      enddo
1480
237
      if (itop.ne.0) ptt=getptv4(p_reco(0,itop))
1481
238
      if (iatop.ne.0) ptat=getptv4(p_reco(0,iatop))
1482
239
      if (itop.ne.0 .and. iatop.ne.0) then
1483
240
         pttt=getptv4_2(p_reco(0,itop),p_reco(0,iatop))
1484
241
         Mtt=getinvm4_2(p_reco(0,itop),p_reco(0,iatop))
1485
242
      endif
1486
243
      if (iem.ne.0) ptem=getptv4(p_reco(0,iem))
1487
244
      if (iep.ne.0) ptep=getptv4(p_reco(0,iep))
1488
245
      if (imm.ne.0) ptmm=getptv4(p_reco(0,imm))
1489
246
      if (imp.ne.0) ptmp=getptv4(p_reco(0,imp))
1490
247
      if (iem.ne.0 .and. iep.ne.0) then
1491
248
         ptepem=getptv4_2(p_reco(0,iem),p_reco(0,iep))
1492
249
         Mepem=getinvm4_2(p_reco(0,iem),p_reco(0,iep))
1493
250
      endif
1494
251
      if (imm.ne.0 .and. imp.ne.0) then
1495
252
         ptmpmm=getptv4_2(p_reco(0,imm),p_reco(0,imp))
1496
253
         Mmpmm=getinvm4_2(p_reco(0,imm),p_reco(0,imp))
1497
254
      endif
1498
255
      if (iep.ne.0 .and. ive.ne.0) then
1499
256
         ptepve=getptv4_2(p_reco(0,iep),p_reco(0,ive))
1500
257
         Mepve=getinvm4_2(p_reco(0,iep),p_reco(0,ive))
1501
258
      endif
1502
259
      if (imm.ne.0 .and. ivm.ne.0) then
1503
260
         ptmmvm=getptv4_2(p_reco(0,imm),p_reco(0,ivm))
1504
261
         Mmmvm=getinvm4_2(p_reco(0,imm),p_reco(0,ivm))
1505
262
      endif
1506
263
      if (iem.ne.0 .and. iep.ne.0 .and. imm.ne.0 .and. imp.ne.0) then
1507
264
         ptepemmpmm=getptv4_4(p_reco(0,iem),p_reco(0,iep),p_reco(0,imm),p_reco(0,imp))
1508
265
         Mepemmpmm=getinvm4_4(p_reco(0,iem),p_reco(0,iep),p_reco(0,imm),p_reco(0,imp))
1509
266
      endif      
1510
267
      if (ive.ne.0 .and. iep.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
1511
268
         ptepvemmvm=getptv4_4(p_reco(0,iep),p_reco(0,ive),p_reco(0,imm),p_reco(0,ivm))
1512
269
         Mepvemmvm=getinvm4_4(p_reco(0,iep),p_reco(0,ive),p_reco(0,imm),p_reco(0,ivm))
1513
270
      endif
1514
271
1515
272
      il=0
1516
273
      if (iem.ne.0) then
1517
274
         il=il+1
1518
275
         ptlep(il)=ptem
1519
276
      endif
1520
277
      if (iep.ne.0) then
1521
278
         il=il+1
1522
279
         ptlep(il)=ptep
1523
280
      endif
1524
281
      if (imm.ne.0) then
1525
282
         il=il+1
1526
283
         ptlep(il)=ptmm
1527
284
      endif
1528
285
      if (imp.ne.0) then
1529
286
         il=il+1
1530
287
         ptlep(il)=ptmp
1531
288
      endif
1532
289
      if (il.gt.1) then
1533
290
         do i=1,il-1
1534
291
            do j=1,il-i
1535
292
               if (ptlep(j).lt.ptlep(j+1)) then
1536
293
                  tmp=ptlep(j)
1537
294
                  ptlep(j)=ptlep(j+1)
1538
295
                  ptlep(j+1)=tmp
1539
296
               endif
1540
297
            enddo
1541
298
         enddo
1542
299
      endif
1543
300
1544
301
      
1545
302
c missing Et
1546
303
      if (ive.ne.0 .and. ivm.ne.0) then
1547
304
         etmiss=getptv4_2(p_reco(0,ivm),p_reco(0,ive))
1548
305
      elseif (ive.ne.0 .or. ivm.ne.0) then
1549
306
         etmiss=getptv4(p_reco(0,ive+ivm))
1550
307
      endif
1551
308
      
1552
309
      if (ih1.ne.0)pth(1)=getptv4(p_reco(0,ih1))
1553
310
      if (ih2.ne.0)pth(2)=getptv4(p_reco(0,ih2))
1554
311
      if (ih1.ne.0 .and. ih2.ne.0) then
1555
312
         pthh=getptv4_2(p_reco(0,ih1),p_reco(0,ih2))
1556
313
         Mhh=getinvm4_2(p_reco(0,ih1),p_reco(0,ih2))
1557
314
c     order the higgs bosons (if there are 2)
1558
315
         if (pth(1).lt.pth(2)) then
1559
316
            tmp=pth(1)
1560
317
            pth(1)=pth(2)
1561
318
            pth(2)=tmp
1562
319
         endif
1563
320
      endif
1564
321
      if (iv1.ne.0) ptv(1)=getptv4(p_reco(0,iv1))
1565
322
      if (iv2.ne.0) ptv(2)=getptv4(p_reco(0,iv2))
1566
323
      if (iv3.ne.0) ptv(3)=getptv4(p_reco(0,iv3))
1567
324
      if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
1568
325
         Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
1569
326
c     order the vector bosons (if there are 3)
1570
327
         do i=1,2
1571
328
            do j=1,3-i
1572
329
               if (ptv(j).lt.ptv(j+1)) then
1573
330
                  tmp=ptv(j)
1574
331
                  ptv(j)=ptv(j+1)
1575
332
                  ptv(j+1)=tmp
1576
333
               endif
1577
334
            enddo
1578
335
         enddo
1579
336
      elseif (iv1.ne.0 .and. iv2.ne.0) then
1580
337
c     order the vector bosons (if there are 2)
1581
338
         if (ptv(1).lt.ptv(2)) then
1582
339
            tmp=ptv(1)
1583
340
            ptv(1)=ptv(2)
1584
341
            ptv(2)=tmp
1585
342
         endif
1586
343
      endif
1587
344
      do i=1,njet
1588
345
         ptjet(i)=getptv4(pjet(0,i))
1589
346
      enddo
1590
347
      if (njet.ge.2) then
1591
348
         Mj1j2=getinvm4_2(pjet(0,1),pjet(0,2))
1592
349
      endif
1593
350
      if (njet.ge.3) then
1594
351
         Mj2j3=getinvm4_2(pjet(0,2),pjet(0,3))
1595
352
         Mj1j3=getinvm4_2(pjet(0,1),pjet(0,3))
1596
353
         Mj1j2j3=getinvm4_3(pjet(0,1),pjet(0,2),pjet(0,3))
1597
354
      endif
1598
355
c
1599
356
      HTparton=0d0
1600
357
      HTreco=0d0
1601
358
      do i=1,nexternal
1602
359
         HTparton=HTparton+getptv4(p(0,i))
1603
360
         if (abs(ipdg_reco(i)).gt.5 .and. ipdg_reco(i).ne.21 .and.
1604
361
     $        ipdg_reco(i).ne.22 .and. abs(ipdg_reco(i)).ne.12 .and.
1605
362
     $        abs(ipdg_reco(i)).ne.14 .and. abs(ipdg_reco(i)).ne.16)
1606
363
     $        then
1607
364
            HTreco=HTreco+getptv4(p_reco(0,i))
1608
365
         endif
1609
366
      enddo
1610
367
      do i=1,njet
1611
368
         HTreco=HTreco+getptv4(pjet(0,i))
1612
369
      enddo
1613
370
      if (ive.ne.0 .or. ivm.ne.0) HTreco=HTreco+etmiss
1614
371
         
1615
372
      do i=1,2
1616
373
         l=(i-1)*55
1617
374
         if (ibody.ne.3 .and.i.eq.2) cycle
1618
375
1619
376
c         How to tag orders (QCD+QED*100)
1620
377
c
1621
378
c         if (i.eq. 3.and.orders_tag_plot.ne.4) cycle
1622
379
c         if (i.eq. 4.and.orders_tag_plot.ne.202) cycle
1623
380
c         if (i.eq. 5.and.orders_tag_plot.ne.400) cycle
1624
381
c         if (i.eq. 6.and.orders_tag_plot.ne.6) cycle
1625
382
c         if (i.eq. 7.and.orders_tag_plot.ne.204) cycle
1626
383
c         if (i.eq. 8.and.orders_tag_plot.ne.402) cycle
1627
384
c         if (i.eq. 9.and.orders_tag_plot.ne.600) cycle
1628
385
c         if (i.eq.10.and.orders_tag_plot.ne.8) cycle
1629
386
c         if (i.eq.11.and.orders_tag_plot.ne.206) cycle
1630
387
c         if (i.eq.12.and.orders_tag_plot.ne.404) cycle
1631
388
c         if (i.eq.13.and.orders_tag_plot.ne.602) cycle
1632
389
c         if (i.eq.14.and.orders_tag_plot.ne.800) cycle
1633
390
1634
391
1635
392
1636
393
c total rate         
1637
394
         call HwU_fill(l+ 1,1d0,wgts)
1638
395
c transverse momenta         
1639
396
         if (il .ge.1) call HwU_fill(l+ 2,l10(ptlep(1)),wgts)
1640
397
         if (il .ge.2) call HwU_fill(l+ 3,l10(ptlep(2)),wgts)
1641
398
         if (il .ge.3) call HwU_fill(l+ 4,l10(ptlep(3)),wgts)
1642
399
         if (il .ge.4) call HwU_fill(l+ 5,l10(ptlep(4)),wgts)
1643
400
         if (iem.ne.0) call HwU_fill(l+ 6,l10(ptem),wgts)
1644
401
         if (iep.ne.0) call HwU_fill(l+ 7,l10(ptep),wgts)
1645
402
         if (imp.ne.0) call HwU_fill(l+ 8,l10(ptmp),wgts)
1646
403
         if (imm.ne.0) call HwU_fill(l+ 9,l10(ptmm),wgts)
1647
404
         if (iep.ne.0 .and. iem.ne.0) call HwU_fill(l+14,l10(ptepem),wgts)
1648
405
         if (imp.ne.0 .and. imm.ne.0) call HwU_fill(l+15,l10(ptmpmm),wgts)
1649
406
         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0) then
1650
407
            if (ptepem.gt.ptmpmm) then
1651
408
               call HwU_fill(l+16,l10(ptepem),wgts)
1652
409
               call HwU_fill(l+17,l10(ptmpmm),wgts)
1653
410
            else 
1654
411
              call HwU_fill(l+17,l10(ptepem),wgts)
1655
412
               call HwU_fill(l+16,l10(ptmpmm),wgts)
1656
413
            endif
1657
414
         elseif (iep.ne.0 .and. iem.ne.0) then
1658
415
            call HwU_fill(l+17,l10(ptepem),wgts)
1659
416
         elseif (imp.ne.0 .and. imm.ne.0) then
1660
417
            call HwU_fill(l+17,l10(ptmpmm),wgts)
1661
418
         endif
1662
419
         if (iep.ne.0 .and. ive.ne.0) call HwU_fill(l+18,l10(ptepve),wgts)
1663
420
         if (imm.ne.0 .and. ivm.ne.0) call HwU_fill(l+19,l10(ptmmvm),wgts)
1664
421
         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
1665
422
            if (ptepve.gt.ptmmvm) then
1666
423
               call HwU_fill(l+20,l10(ptepve),wgts)
1667
424
               call HwU_fill(l+21,l10(ptmmvm),wgts)
1668
425
            else
1669
426
               call HwU_fill(l+21,l10(ptepve),wgts)
1670
427
               call HwU_fill(l+20,l10(ptmmvm),wgts)
1671
428
            endif
1672
429
         elseif (iep.ne.0 .and. ive.ne.0) then
1673
430
            call HwU_fill(l+20,l10(ptepve),wgts)
1674
431
         elseif (imm.ne.0 .and. ivm.ne.0) then
1675
432
            call HwU_fill(l+20,l10(ptmmvm),wgts)
1676
433
         endif
1677
434
         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0)
1678
435
     $        call HwU_fill(l+22,l10(ptepemmpmm),wgts)
1679
436
         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0)
1680
437
     $        call HwU_fill(l+23,l10(ptepvemmvm),wgts)
1681
438
         if (itop.ne.0) call HwU_fill(l+24,l10(ptt),wgts)
1682
439
         if (iatop.ne.0) call HwU_fill(l+25,l10(ptat),wgts)
1683
440
         if (itop.ne.0 .and. iatop.ne.0) call HwU_fill(l+26,l10(pttt),wgts)
1684
441
         if (ive.ne.0 .or. ivm .ne.0) call HwU_fill(l+27,l10(etmiss),wgts)
1685
442
         if (ih1.ne.0) call HwU_fill(l+28,l10(pth(1)),wgts)
1686
443
         if (ih2.ne.0) call HwU_fill(l+29,l10(pth(2)),wgts)
1687
444
         if (ih1.ne.0 .and. ih2.ne.0) call HwU_fill(l+30,l10(pthh),wgts)
1688
445
         if (iv1.ne.0) call HwU_fill(l+31,l10(ptv(1)),wgts)
1689
446
         if (iv2.ne.0) call HwU_fill(l+32,l10(ptv(2)),wgts)
1690
447
         if (iv3.ne.0) call HwU_fill(l+33,l10(ptv(3)),wgts)
1691
448
         if (njet.ge.1) call HwU_fill(l+34,l10(ptjet(1)),wgts)
1692
449
         if (njet.ge.2) call HwU_fill(l+35,l10(ptjet(2)),wgts)
1693
450
         if (njet.ge.3) call HwU_fill(l+36,l10(ptjet(3)),wgts)
1694
451
c invariant masses
1695
452
         if (iep.ne.0 .and. iem.ne.0) call HwU_fill(l+37,l10(Mepem),wgts)
1696
453
         if (imp.ne.0 .and. imm.ne.0) call HwU_fill(l+38,l10(Mmpmm),wgts)
1697
454
         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0) then
1698
455
            if (Mepem.gt.Mmpmm) then
1699
456
               call HwU_fill(l+39,l10(Mepem),wgts)
1700
457
               call HwU_fill(l+40,l10(Mmpmm),wgts)
1701
458
            else
1702
459
               call HwU_fill(l+40,l10(Mepem),wgts)
1703
460
               call HwU_fill(l+39,l10(Mmpmm),wgts)
1704
461
            endif
1705
462
         elseif (iep.ne.0 .and. iem.ne.0) then
1706
463
            call HwU_fill(l+39,l10(Mepem),wgts)
1707
464
         elseif (imp.ne.0 .and. imm.ne.0) then
1708
465
            call HwU_fill(l+39,l10(Mmpmm),wgts)
1709
466
         endif
1710
467
         if (iep.ne.0 .and. ive.ne.0) call HwU_fill(l+41,l10(Mepve),wgts)
1711
468
         if (imm.ne.0 .and. ivm.ne.0) call HwU_fill(l+42,l10(Mmmvm),wgts)
1712
469
         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0) then
1713
470
            if (Mepve.gt.Mmmvm) then
1714
471
               call HwU_fill(l+43,l10(Mepve),wgts)
1715
472
               call HwU_fill(l+44,l10(Mmmvm),wgts)
1716
473
            else
1717
474
               call HwU_fill(l+44,l10(Mepve),wgts)
1718
475
               call HwU_fill(l+43,l10(Mmmvm),wgts)
1719
476
            endif
1720
477
         elseif (iep.ne.0 .and. ive.ne.0) then
1721
478
            call HwU_fill(l+43,l10(Mepve),wgts)
1722
479
         elseif (imm.ne.0 .and. ivm.ne.0) then
1723
480
            call HwU_fill(l+43,l10(Mmmvm),wgts)
1724
481
         endif
1725
482
         if (iep.ne.0 .and. iem.ne.0 .and. imp.ne.0 .and. imm.ne.0)
1726
483
     $        call HwU_fill(l+45,l10(Mepemmpmm),wgts)
1727
484
         if (iep.ne.0 .and. ive.ne.0 .and. imm.ne.0 .and. ivm.ne.0)
1728
485
     $        call HwU_fill(l+46,l10(Mepvemmvm),wgts)
1729
486
         if (itop.ne.0 .and. iatop.ne.0) call HwU_fill(l+47,l10(Mtt),wgts)
1730
487
         if (ih1.ne.0 .and. ih2.ne.0) call HwU_fill(l+48,l10(Mhh),wgts)
1731
488
         if (njet.ge.2) call HwU_fill(l+49,l10(Mj1j2),wgts)
1732
489
         if (njet.ge.3) call HwU_fill(l+50,l10(Mj2j3),wgts)
1733
490
         if (njet.ge.3) call HwU_fill(l+51,l10(Mj1j3),wgts)
1734
491
         if (njet.ge.3) call HwU_fill(l+52,l10(Mj1j2j3),wgts)
1735
492
         if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) call HwU_fill(l+53,l10(Mvvv),wgts)
1736
493
c HT
1737
494
         call HwU_fill(l+54,l10(HTparton),wgts)
1738
495
         call HwU_fill(l+55,l10(HTreco),wgts)
1739
496
      enddo
1740
497
1741
498
 999  return      
1742
499
      end
1743
500
1744
501
1745
502
      double precision function l10(var)
1746
503
      implicit none
1747
504
      double precision var
1748
505
      if (var.gt.0) then
1749
506
         l10=log10(var)
1750
507
      else
1751
508
         l10=-99d99
1752
509
      endif
1753
510
      return
1754
511
      end
1755
512
      
1756
513
      function getrapidity(en,pl)
1757
514
      implicit none
1758
515
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
1759
516
      parameter (tiny=1.d-8)
1760
517
c
1761
518
      xplus=en+pl
1762
519
      xminus=en-pl
1763
520
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
1764
521
        if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
1765
522
          y=0.5d0*log( xplus/xminus )
1766
523
        else
1767
524
          y=sign(1.d0,pl)*1.d8
1768
525
        endif
1769
526
      else
1770
527
        y=sign(1.d0,pl)*1.d8
1771
528
      endif
1772
529
      getrapidity=y
1773
530
      return
1774
531
      end
1775
532
1776
533
1777
534
      function getpseudorap(en,ptx,pty,pl)
1778
535
      implicit none
1779
536
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
1780
537
      parameter (tiny=1.d-5)
1781
538
c
1782
539
      pt=sqrt(ptx**2+pty**2)
1783
540
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
1784
541
        eta=sign(1.d0,pl)*1.d8
1785
542
      else
1786
543
        th=atan2(pt,pl)
1787
544
        eta=-log(tan(th/2.d0))
1788
545
      endif
1789
546
      getpseudorap=eta
1790
547
      return
1791
548
      end
1792
549
1793
550
1794
551
      function getinvm(en,ptx,pty,pl)
1795
552
      implicit none
1796
553
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
1797
554
      parameter (tiny=1.d-5)
1798
555
c
1799
556
      tmp=en**2-ptx**2-pty**2-pl**2
1800
557
      if(tmp.gt.0.d0)then
1801
558
        tmp=sqrt(tmp)
1802
559
      elseif(tmp.gt.-tiny)then
1803
560
        tmp=0.d0
1804
561
      else
1805
562
        write(*,*)'Attempt to compute a negative mass'
1806
563
        stop
1807
564
      endif
1808
565
      getinvm=tmp
1809
566
      return
1810
567
      end
1811
568
1812
569
1813
570
      function getdelphi(ptx1,pty1,ptx2,pty2)
1814
571
      implicit none
1815
572
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
1816
573
      parameter (tiny=1.d-5)
1817
574
c
1818
575
      pt1=sqrt(ptx1**2+pty1**2)
1819
576
      pt2=sqrt(ptx2**2+pty2**2)
1820
577
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
1821
578
        tmp=ptx1*ptx2+pty1*pty2
1822
579
        tmp=tmp/(pt1*pt2)
1823
580
        if(abs(tmp).gt.1.d0+tiny)then
1824
581
          write(*,*)'Cosine larger than 1'
1825
582
          stop
1826
583
        elseif(abs(tmp).ge.1.d0)then
1827
584
          tmp=sign(1.d0,tmp)
1828
585
        endif
1829
586
        tmp=acos(tmp)
1830
587
      else
1831
588
        tmp=1.d8
1832
589
      endif
1833
590
      getdelphi=tmp
1834
591
      return
1835
592
      end
1836
593
1837
594
1838
595
      function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
1839
596
      implicit none
1840
597
      real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
1841
598
     # getpseudorap,getdelphi
1842
599
c
1843
600
      deta=getpseudorap(en1,ptx1,pty1,pl1)-
1844
601
     #     getpseudorap(en2,ptx2,pty2,pl2)
1845
602
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
1846
603
      getdr=sqrt(dphi**2+deta**2)
1847
604
      return
1848
605
      end
1849
606
1850
607
1851
608
      function getdry(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
1852
609
      implicit none
1853
610
      real*8 getdry,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
1854
611
     # getrapidity,getdelphi
1855
612
c
1856
613
      deta=getrapidity(en1,pl1)-
1857
614
     #     getrapidity(en2,pl2)
1858
615
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
1859
616
      getdry=sqrt(dphi**2+deta**2)
1860
617
      return
1861
618
      end
1862
619
1863
620
1864
621
      function getptv(p)
1865
622
      implicit none
1866
623
      real*8 getptv,p(5)
1867
624
c
1868
625
      getptv=sqrt(p(1)**2+p(2)**2)
1869
626
      return
1870
627
      end
1871
628
1872
629
1873
630
      function getpseudorapv(p)
1874
631
      implicit none
1875
632
      real*8 getpseudorapv,p(5)
1876
633
      real*8 getpseudorap
1877
634
c
1878
635
      getpseudorapv=getpseudorap(p(4),p(1),p(2),p(3))
1879
636
      return
1880
637
      end
1881
638
1882
639
1883
640
      function getrapidityv(p)
1884
641
      implicit none
1885
642
      real*8 getrapidityv,p(5)
1886
643
      real*8 getrapidity
1887
644
c
1888
645
      getrapidityv=getrapidity(p(4),p(3))
1889
646
      return
1890
647
      end
1891
648
1892
649
1893
650
      function getdrv(p1,p2)
1894
651
      implicit none
1895
652
      real*8 getdrv,p1(5),p2(5)
1896
653
      real*8 getdr
1897
654
c
1898
655
      getdrv=getdr(p1(4),p1(1),p1(2),p1(3),
1899
656
     #             p2(4),p2(1),p2(2),p2(3))
1900
657
      return
1901
658
      end
1902
659
1903
660
1904
661
      function getinvmv(p)
1905
662
      implicit none
1906
663
      real*8 getinvmv,p(5)
1907
664
      real*8 getinvm
1908
665
c
1909
666
      getinvmv=getinvm(p(4),p(1),p(2),p(3))
1910
667
      return
1911
668
      end
1912
669
1913
670
1914
671
      function getdelphiv(p1,p2)
1915
672
      implicit none
1916
673
      real*8 getdelphiv,p1(5),p2(5)
1917
674
      real*8 getdelphi
1918
675
c
1919
676
      getdelphiv=getdelphi(p1(1),p1(2),
1920
677
     #                     p2(1),p2(2))
1921
678
      return
1922
679
      end
1923
680
1924
681
1925
682
      function getptv4(p)
1926
683
      implicit none
1927
684
      real*8 getptv4,p(0:3)
1928
685
c
1929
686
      getptv4=sqrt(max(p(1)**2+p(2)**2,0d0))
1930
687
      return
1931
688
      end
1932
689
1933
690
      function getptv4_2(p1,p2)
1934
691
      implicit none
1935
692
      real*8 getptv4_2,p1(0:3),p2(0:3),psum(0:3)
1936
693
      integer i
1937
694
      do i=0,3
1938
695
         psum(i)=p1(i)+p2(i)
1939
696
      enddo
1940
697
      getptv4_2=sqrt(max(psum(1)**2+psum(2)**2,0d0))
1941
698
      return
1942
699
      end
1943
700
1944
701
      function getptv4_4(p1,p2,p3,p4)
1945
702
      implicit none
1946
703
      real*8 getptv4_4,p1(0:3),p2(0:3),p3(0:3),p4(0:3),psum(0:3)
1947
704
      integer i
1948
705
      do i=0,3
1949
706
         psum(i)=p1(i)+p2(i)+p3(i)+p4(i)
1950
707
      enddo
1951
708
      getptv4_4=sqrt(max(psum(1)**2+psum(2)**2,0d0))
1952
709
      return
1953
710
      end
1954
711
1955
712
1956
713
      function getpseudorapv4(p)
1957
714
      implicit none
1958
715
      real*8 getpseudorapv4,p(0:3)
1959
716
      real*8 getpseudorap
1960
717
c
1961
718
      getpseudorapv4=getpseudorap(p(0),p(1),p(2),p(3))
1962
719
      return
1963
720
      end
1964
721
1965
722
1966
723
      function getrapidityv4(p)
1967
724
      implicit none
1968
725
      real*8 getrapidityv4,p(0:3)
1969
726
      real*8 getrapidity
1970
727
c
1971
728
      getrapidityv4=getrapidity(p(0),p(3))
1972
729
      return
1973
730
      end
1974
731
1975
732
1976
733
      function getdrv4(p1,p2)
1977
734
      implicit none
1978
735
      real*8 getdrv4,p1(0:3),p2(0:3)
1979
736
      real*8 getdr
1980
737
c
1981
738
      getdrv4=getdr(p1(0),p1(1),p1(2),p1(3),
1982
739
     #              p2(0),p2(1),p2(2),p2(3))
1983
740
      return
1984
741
      end
1985
742
1986
743
1987
744
      function getinvm4(p)
1988
745
      implicit none
1989
746
      real*8 getinvm4,p(0:3)
1990
747
      real*8 getinvm
1991
748
c
1992
749
      getinvm4=getinvm(p(0),p(1),p(2),p(3))
1993
750
      return
1994
751
      end
1995
752
1996
753
      function getinvm4_2(p1,p2)
1997
754
      implicit none
1998
755
      real*8 getinvm4_2,p1(0:3),p2(0:3),p(0:3)
1999
756
      real*8 getinvm
2000
757
      integer i
2001
758
      do i=0,3
2002
759
         p(i)=p1(i)+p2(i)
2003
760
      enddo
2004
761
      getinvm4_2=getinvm(p(0),p(1),p(2),p(3))
2005
762
      return
2006
763
      end
2007
764
2008
765
      function getinvm4_3(p1,p2,p3)
2009
766
      implicit none
2010
767
      real*8 getinvm4_3,p1(0:3),p2(0:3),p3(0:3),p(0:3)
2011
768
      real*8 getinvm
2012
769
      integer i
2013
770
      do i=0,3
2014
771
         p(i)=p1(i)+p2(i)+p3(i)
2015
772
      enddo
2016
773
      getinvm4_3=getinvm(p(0),p(1),p(2),p(3))
2017
774
      return
2018
775
      end
2019
776
2020
777
      function getinvm4_4(p1,p2,p3,p4)
2021
778
      implicit none
2022
779
      real*8 getinvm4_4,p1(0:3),p2(0:3),p3(0:3),p4(0:3),p(0:3)
2023
780
      real*8 getinvm
2024
781
      integer i
2025
782
      do i=0,3
2026
783
         p(i)=p1(i)+p2(i)+p3(i)+p4(i)
2027
784
      enddo
2028
785
      getinvm4_4=getinvm(p(0),p(1),p(2),p(3))
2029
786
      return
2030
787
      end
2031
788
2032
789
2033
790
      function getdelphiv4(p1,p2)
2034
791
      implicit none
2035
792
      real*8 getdelphiv4,p1(0:3),p2(0:3)
2036
793
      real*8 getdelphi
2037
794
c
2038
795
      getdelphiv4=getdelphi(p1(1),p1(2),
2039
796
     #                      p2(1),p2(2))
2040
797
      return
2041
798
      end
2042
799
2043
800
2044
801
      function getcosv4(q1,q2)
2045
802
      implicit none
2046
803
      real*8 getcosv4,q1(0:3),q2(0:3)
2047
804
      real*8 xnorm1,xnorm2,tmp
2048
805
c
2049
806
      if(q1(0).lt.0.d0.or.q2(0).lt.0.d0)then
2050
807
        getcosv4=-1.d10
2051
808
        return
2052
809
      endif
2053
810
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
2054
811
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
2055
812
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
2056
813
        tmp=-1.d10
2057
814
      else
2058
815
        tmp=q1(1)*q2(1)+q1(2)*q2(2)+q1(3)*q2(3)
2059
816
        tmp=tmp/(xnorm1*xnorm2)
2060
817
        if(abs(tmp).gt.1.d0.and.abs(tmp).le.1.001d0)then
2061
818
          tmp=sign(1.d0,tmp)
2062
819
        elseif(abs(tmp).gt.1.001d0)then
2063
820
          write(*,*)'Error in getcosv4',tmp
2064
821
          stop
2065
822
        endif
2066
823
      endif
2067
824
      getcosv4=tmp
2068
825
      return
2069
826
      end
2070
827
2071
828
2072
829
2073
830
      function getmod(p)
2074
831
      implicit none
2075
832
      double precision p(0:3),getmod
2076
833
2077
834
      getmod=sqrt(p(1)**2+p(2)**2+p(3)**2)
2078
835
2079
836
      return
2080
837
      end
2081
838
2082
839
2083
840
2084
841
      subroutine getperpenv4(q1,q2,qperp)
2085
842
c Normal to the plane defined by \vec{q1},\vec{q2}
2086
843
      implicit none
2087
844
      real*8 q1(0:3),q2(0:3),qperp(0:3)
2088
845
      real*8 xnorm1,xnorm2
2089
846
      integer i
2090
847
c
2091
848
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
2092
849
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
2093
850
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
2094
851
        do i=1,4
2095
852
          qperp(i)=-1.d10
2096
853
        enddo
2097
854
      else
2098
855
        qperp(1)=q1(2)*q2(3)-q1(3)*q2(2)
2099
856
        qperp(2)=q1(3)*q2(1)-q1(1)*q2(3)
2100
857
        qperp(3)=q1(1)*q2(2)-q1(2)*q2(1)
2101
858
        do i=1,3
2102
859
          qperp(i)=qperp(i)/(xnorm1*xnorm2)
2103
860
        enddo
2104
861
        qperp(0)=1.d0
2105
862
      endif
2106
863
      return
2107
864
      end
2108
865
2109
866
2110
867
2111
868
2112
869
      subroutine getwedge(p1,p2,pout)
2113
870
      implicit none
2114
871
      real*8 p1(0:3),p2(0:3),pout(0:3)
2115
872
2116
873
      pout(1)=p1(2)*p2(3)-p1(3)*p2(2)
2117
874
      pout(2)=p1(3)*p2(1)-p1(1)*p2(3)
2118
875
      pout(3)=p1(1)*p2(2)-p1(2)*p2(1)
2119
876
      pout(0)=0d0
2120
877
2121
878
      return
2122
879
      end
2123
880
2124
0
881
2125
=== modified file 'Template/NLO/MCatNLO/Makefile_MadFKS'
2126
--- Template/NLO/MCatNLO/Makefile_MadFKS	2016-02-24 00:19:50 +0000
2127
+++ Template/NLO/MCatNLO/Makefile_MadFKS	2018-12-20 08:18:55 +0000
2128
@@ -1,7 +1,7 @@
2129
1
-include ../../Source/make_opts
1
-include ../../Source/make_opts
2130
2
2
2131
3
DEBUG=
3
DEBUG=
2133
4
FF=$(FC) $(FFLAGS) $(DEBUG)
4
FF=$(FC) $(FFLAGS) $(DEBUG) -std=legacy
2134
5
CPP=$(CXX) $(CXXFLAGS) $(DEBUG) $(INCLOPTION)
5
CPP=$(CXX) $(CXXFLAGS) $(DEBUG) $(INCLOPTION)
2135
6
CC=$(CXX) $(CFLAGS) $(DEBUG) $(INCLOPTION)
6
CC=$(CXX) $(CFLAGS) $(DEBUG) $(INCLOPTION)
2136
7
7
2137
8
8
2138
=== modified file 'Template/NLO/SubProcesses/MC_integer.f'
2139
--- Template/NLO/SubProcesses/MC_integer.f	2015-02-09 18:56:09 +0000
2140
+++ Template/NLO/SubProcesses/MC_integer.f	2018-12-20 08:18:55 +0000
2141
@@ -68,6 +68,7 @@
2142
68
      common/to_readgrid/flat_grid                !Tells if grid read from file
68
      common/to_readgrid/flat_grid                !Tells if grid read from file
2143
69
      if (this_dim.lt.1.or.this_dim.gt.maxdim) then
69
      if (this_dim.lt.1.or.this_dim.gt.maxdim) then
2144
70
         write (*,*) 'Increase maxdim in MC_integer.f',maxdim,this_dim
70
         write (*,*) 'Increase maxdim in MC_integer.f',maxdim,this_dim
2145
71
         stop 1
2146
71
      endif
72
      endif
2147
72
c Set the number of intervales for all the dimensions to zero the very
73
c Set the number of intervales for all the dimensions to zero the very
2148
73
c first time this subroutine is called
74
c first time this subroutine is called
2149
@@ -283,12 +284,8 @@
2150
283
      enddo
284
      enddo
2151
284
c
285
c
2152
285
c Reset the accumulated results because we start new iteration.
286
c Reset the accumulated results because we start new iteration.
2159
286
      do this_dim=1,maxdim
287
      call empty_MC_integer
2160
287
         do i=0,nintervals(this_dim)
288
c
2155
288
            acc(i,this_dim)=0d0
2156
289
            ncall(i,this_dim)=0
2157
290
         enddo
2158
291
      enddo
2161
292
      return
289
      return
2162
293
 999  write (*,*) 'Cannot open "grid.MC_integer" file'
290
 999  write (*,*) 'Cannot open "grid.MC_integer" file'
2163
294
      stop
291
      stop
2164
295
292
2165
=== modified file 'Template/NLO/SubProcesses/analysis_lhe.f'
2166
--- Template/NLO/SubProcesses/analysis_lhe.f	2017-06-19 16:05:33 +0000
2167
+++ Template/NLO/SubProcesses/analysis_lhe.f	2018-12-20 08:18:55 +0000
2168
@@ -7,9 +7,10 @@
2169
7
      implicit none
7
      implicit none
2170
8
      integer nwgt
8
      integer nwgt
2171
9
      character*(*) weights_info(*)
9
      character*(*) weights_info(*)
2172
10
      logical passed_unwgt
2173
10
      integer nwgts,nevents
11
      integer nwgts,nevents
2174
11
      double precision sum_of_wgts
12
      double precision sum_of_wgts
2176
12
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
13
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
2177
13
      logical lopen
14
      logical lopen
2178
14
      nwgts=nwgt
15
      nwgts=nwgt
2179
15
      sum_of_wgts=0d0
16
      sum_of_wgts=0d0
2180
@@ -23,9 +24,10 @@
2181
23
      subroutine analysis_end(xnorm)
24
      subroutine analysis_end(xnorm)
2182
24
      implicit none
25
      implicit none
2183
25
      double precision xnorm
26
      double precision xnorm
2184
27
      logical passed_unwgt
2185
26
      integer nwgts,nevents
28
      integer nwgts,nevents
2186
27
      double precision sum_of_wgts
29
      double precision sum_of_wgts
2188
28
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
30
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
2189
29
      logical lopen
31
      logical lopen
2190
30
      integer                                   npoints
32
      integer                                   npoints
2191
31
      double precision            cross_section
33
      double precision            cross_section
2192
@@ -55,9 +57,10 @@
2193
55
      double precision p(0:4,nexternal)
57
      double precision p(0:4,nexternal)
2194
56
      double precision wgts(*)
58
      double precision wgts(*)
2195
57
      integer ibody
59
      integer ibody
2196
60
      logical passed_unwgt
2197
58
      integer nwgts,nevents
61
      integer nwgts,nevents
2198
59
      double precision sum_of_wgts
62
      double precision sum_of_wgts
2200
60
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
63
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
2201
61
c
64
c
2202
62
      integer i,j,npart
65
      integer i,j,npart
2203
63
      double precision ran2
66
      double precision ran2
2204
@@ -95,23 +98,33 @@
2205
95
      integer                                   npoints
98
      integer                                   npoints
2206
96
      double precision            cross_section
99
      double precision            cross_section
2207
97
      common /for_FixedOrder_lhe/ cross_section,npoints
100
      common /for_FixedOrder_lhe/ cross_section,npoints
2208
101
      double precision,allocatable :: wwgts(:)
2209
98
c --- do the partial unweighting (don't do it when to close to singular
102
c --- do the partial unweighting (don't do it when to close to singular
2210
99
c --- region)
103
c --- region)
2211
104
      if (.not. allocated(wwgts)) then
2212
105
         allocate(wwgts(nwgts))
2213
106
      endif
2214
107
      passed_unwgt=.true.
2215
100
      twgt=abs(cross_section)*FO_LHE_weight_ratio !/npoints
108
      twgt=abs(cross_section)*FO_LHE_weight_ratio !/npoints
2216
101
      if(abs(wgts(1)).lt.abs(twgt) .and. xi_i_fks_ev.gt.1d-3 .and.
109
      if(abs(wgts(1)).lt.abs(twgt) .and. xi_i_fks_ev.gt.1d-3 .and.
2217
102
     $     1d0-y_ij_fks_ev.gt.1d-2)then
110
     $     1d0-y_ij_fks_ev.gt.1d-2)then
2218
103
         R = ran2()*abs(twgt)
111
         R = ran2()*abs(twgt)
2219
104
         if (R.gt.abs(wgts(1)))then
112
         if (R.gt.abs(wgts(1)))then
2220
113
            passed_unwgt=.false.
2221
105
            return
114
            return
2222
106
         else
115
         else
2223
107
            do i=2,nwgts
116
            do i=2,nwgts
2225
108
               wgts(i) = wgts(i)*abs(twgt/wgts(1))
117
               wwgts(i) = wgts(i)*abs(twgt/wgts(1))
2226
109
            enddo
118
            enddo
2228
110
            wgts(1) = sign(twgt,wgts(1))
119
            wwgts(1) = sign(twgt,wgts(1))
2229
111
         endif
120
         endif
2230
121
      else
2231
122
         do i=1,nwgts
2232
123
            wwgts(i)=wgts(i)
2233
124
         enddo
2234
112
      endif
125
      endif
2235
113
c --- accumulate the total weights of all the events in the event file
126
c --- accumulate the total weights of all the events in the event file
2237
114
      sum_of_wgts=sum_of_wgts+wgts(1)
127
      sum_of_wgts=sum_of_wgts+wwgts(1)
2238
115
128
2239
116
c --- fill the multi-weight common blocks with the scale and PDF
129
c --- fill the multi-weight common blocks with the scale and PDF
2240
117
c --- variation weights
130
c --- variation weights
2241
@@ -122,12 +135,12 @@
2242
122
               do ii=1,nint(scalevarF(0))
135
               do ii=1,nint(scalevarF(0))
2243
123
                  do jj=1,nint(scalevarR(0))
136
                  do jj=1,nint(scalevarR(0))
2244
124
                     i_wgt=i_wgt+1
137
                     i_wgt=i_wgt+1
2246
125
                     wgtxsecmu(jj,ii,kk)= wgts(i_wgt)
138
                     wgtxsecmu(jj,ii,kk)= wwgts(i_wgt)
2247
126
                  enddo
139
                  enddo
2248
127
               enddo
140
               enddo
2249
128
            else
141
            else
2250
129
               i_wgt=i_wgt+1
142
               i_wgt=i_wgt+1
2252
130
               wgtxsecmu(1,1,kk)= wgts(i_wgt)
143
               wgtxsecmu(1,1,kk)= wwgts(i_wgt)
2253
131
            endif
144
            endif
2254
132
         enddo
145
         enddo
2255
133
      endif
146
      endif
2256
@@ -136,11 +149,11 @@
2257
136
            if (lpdfvar(nn)) then
149
            if (lpdfvar(nn)) then
2258
137
               do n=0,nmemPDF(nn)
150
               do n=0,nmemPDF(nn)
2259
138
                  i_wgt=i_wgt+1
151
                  i_wgt=i_wgt+1
2261
139
                  wgtxsecPDF(n,nn) = wgts(i_wgt)
152
                  wgtxsecPDF(n,nn) = wwgts(i_wgt)
2262
140
               enddo
153
               enddo
2263
141
            else
154
            else
2264
142
               i_wgt=i_wgt+1
155
               i_wgt=i_wgt+1
2266
143
               wgtxsecPDF(0,nn) = wgts(i_wgt)
156
               wgtxsecPDF(0,nn) = wwgts(i_wgt)
2267
144
            endif
157
            endif
2268
145
         enddo
158
         enddo
2269
146
      endif
159
      endif
2270
@@ -187,7 +200,7 @@
2271
187
200
2272
188
c --- write the event
201
c --- write the event
2273
189
      call write_lhef_event(41,
202
      call write_lhef_event(41,
2275
190
     &     npart,IDPRUP,wgts(1),0d0,0d0,0d0,
203
     &     npart,IDPRUP,wwgts(1),0d0,0d0,0d0,
2276
191
     &     IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
204
     &     IDUP,ISTUP,MOTHUP,ICOLUP,PUP,VTIMUP,SPINUP,buff)
2277
192
205
2278
193
 201  format(a9,1x,i1,4(1x,i2),2(1x,d14.8),2x,i2,2(1x,i2),5(1x,d14.8))
206
 201  format(a9,1x,i1,4(1x,i2),2(1x,d14.8),2x,i2,2(1x,i2),5(1x,d14.8))
2279
@@ -196,12 +209,15 @@
2280
196
c This we can use for the event grouping!      
209
c This we can use for the event grouping!      
2281
197
      subroutine HwU_add_points
210
      subroutine HwU_add_points
2282
198
      implicit none
211
      implicit none
2283
212
      logical passed_unwgt
2284
199
      integer nwgts,nevents
213
      integer nwgts,nevents
2285
200
      double precision sum_of_wgts
214
      double precision sum_of_wgts
2290
201
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents
215
      common/to_lhe_analysis/sum_of_wgts,nwgts,nevents,passed_unwgt
2291
202
      nevents=nevents+1
216
      if (passed_unwgt) then
2292
203
      write (41,'(a)') '</eventgroup>'
217
         nevents=nevents+1
2293
204
      write (41,'(a)') '<eventgroup>'
218
         write (41,'(a)') '</eventgroup>'
2294
219
         write (41,'(a)') '<eventgroup>'
2295
220
      endif
2296
205
      end
221
      end
2297
206
222
2298
207
c Dummy routines
223
c Dummy routines
2299
208
224
2300
=== modified file 'Template/NLO/SubProcesses/check_poles.f'
2301
--- Template/NLO/SubProcesses/check_poles.f	2017-07-31 08:48:20 +0000
2302
+++ Template/NLO/SubProcesses/check_poles.f	2018-12-20 08:18:55 +0000
2303
@@ -18,7 +18,7 @@
2304
18
      include 'genps.inc'
18
      include 'genps.inc'
2305
19
      include 'nexternal.inc'
19
      include 'nexternal.inc'
2306
20
      include 'nFKSconfigs.inc'
20
      include 'nFKSconfigs.inc'
2308
21
      double precision p(0:3, nexternal), prambo(0:3, nexternal)
21
      double precision p(0:3, nexternal), prambo(0:3,100)
2309
22
      double precision p_born(0:3,nexternal-1)
22
      double precision p_born(0:3,nexternal-1)
2310
23
      common/pborn/p_born
23
      common/pborn/p_born
2311
24
      double precision pswgt
24
      double precision pswgt
2312
@@ -48,7 +48,7 @@
2313
48
      include 'coupl.inc'
48
      include 'coupl.inc'
2314
49
      include 'q_es.inc'
49
      include 'q_es.inc'
2315
50
      integer nsqso,MLResArrayDim
50
      integer nsqso,MLResArrayDim
2317
51
      double precision pmass(nexternal), pmass_rambo(nexternal)
51
      double precision pmass(nexternal), pmass_rambo(100)
2318
52
      integer nfail
52
      integer nfail
2319
53
      logical first_time
53
      logical first_time
2320
54
      data first_time/.TRUE./
54
      data first_time/.TRUE./
2321
55
55
2322
=== modified file 'Template/NLO/SubProcesses/combine_plots_FO.sh'
2323
--- Template/NLO/SubProcesses/combine_plots_FO.sh	2013-05-25 12:52:41 +0000
2324
+++ Template/NLO/SubProcesses/combine_plots_FO.sh	2018-12-20 08:18:55 +0000
2325
@@ -1,6 +1,4 @@
2326
1
#!/bin/bash
1
#!/bin/bash
2327
2
TD=td_mac_intel
2328
3
PSOPEN=open
2329
4
2
2330
5
# $? is the value of last executed command. A call to this function
3
# $? is the value of last executed command. A call to this function
2331
6
# after a failure will cause the program to quit the script
4
# after a failure will cause the program to quit the script
2332
@@ -15,126 +13,74 @@
2333
15
fi
13
fi
2334
16
}
14
}
2335
17
15
2343
18
thisdir=`pwd`
16
function combine {
2344
19
if [ -f $thisdir/MADatNLO.top ]
17
counter=0
2345
20
then
18
i=1
2339
21
rm -f $thisdir/MADatNLO*.top
2340
22
fi
2341
23
make read40
2342
24
EXENAME=$thisdir/read40
2346
25
echo -n "" > dir
19
echo -n "" > dir
2380
26
counterp=0
20
for plot_file in "${@:2}" ; do
2381
27
for p in P* ; do
21
    counter=$(($counter + 1))
2382
28
    cd $thisdir
22
    echo $plot_file >> dir
2383
29
    cd $p
23
    if [[ $(($counter % 40)) == 0 ]]; then
2384
30
    rm -f MADatNLO*
24
        echo -n "" > dir2
2385
31
    rm -f *read40*.out
25
	echo "40" > dir2
2353
32
    echo -n "" > dir
2354
33
    counter=0
2355
34
    i=1
2356
35
    for g in $* ; do
2357
36
        if [ ! -f $g/MADatNLO*.top ]
2358
37
        then
2359
38
          echo "No topdrawer file in $p/$g"
2360
39
	else
2361
40
	counter=$(($counter + 1))
2362
41
	echo $g"/MADatNLO.top" >> dir
2363
42
	if [[ $(($counter % 40)) == 0 ]]; then
2364
43
            echo -n "" > dir2
2365
44
	    echo "40" > dir2
2366
45
	    cat dir >> dir2
2367
46
	    $EXENAME < dir2
2368
47
            teststatus Failure in step 1
2369
48
	    rm -f dir
2370
49
	    rm -f dir2
2371
50
	    mv fort.88 MADatNLO$i\.top 
2372
51
	    mv read40.out "S1read40_"$i\.out
2373
52
	    i=$(($i + 1))
2374
53
	    echo -n "" > dir
2375
54
        fi
2376
55
	fi
2377
56
    done
2378
57
    if [[ $(($counter % 40)) != 0 ]] ; then
2379
58
	echo $(($counter % 40)) > dir2
2386
59
	cat dir >> dir2
26
	cat dir >> dir2
2392
60
	while [[ $(($counter % 40)) -ne 0 ]]; do
27
	./read40 < dir2
2393
61
	    counter=$[$counter-1]
28
        teststatus Failure in step 1
2389
62
	done
2390
63
	$EXENAME < dir2
2391
64
	teststatus Failure in step 2
2394
65
	rm -f dir
29
	rm -f dir
2395
66
	rm -f dir2
30
	rm -f dir2
2446
67
	mv fort.88 MADatNLO$i\.top 
31
	mv fort.88 $1_$i\.top 
2447
68
	mv read40.out "S2read40_"$i\.out
32
	mv read40.out "S1read40_"$1_$i\.out
2448
69
    fi
33
	i=$(($i + 1))
2449
70
	
34
	echo -n "" > dir
2450
71
    counter=0
35
    fi
2451
72
    for m in MADatNLO*.top ; do
36
done
2452
73
	counter=$[$counter+1]
37
if [[ $(($counter % 40)) != 0 ]] ; then
2453
74
	echo $m >> dir
38
    echo $(($counter % 40)) > dir2
2454
75
    done
39
    cat dir >> dir2
2455
76
    if [[ $counter -gt 1 ]] ; then
40
    ./read40 < dir2
2456
77
	echo $counter > dir2
41
    teststatus Failure in step 2
2407
78
	cat dir >> dir2
2408
79
	$EXENAME < dir2
2409
80
	teststatus Failure in step 3
2410
81
	rm -f dir2
2411
82
	mv fort.88 MADatNLO.top 
2412
83
	mv read40.out S3read40.out
2413
84
    else
2414
85
	mv MADatNLO1.top MADatNLO.top
2415
86
    fi
2416
87
    rm -f dir
2417
88
    cd $thisdir
2418
89
    counterp=$(($counterp + 1))
2419
90
    echo $p"/MADatNLO.top" >> dir
2420
91
    echo $counterp $p "done"
2421
92
done
2422
93
2423
94
cd $thisdir
2424
95
make split40
2425
96
./split40
2426
97
2427
98
i=0
2428
99
#pdir* are created by split
2429
100
for p in pdir* ; do
2430
101
    i=$(($i+1))
2431
102
    $EXENAME < $p    
2432
103
    teststatus Failure in step 4
2433
104
    mv fort.88 MADatNLO$i\.top
2434
105
    mv read40.out "S4read40_"$i\.out
2435
106
    rm $p
2436
107
    echo $i $p "done"
2437
108
done
2438
109
2439
110
if [[ $i -gt 1 ]] ; then
2440
111
    echo $i > dir2
2441
112
    for p in MADatNLO*.top ; do
2442
113
	echo $p >>dir2
2443
114
    done
2444
115
    $EXENAME < dir2
2445
116
    teststatus Failure in step 5
2457
117
    rm -f dir
42
    rm -f dir
2458
118
    rm -f dir2
43
    rm -f dir2
2461
119
    mv fort.88 MADatNLO.top 
44
    mv fort.88 $1_$i\.top 
2462
120
    mv read40.out S5read40.out
45
    mv read40.out "S2read40_"$1_$i\.out
2463
46
fi
2464
47
}
2465
48
2466
49
make read40
2467
50
rm -f MADatNLO*.top
2468
51
2469
52
if [ "$#" -gt 64000 ]; then
2470
53
    combine "MADatNLO_A" "$@"
2471
54
    combine "MADatNLO_B" MADatNLO_A*.top
2472
55
    combine "MADatNLO_C" MADatNLO_B*.top
2473
56
    combine "MADatNLO_D" MADatNLO_C*.top
2474
57
    mv MADatNLO_D_1.top MADatNLO.top
2475
58
elif [ "$#" -gt 1600 ]; then
2476
59
    combine "MADatNLO_A" "$@"
2477
60
    combine "MADatNLO_B" MADatNLO_A*.top
2478
61
    combine "MADatNLO_C" MADatNLO_B*.top
2479
62
    mv MADatNLO_C_1.top MADatNLO.top
2480
63
elif [ "$#" -gt 40 ]; then
2481
64
    combine "MADatNLO_A" "$@"
2482
65
    combine "MADatNLO_B" MADatNLO_A*.top
2483
66
    mv MADatNLO_B_1.top MADatNLO.top
2484
67
elif [ "$#" -gt 1 ]; then
2485
68
    combine "MADatNLO_A" "$@"
2486
69
    mv MADatNLO_A_1.top MADatNLO.top
2487
121
else
70
else
2489
122
    mv MADatNLO1.top MADatNLO.top
71
    cp $1 MADatNLO.top
2490
123
fi
72
fi
2491
124
73
2493
125
if [ -f $thisdir/read40.errors ]
74
2494
75
if [ -f read40.errors ]
2495
126
then
76
then
2497
127
rm -f $thisdir/read40.errors
77
rm -f read40.errors
2498
128
fi
78
fi
2502
129
find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > $thisdir/read40.errors
79
2503
130
80
find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > read40.errors
2504
131
if [ ! -s $thisdir/read40.errors ] ; then
81
2505
82
if [ ! -s read40.errors ] ; then
2506
132
    echo 'no errors found'
83
    echo 'no errors found'
2507
133
else
84
else
2508
134
    echo 'ERRORS! see read40.errors for details'
85
    echo 'ERRORS! see read40.errors for details'
2509
135
fi
86
fi
2510
136
2511
137
#./NLO_Born3.py
2512
138
#$TD MADatNLO_combined.top "POSTSCR,ORIENT=3"
2513
139
#$PSOPEN MADatNLO_combined.ps
2514
140
##td2pdf MADatNLO_con
2515
141
87
2516
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
2517
--- Template/NLO/SubProcesses/driver_mintFO.f	2018-04-16 14:08:47 +0000
2518
+++ Template/NLO/SubProcesses/driver_mintFO.f	2018-12-20 08:18:55 +0000
2519
@@ -48,7 +48,7 @@
2520
48
      include "mint.inc"
48
      include "mint.inc"
2521
49
      integer nhits_in_grids(maxchannels)
49
      integer nhits_in_grids(maxchannels)
2522
50
      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
50
      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
2524
51
     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
51
     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
2525
52
     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
52
     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
2526
53
     $     ,0:maxchannels),x(ndimmax),itmax_fl
53
     $     ,0:maxchannels),x(ndimmax),itmax_fl
2527
54
      integer ixi_i,iphi_i,iy_ij,vn
54
      integer ixi_i,iphi_i,iy_ij,vn
2528
@@ -402,9 +402,10 @@
2529
402
      double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
402
      double precision xx(ndimmax),vegas_wgt,f(nintegrals),jac,p(0:3
2530
403
     $     ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
403
     $     ,nexternal),rwgt,vol,sig,x(99),MC_int_wgt
2531
404
      integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min,iamp
404
      integer ifl,nFKS_born,nFKS_picked,iFKS,nFKS_min,iamp
2533
405
     $     ,nFKS_max,izero,ione,itwo,mohdr
405
     $     ,nFKS_max,izero,ione,itwo,mohdr,i,iran_picked
2534
406
      parameter (izero=0,ione=1,itwo=2,mohdr=-100)
406
      parameter (izero=0,ione=1,itwo=2,mohdr=-100)
2536
407
      logical passcuts,passcuts_nbody,passcuts_n1body,sum
407
      logical passcuts,passcuts_nbody,passcuts_n1body,sum,firsttime
2537
408
      data firsttime/.true./
2538
408
      external passcuts
409
      external passcuts
2539
409
      integer             ini_fin_fks(maxchannels)
410
      integer             ini_fin_fks(maxchannels)
2540
410
      common/fks_channels/ini_fin_fks
411
      common/fks_channels/ini_fin_fks
2541
@@ -431,6 +432,16 @@
2542
431
      common /for_applgrid/ iappl
432
      common /for_applgrid/ iappl
2543
432
      double precision       wgt_ME_born,wgt_ME_real
433
      double precision       wgt_ME_born,wgt_ME_real
2544
433
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
434
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
2545
435
      integer ini_fin_fks_map(0:2,0:fks_configs)
2546
436
      save ini_fin_fks_map
2547
437
      if (firsttime) then
2548
438
         firsttime=.false.
2549
439
         call setup_ini_fin_fks_map(ini_fin_fks_map)
2550
440
         write (*,*) 'initial-final FKS maps:'
2551
441
         write (*,*) 0 ,':',ini_fin_fks_map(0,:)
2552
442
         write (*,*) 1 ,':',ini_fin_fks_map(1,:)
2553
443
         write (*,*) 2 ,':',ini_fin_fks_map(2,:)
2554
444
      endif
2555
434
      if (ifl.ne.0) then
445
      if (ifl.ne.0) then
2556
435
         write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
446
         write (*,*) 'ERROR ifl not equal to zero in sigint',ifl
2557
436
         stop 1
447
         stop 1
2558
@@ -453,20 +464,9 @@
2559
453
      if (ickkw.eq.-1) H1_factor_virt=0d0
464
      if (ickkw.eq.-1) H1_factor_virt=0d0
2560
454
      if (ickkw.eq.3) call set_FxFx_scale(0,p)
465
      if (ickkw.eq.3) call set_FxFx_scale(0,p)
2561
455
      call update_vegas_x(xx,x)
466
      call update_vegas_x(xx,x)
2576
456
      call get_MC_integer(max(ini_fin_fks(ichan),1),fks_configs
467
      call get_MC_integer(max(ini_fin_fks(ichan),1)
2577
457
     $     ,nFKS_picked,vol)
468
     $     ,ini_fin_fks_map(ini_fin_fks(ichan),0),iran_picked,vol)
2578
458
      if (ini_fin_fks(ichan).eq.1) then
469
      nFKS_picked=ini_fin_fks_map(ini_fin_fks(ichan),iran_picked)
2565
459
         do while (fks_j_d(nFKS_picked).le.nincoming) 
2566
460
            call get_MC_integer(ini_fin_fks(ichan),fks_configs
2567
461
     $           ,nFKS_picked,vol)
2568
462
         enddo
2569
463
      elseif (ini_fin_fks(ichan).eq.2) then
2570
464
         do while (fks_j_d(nFKS_picked).gt.nincoming) 
2571
465
            call get_MC_integer(ini_fin_fks(ichan),fks_configs
2572
466
     $           ,nFKS_picked,vol)
2573
467
         enddo
2574
468
      endif
2575
469
2579
470
      
470
      
2580
471
c The nbody contributions
471
c The nbody contributions
2581
472
      if (abrv.eq.'real') goto 11
472
      if (abrv.eq.'real') goto 11
2582
@@ -503,14 +503,15 @@
2583
503
      nbody=.false.
503
      nbody=.false.
2584
504
      if (sum) then
504
      if (sum) then
2585
505
         nFKS_min=1
505
         nFKS_min=1
2587
506
         nFKS_max=fks_configs
506
         nFKS_max=ini_fin_fks_map(ini_fin_fks(ichan),0)
2588
507
         MC_int_wgt=1d0
507
         MC_int_wgt=1d0
2589
508
      else
508
      else
2592
509
         nFKS_min=nFKS_picked
509
         nFKS_min=iran_picked
2593
510
         nFKS_max=nFKS_picked
510
         nFKS_max=iran_picked
2594
511
         MC_int_wgt=1d0/vol
511
         MC_int_wgt=1d0/vol
2595
512
      endif
512
      endif
2597
513
      do iFKS=nFKS_min,nFKS_max
513
      do i=nFKS_min,nFKS_max
2598
514
         iFKS=ini_fin_fks_map(ini_fin_fks(ichan),i)
2599
514
         calculatedBorn=.false. 
515
         calculatedBorn=.false. 
2600
515
         ! MZ this is a temporary fix for processes without
516
         ! MZ this is a temporary fix for processes without
2601
516
         ! soft singularities associated to the initial state
517
         ! soft singularities associated to the initial state
2602
@@ -566,11 +567,11 @@
2603
566
c Importance sampling for FKS configurations
567
c Importance sampling for FKS configurations
2604
567
      if (sum) then
568
      if (sum) then
2605
568
         call get_wgt_nbody(sig)
569
         call get_wgt_nbody(sig)
2607
569
         call fill_MC_integer(max(ini_fin_fks(ichan),1),nFKS_picked
570
         call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
2608
570
     $        ,abs(sig))
571
     $        ,abs(sig))
2609
571
      else
572
      else
2610
572
         call get_wgt_no_nbody(sig)
573
         call get_wgt_no_nbody(sig)
2612
573
         call fill_MC_integer(max(ini_fin_fks(ichan),1),nFKS_picked
574
         call fill_MC_integer(max(ini_fin_fks(ichan),1),iran_picked
2613
574
     $        ,abs(sig)*vol)
575
     $        ,abs(sig)*vol)
2614
575
      endif
576
      endif
2615
576
577
2616
@@ -595,6 +596,35 @@
2617
595
      return
596
      return
2618
596
      end
597
      end
2619
597
      
598
      
2620
599
      subroutine setup_ini_fin_FKS_map(ini_fin_FKS_map)
2621
600
      implicit none
2622
601
      include 'nexternal.inc'
2623
602
      include 'nFKSconfigs.inc'
2624
603
      include 'fks_info.inc'
2625
604
      integer ini_fin_FKS_map(0:2,0:fks_configs),iFKS
2626
605
      ini_fin_FKS_map(0,0)=0
2627
606
      ini_fin_FKS_map(1,0)=0
2628
607
      ini_fin_FKS_map(2,0)=0
2629
608
      do iFKS=1,fks_configs
2630
609
         ini_fin_FKS_map(0,0)=ini_fin_FKS_map(0,0)+1
2631
610
         ini_fin_FKS_map(0,ini_fin_FKS_map(0,0))=iFKS
2632
611
         if (fks_j_d(iFKS).le.nincoming .and.
2633
612
     $       fks_j_d(iFKS).gt.0) then
2634
613
            ini_fin_FKS_map(2,0)=ini_fin_FKS_map(2,0)+1
2635
614
            ini_fin_FKS_map(2,ini_fin_FKS_map(2,0))=iFKS
2636
615
         elseif (fks_j_d(iFKS).gt.nincoming .and.
2637
616
     $           fks_j_d(iFKS).le.nexternal) then
2638
617
            ini_fin_FKS_map(1,0)=ini_fin_FKS_map(1,0)+1
2639
618
            ini_fin_FKS_map(1,ini_fin_FKS_map(1,0))=iFKS
2640
619
         else
2641
620
            write (*,*) 'ERROR in setup_ini_fin_FKS_map',fks_j_d(iFKS)
2642
621
     $           ,nincoming,iFKS
2643
622
            stop 1
2644
623
         endif
2645
624
      enddo
2646
625
      return
2647
626
      end
2648
627
2649
598
      subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)
628
      subroutine get_born_nFKSprocess(nFKS_in,nFKS_out)
2650
599
      implicit none
629
      implicit none
2651
600
      include 'nexternal.inc'
630
      include 'nexternal.inc'
2652
601
631
2653
=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
2654
--- Template/NLO/SubProcesses/driver_mintMC.f	2018-04-16 14:08:47 +0000
2655
+++ Template/NLO/SubProcesses/driver_mintMC.f	2018-12-20 08:18:55 +0000
2656
@@ -59,7 +59,7 @@
2657
59
      common /event_normalisation/event_norm
59
      common /event_normalisation/event_norm
2658
60
c For MINT:
60
c For MINT:
2659
61
      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
61
      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
2661
62
     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
62
     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
2662
63
     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
63
     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
2663
64
     $     ,0:maxchannels),x(ndimmax)
64
     $     ,0:maxchannels),x(ndimmax)
2664
65
      integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)
65
      integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)
2665
@@ -842,7 +842,7 @@
2666
842
         firsttime=.false.
842
         firsttime=.false.
2667
843
c Determines the proc_map that sets which FKS configuration can be
843
c Determines the proc_map that sets which FKS configuration can be
2668
844
c summed explicitly and which by MC-ing.
844
c summed explicitly and which by MC-ing.
2670
845
         call setup_proc_map(sum,proc_map)
845
         call setup_proc_map(sum,proc_map,ini_fin_fks)
2671
846
c For the S-events, we can combine processes when they give identical
846
c For the S-events, we can combine processes when they give identical
2672
847
c processes at the Born. Make sure we check that we get indeed identical
847
c processes at the Born. Make sure we check that we get indeed identical
2673
848
c IRPOC's
848
c IRPOC's
2674
@@ -869,16 +869,6 @@
2675
869
         if (ickkw.eq.3) call set_FxFx_scale(0,p)
869
         if (ickkw.eq.3) call set_FxFx_scale(0,p)
2676
870
         call update_vegas_x(xx,x)
870
         call update_vegas_x(xx,x)
2677
871
         call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
871
         call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
2678
872
         if (ini_fin_fks.eq.1) then
2679
873
            do while (fks_j_d(proc_map(proc_map(0,1),1)).le.nincoming) 
2680
874
               call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
2681
875
            enddo
2682
876
         elseif (ini_fin_fks.eq.2) then
2683
877
            do while (fks_j_d(proc_map(proc_map(0,1),1)).gt.nincoming) 
2684
878
               call get_MC_integer(1,proc_map(0,0),proc_map(0,1),vol1)
2685
879
            enddo
2686
880
         endif
2687
881
2688
882
c The nbody contributions
872
c The nbody contributions
2689
883
         if (abrv.eq.'real') goto 11
873
         if (abrv.eq.'real') goto 11
2690
884
         nbody=.true.
874
         nbody=.true.
2691
@@ -1043,7 +1033,7 @@
2692
1043
      end
1033
      end
2693
1044
1034
2694
1045
1035
2696
1046
      subroutine setup_proc_map(sum,proc_map)
1036
      subroutine setup_proc_map(sum,proc_map,ini_fin_fks)
2697
1047
c Determines the proc_map that sets which FKS configuration can be
1037
c Determines the proc_map that sets which FKS configuration can be
2698
1048
c summed explicitly and which by MC-ing.
1038
c summed explicitly and which by MC-ing.
2699
1049
      implicit none
1039
      implicit none
2700
@@ -1056,7 +1046,7 @@
2701
1056
      logical found_ini1,found_ini2,found_fnl
1046
      logical found_ini1,found_ini2,found_fnl
2702
1057
      integer proc_map(0:fks_configs,0:fks_configs)
1047
      integer proc_map(0:fks_configs,0:fks_configs)
2703
1058
     $     ,j_fks_proc(fks_configs),i_fks_pdg_proc(fks_configs)
1048
     $     ,j_fks_proc(fks_configs),i_fks_pdg_proc(fks_configs)
2705
1059
     $     ,j_fks_pdg_proc(fks_configs),i,sum,j
1049
     $     ,j_fks_pdg_proc(fks_configs),i,sum,j,ini_fin_fks
2706
1060
      integer              nFKSprocess
1050
      integer              nFKSprocess
2707
1061
      common/c_nFKSprocess/nFKSprocess
1051
      common/c_nFKSprocess/nFKSprocess
2708
1062
      INTEGER              IPROC
1052
      INTEGER              IPROC
2709
@@ -1097,6 +1087,8 @@
2710
1097
c them in the process map
1087
c them in the process map
2711
1098
         do nFKSprocess=1,fks_configs
1088
         do nFKSprocess=1,fks_configs
2712
1099
            call fks_inc_chooser()
1089
            call fks_inc_chooser()
2713
1090
            if (ini_fin_fks.eq.1 .and. j_fks.le.nincoming) cycle
2714
1091
            if (ini_fin_fks.eq.2 .and. j_fks.gt.nincoming) cycle
2715
1100
            if (need_color_links.or.need_charge_links) then
1092
            if (need_color_links.or.need_charge_links) then
2716
1101
               proc_map(0,0)=proc_map(0,0)+1
1093
               proc_map(0,0)=proc_map(0,0)+1
2717
1102
               proc_map(proc_map(0,0),0)=proc_map(proc_map(0,0),0)+1
1094
               proc_map(proc_map(0,0),0)=proc_map(proc_map(0,0),0)+1
2718
@@ -1156,6 +1148,8 @@
2719
1156
c gluons splitting
1148
c gluons splitting
2720
1157
         do nFKSprocess=1,fks_configs
1149
         do nFKSprocess=1,fks_configs
2721
1158
            call fks_inc_chooser()
1150
            call fks_inc_chooser()
2722
1151
            if (ini_fin_fks.eq.1 .and. j_fks.le.nincoming) cycle
2723
1152
            if (ini_fin_fks.eq.2 .and. j_fks.gt.nincoming) cycle
2724
1159
            if (.not.(need_color_links.or.need_charge_links)) then
1153
            if (.not.(need_color_links.or.need_charge_links)) then
2725
1160
               if (j_fks.eq.1 .and. found_ini1) then
1154
               if (j_fks.eq.1 .and. found_ini1) then
2726
1161
                  do i=1,proc_map(0,0)
1155
                  do i=1,proc_map(0,0)
2727
1162
1156
2728
=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
2729
--- Template/NLO/SubProcesses/fks_singular.f	2018-04-16 14:08:47 +0000
2730
+++ Template/NLO/SubProcesses/fks_singular.f	2018-12-20 08:18:55 +0000
2731
@@ -1115,8 +1115,6 @@
2732
1115
      double precision    xi_i_fks_ev,y_ij_fks_ev
1115
      double precision    xi_i_fks_ev,y_ij_fks_ev
2733
1116
      double precision    p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
1116
      double precision    p_i_fks_ev(0:3),p_i_fks_cnt(0:3,-2:2)
2734
1117
      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
1117
      common/fksvariables/xi_i_fks_ev,y_ij_fks_ev,p_i_fks_ev,p_i_fks_cnt
2735
1118
      double precision  xi_i_hat_ev,xi_i_hat_cnt(-2:2)
2736
1119
      common /cxi_i_hat/xi_i_hat_ev,xi_i_hat_cnt
2737
1120
      integer            i_fks,j_fks
1118
      integer            i_fks,j_fks
2738
1121
      common/fks_indices/i_fks,j_fks
1119
      common/fks_indices/i_fks,j_fks
2739
1122
      double precision   xi_i_fks_cnt(-2:2)
1120
      double precision   xi_i_fks_cnt(-2:2)
2740
@@ -1249,7 +1247,7 @@
2741
1249
      endif
1247
      endif
2742
1250
1248
2743
1251
c f_* multiplication factors for real-emission, soft counter, ... etc.       
1249
c f_* multiplication factors for real-emission, soft counter, ... etc.       
2745
1252
      prefact=1d0/xi_i_hat_ev/(1-y_ij_fks_ev)
1250
      prefact=xinorm_ev/xi_i_fks_ev/(1-y_ij_fks_ev)
2746
1253
      f_r=prefact*jac_ev*enhance_real*fkssymmetryfactor*vegas_wgt
1251
      f_r=prefact*jac_ev*enhance_real*fkssymmetryfactor*vegas_wgt
2747
1254
      f_MC_S=f_r
1252
      f_MC_S=f_r
2748
1255
      f_MC_H=f_r
1253
      f_MC_H=f_r
2749
@@ -1736,11 +1734,16 @@
2750
1736
c coefficients for PDF and scale computations. 
1734
c coefficients for PDF and scale computations. 
2751
1737
      use weight_lines
1735
      use weight_lines
2752
1738
      implicit none
1736
      implicit none
2754
1739
      integer i,j
1737
      include 'orders.inc'
2755
1738
      include 'mint.inc'
2756
1739
      integer orders(nsplitorders)
2757
1740
      integer i,j,iamp
2758
1741
      logical virt_found
2759
1740
      double precision bias
1742
      double precision bias
2760
1741
      character*7 event_norm
1743
      character*7 event_norm
2761
1742
      common /event_normalisation/event_norm
1744
      common /event_normalisation/event_norm
2763
1743
      double precision           virt_wgt_mint,born_wgt_mint
1745
      double precision           virt_wgt_mint(0:n_ave_virt),
2764
1746
     &                           born_wgt_mint(0:n_ave_virt)
2765
1744
      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
1747
      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
2766
1745
c Set the bias_wgt to 1 in case we do not have to do any biassing
1748
c Set the bias_wgt to 1 in case we do not have to do any biassing
2767
1746
      if (event_norm(1:4).ne.'bias') then
1749
      if (event_norm(1:4).ne.'bias') then
2768
@@ -1749,6 +1752,7 @@
2769
1749
         enddo
1752
         enddo
2770
1750
         return
1753
         return
2771
1751
      endif
1754
      endif
2772
1755
      virt_found=.false.
2773
1752
c loop over all contributions
1756
c loop over all contributions
2774
1753
      do i=1,icontr
1757
      do i=1,icontr
2775
1754
         if (itype(i).eq.1) then
1758
         if (itype(i).eq.1) then
2776
@@ -1771,9 +1775,15 @@
2777
1771
         do j=1,3
1775
         do j=1,3
2778
1772
            wgt(j,i)=wgt(j,i)*bias_wgt(i)
1776
            wgt(j,i)=wgt(j,i)*bias_wgt(i)
2779
1773
         enddo
1777
         enddo
2783
1774
         if (itype(i).eq.14) then
1778
         if (itype(i).eq.14 .and. .not.virt_found) then
2784
1775
            virt_wgt_mint=virt_wgt_mint*bias_wgt(i)
1779
            virt_found=.true.
2785
1776
            born_wgt_mint=born_wgt_mint*bias_wgt(i)
1780
            virt_wgt_mint(0)=virt_wgt_mint(0)*bias_wgt(i)
2786
1781
            born_wgt_mint(0)=born_wgt_mint(0)*bias_wgt(i)
2787
1782
            do iamp=1,amp_split_size
2788
1783
               call amp_split_pos_to_orders(iamp, orders)
2789
1784
               virt_wgt_mint(iamp)=virt_wgt_mint(iamp)*bias_wgt(i)
2790
1785
               born_wgt_mint(iamp)=born_wgt_mint(iamp)*bias_wgt(i)
2791
1786
            enddo
2792
1777
         endif
1787
         endif
2793
1778
      enddo
1788
      enddo
2794
1779
      return
1789
      return
2795
@@ -3805,7 +3815,7 @@
2796
3805
      if (iextra_cnt.gt.0)
3815
      if (iextra_cnt.gt.0)
2797
3806
     1    call extra_cnt(p_born, iextra_cnt, ans_extra_cnt)
3816
     1    call extra_cnt(p_born, iextra_cnt, ans_extra_cnt)
2798
3807
      call AP_reduced(j_type,i_type,ch_j,ch_i,t,z,ap)
3817
      call AP_reduced(j_type,i_type,ch_j,ch_i,t,z,ap)
2800
3808
      call Qterms_reduced_timelike(j_type,i_type,ch_i,ch_j,t,z,Q)
3818
      call Qterms_reduced_timelike(j_type,i_type,ch_j,ch_i,t,z,Q)
2801
3809
      wgt=0d0
3819
      wgt=0d0
2802
3810
      do iord = 1, nsplitorders
3820
      do iord = 1, nsplitorders
2803
3811
         if (.not.split_type(iord) .or. (iord.ne.qed_pos .and.
3821
         if (.not.split_type(iord) .or. (iord.ne.qed_pos .and.
2804
3812
3822
2805
=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
2806
--- Template/NLO/SubProcesses/genps_fks.f	2018-01-29 05:41:33 +0000
2807
+++ Template/NLO/SubProcesses/genps_fks.f	2018-12-20 08:18:55 +0000
2808
@@ -2557,7 +2557,7 @@
2809
2557
      integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
2557
      integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
2810
2558
      double precision cBW_mass(-1:1,-nexternal:-1),
2558
      double precision cBW_mass(-1:1,-nexternal:-1),
2811
2559
     &     cBW_width(-1:1,-nexternal:-1)
2559
     &     cBW_width(-1:1,-nexternal:-1)
2813
2560
      double precision s_mass(-nexternal:-1)
2560
      double precision s_mass(-nexternal:nexternal)
2814
2561
      common/to_phase_space_s_channel/s_mass
2561
      common/to_phase_space_s_channel/s_mass
2815
2562
c Common block with granny information
2562
c Common block with granny information
2816
2563
      logical granny_is_res
2563
      logical granny_is_res
2817
2564
2564
2818
=== modified file 'Template/NLO/SubProcesses/makefile'
2819
--- Template/NLO/SubProcesses/makefile	2017-06-15 16:22:23 +0000
2820
+++ Template/NLO/SubProcesses/makefile	2018-12-20 08:18:55 +0000
2821
@@ -7,9 +7,6 @@
2822
7
# Files for the read40 combiner of top drawer files
7
# Files for the read40 combiner of top drawer files
2823
8
READ40=read40.o
8
READ40=read40.o
2824
9
9
2825
10
# Files for the split40 splitter of top drawer files
2826
11
SPLIT40=split40.o
2827
12
2828
13
# Files to collect all the events in the separate integration channels into a single file
10
# Files to collect all the events in the separate integration channels into a single file
2829
14
COLLECT_EVENTS=collect_events.o handling_lhe_events.o fill_MC_mshell.o
11
COLLECT_EVENTS=collect_events.o handling_lhe_events.o fill_MC_mshell.o
2830
15
12
2831
@@ -21,9 +18,6 @@
2832
21
read40: $(READ40)
18
read40: $(READ40)
2833
22
	$(FC) $(LDFLAGS) -o read40 $(READ40)
19
	$(FC) $(LDFLAGS) -o read40 $(READ40)
2834
23
20
2835
24
split40: $(SPLIT40)
2836
25
	$(FC) $(LDFLAGS) -o split40 $(SPLIT40)
2837
26
2838
27
collect_events: $(COLLECT_EVENTS) $(LIBS)
21
collect_events: $(COLLECT_EVENTS) $(LIBS)
2839
28
	$(FC) -o collect_events $(COLLECT_EVENTS) $(LINKLIBS) $(LDFLAGS)
22
	$(FC) -o collect_events $(COLLECT_EVENTS) $(LINKLIBS) $(LDFLAGS)
2840
29
	rm handling_lhe_events.o
23
	rm handling_lhe_events.o
2841
30
24
2842
=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
2843
--- Template/NLO/SubProcesses/mint-integrator2.f	2017-07-31 08:48:20 +0000
2844
+++ Template/NLO/SubProcesses/mint-integrator2.f	2018-12-20 08:18:55 +0000
2845
@@ -77,7 +77,7 @@
2846
77
     $     ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
77
     $     ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
2847
78
     $     ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
78
     $     ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
2848
79
     $     ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan,error_virt
79
     $     ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan,error_virt
2850
80
      real * 8 xint_virt(maxchannels),ymax_virt(maxchannels)
80
      real * 8 xint_virt(maxchannels),ymax_virt(0:maxchannels)
2851
81
     $     ,ans_chan(0:maxchannels)
81
     $     ,ans_chan(0:maxchannels)
2852
82
      real * 8 x(ndimmax),vol
82
      real * 8 x(ndimmax),vol
2853
83
      real * 8 xacc(0:nintervals,ndimmax,maxchannels)
83
      real * 8 xacc(0:nintervals,ndimmax,maxchannels)
2854
@@ -216,7 +216,7 @@
2855
216
            ans_chan(0)=ans_chan(0)+ans(1,kchan)
216
            ans_chan(0)=ans_chan(0)+ans(1,kchan)
2856
217
         enddo
217
         enddo
2857
218
      endif
218
      endif
2859
219
      cross_section=ans_chan(0)
219
      cross_section=ans_chan(0) * wgt_mult
2860
220
      nit=0
220
      nit=0
2861
221
      nit_included=0
221
      nit_included=0
2862
222
      do i=1,nintegrals
222
      do i=1,nintegrals
2863
@@ -1031,8 +1031,8 @@
2864
1031
      integer ndim,imode
1031
      integer ndim,imode
2865
1032
      include "mint.inc"
1032
      include "mint.inc"
2866
1033
      real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
1033
      real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
2869
1034
     $     ,ymax(nintervals,ndimmax,maxchannels),ymax_virt(maxchannels)
1034
     $     ,ymax(nintervals,ndimmax,maxchannels)
2870
1035
     $     ,x(ndimmax)
1035
     $     ,ymax_virt(0:maxchannels),x(ndimmax)
2871
1036
      real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
1036
      real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
2872
1037
      integer icell(ndimmax),ncell(ndimmax),ncell_virt
1037
      integer icell(ndimmax),ncell(ndimmax),ncell_virt
2873
1038
      integer ifold(ndimmax),kfold(ndimmax)
1038
      integer ifold(ndimmax),kfold(ndimmax)
2874
1039
1039
2875
=== modified file 'Template/NLO/SubProcesses/montecarlocounter.f'
2876
--- Template/NLO/SubProcesses/montecarlocounter.f	2018-04-16 14:08:47 +0000
2877
+++ Template/NLO/SubProcesses/montecarlocounter.f	2018-12-20 08:18:55 +0000
2878
@@ -1,3231 +1,3259 @@
2879
1
      subroutine set_mc_matrices
2880
2
      implicit none
2881
3
      include "genps.inc"
2882
4
      include 'nexternal.inc'
2883
5
      include "born_nhel.inc"
2884
6
c Nexternal is the number of legs (initial and final) al NLO, while max_bcol
2885
7
c is the number of color flows at Born level
2886
8
      integer i,j,k,l,k0,mothercol(2),i1(2)
2887
9
      integer idup(nexternal-1,maxproc)
2888
10
      integer mothup(2,nexternal-1,maxproc)
2889
11
      integer icolup(2,nexternal-1,max_bcol)
2890
12
      include 'born_leshouche.inc'
2891
13
      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
2892
14
      common /MC_info/ ipartners,colorflow
2893
15
      integer i_fks,j_fks
2894
16
      common/fks_indices/i_fks,j_fks
2895
17
      integer fksfather
2896
18
      common/cfksfather/fksfather
2897
19
      logical notagluon,found
2898
20
      common/cnotagluon/notagluon
2899
21
      integer nglu,nsngl
2900
22
      logical isspecial,isspecial0
2901
23
      common/cisspecial/isspecial
2902
24
      logical spec_case
2903
25
2904
26
      include 'orders.inc'
2905
27
      logical split_type(nsplitorders) 
2906
28
      common /c_split_type/split_type
2907
29
      double precision particle_charge(nexternal)
2908
30
      common /c_charges/particle_charge
2909
31
2910
32
c
2911
33
      logical is_leading_cflow(max_bcol)
2912
34
      integer num_leading_cflows
2913
35
      common/c_leading_cflows/is_leading_cflow,num_leading_cflows
2914
36
      include 'born_conf.inc'
2915
37
      include 'born_coloramps.inc'
2916
38
c
2917
39
      ipartners(0)=0
2918
40
      do i=1,nexternal-1
2919
41
         colorflow(i,0)=0
2920
42
      enddo
2921
43
2922
44
C What follows is true for QCD-type splittings.
2923
45
C For QED-type splittings, ipartner is simply all the charged particles
2924
46
C in the event except for FKSfather. In this case, all the born color
2925
47
C flows are allowed
2926
48
2927
49
c ipartners(0): number of particles that can be colour or anticolour partner 
2928
50
c   of the father, the Born-level particle to which i_fks and j_fks are 
2929
51
c   attached. If one given particle is the colour/anticolour partner of
2930
52
c   the father in more than one colour flow, it is counted only once
2931
53
c   in ipartners(0)
2932
54
c ipartners(i), 1<=i<=nexternal-1: the label (according to Born-level
2933
55
c   labelling) of the i^th colour partner of the father
2934
56
c
2935
57
c colorflow(i,0), 1<=i<=nexternal-1: number of colour flows in which
2936
58
c   the particle ipartners(i) is a colour partner of the father
2937
59
c colorflow(i,j): the actual label (according to born_leshouche.inc)
2938
60
c   of the j^th colour flow in which the father and ipartners(i) are
2939
61
c   colour partners
2940
62
c
2941
63
c Example: in the process q(1) qbar(2) -> g(3) g(4), the two color flows are
2942
64
c
2943
65
c j=1    i    icolup(1)    icolup(2)       j=2    i    icolup(1)    icolup(2)
2944
66
c        1      500           0                   1      500           0
2945
67
c        2       0           501                  2       0           501
2946
68
c        3      500          502                  3      502          501
2947
69
c        4      502          501                  4      500          502
2948
70
c
2949
71
c and if one fixes for example fksfather=3, then ipartners = 3
2950
72
c while colorflow =  0  0  0                                 1
2951
73
c                    1  1  0                                 4
2952
74
c                    2  1  2                                 2
2953
75
c                    1  2  0
2954
76
c
2955
77
2956
78
      fksfather=min(i_fks,j_fks)
2957
79
2958
80
c isspecial will be set equal to .true. only if the father is a gluon,
2959
81
c and another gluon will be found which is connected to it by both
2960
82
c colour and anticolour
2961
83
      isspecial=.false.
2962
84
c
2963
85
      if (split_type(qcd_pos)) then
2964
86
        ! identify the color partners 
2965
87
c consider only leading colour flows
2966
88
        num_leading_cflows=0
2967
89
        do i=1,max_bcol
2968
90
          is_leading_cflow(i)=.false.
2969
91
          do j=1,mapconfig(0)
2970
92
            if(icolamp(i,j,1))then
2971
93
               is_leading_cflow(i)=.true.
2972
94
               num_leading_cflows=num_leading_cflows+1
2973
95
               exit
2974
96
            endif
2975
97
          enddo
2976
98
        enddo
2977
99
c
2978
100
        do i=1,max_bcol
2979
101
          if(.not.is_leading_cflow(i))cycle
2980
102
c Loop over Born-level colour flows
2981
103
          isspecial0=.false.
2982
104
c nglu and nsngl are the number of gluons (except for the father) and of 
2983
105
c colour singlets in the Born process, according to the information 
2984
106
c stored in ICOLUP
2985
107
          nglu=0
2986
108
          nsngl=0
2987
109
          mothercol(1)=ICOLUP(1,fksfather,i)
2988
110
          mothercol(2)=ICOLUP(2,fksfather,i)
2989
111
          notagluon=(mothercol(1).eq.0 .or. mothercol(2).eq.0)
2990
112
c
2991
113
          do j=1,nexternal-1
2992
114
c Loop over Born-level particles; j is the possible colour partner of father,
2993
115
c and whether this is the case is determined inside this loop
2994
116
            if (j.ne.fksfather) then
2995
117
c Skip father (it cannot be its own colour partner)
2996
118
               if(ICOLUP(1,j,i).eq.0.and.ICOLUP(2,j,i).eq.0)
2997
119
     #           nsngl=nsngl+1
2998
120
               if(ICOLUP(1,j,i).ne.0.and.ICOLUP(2,j,i).ne.0)
2999
121
     #           nglu=nglu+1
3000
122
               if ( (j.le.nincoming.and.fksfather.gt.nincoming) .or.
3001
123
     #              (j.gt.nincoming.and.fksfather.le.nincoming) ) then
3002
124
c father and j not both in the initial or in the final state -- connect
3003
125
c colour (1) with colour (i1(1)), and anticolour (2) with anticolour (i1(2))
3004
126
                  i1(1)=1
3005
127
                  i1(2)=2
3006
128
               else
3007
129
c father and j both in the initial or in the final state -- connect
3008
130
c colour (1) with anticolour (i1(2)), and anticolour (2) with colour (i1(1))
3009
131
                  i1(1)=2
3010
132
                  i1(2)=1
3011
133
               endif
3012
134
               do l=1,2
3013
135
c Loop over colour and anticolour of father
3014
136
                  found=.false.
3015
137
                  if( ICOLUP(i1(l),j,i).eq.mothercol(l) .and.
3016
138
     &                ICOLUP(i1(l),j,i).ne.0 ) then
3017
139
c When ICOLUP(i1(l),j,i) = mothercol(l), the colour (if i1(l)=1) or
3018
140
c the anticolour (if i1(l)=2) of particle j is connected to the
3019
141
c colour (if l=1) or the anticolour (if l=2) of the father
3020
142
                     k0=-1
3021
143
                     do k=1,ipartners(0)
3022
144
c Loop over previously-found colour/anticolour partners of father
3023
145
                        if(ipartners(k).eq.j)then
3024
146
                           if(found)then
3025
147
c Safety measure: if this condition is met, it means that there exist
3026
148
c k1 and k2 such that ipartners(k1)=ipartners(k2). This is thus a bug,
3027
149
c since ipartners() is the list of possible partners of father, where each
3028
150
c Born-level particle must appears at most once
3029
151
                              write(*,*)'Error #1 in set_matrices'
3030
152
                              write(*,*)i,j,l,k
3031
153
                              stop
3032
154
                           endif
3033
155
                           found=.true.
3034
156
                           k0=k
3035
157
                        endif
3036
158
                     enddo
3037
159
                     if (.not.found) then
3038
160
                        ipartners(0)=ipartners(0)+1
3039
161
                        ipartners(ipartners(0))=j
3040
162
                        k0=ipartners(0)
3041
163
                     endif
3042
164
c At this point, k0 is the k0^th colour/anticolour partner of father.
3043
165
c Therefore, ipartners(k0)=j
3044
166
                     if(k0.le.0.or.ipartners(k0).ne.j)then
3045
167
                        write(*,*)'Error #2 in set_matrices'
3046
168
                        write(*,*)i,j,l,k0,ipartners(k0)
3047
169
                        stop
3048
170
                     endif
3049
171
                     spec_case=l.eq.2 .and. colorflow(k0,0).ge.1 .and.
3050
172
     &                    colorflow(k0,colorflow(k0,0)).eq.i 
3051
173
                     if (.not.spec_case)then
3052
174
c Increase by one the number of colour flows in which the father is
3053
175
c (anti)colour-connected with its k0^th partner (according to the
3054
176
c list defined by ipartners)
3055
177
                        colorflow(k0,0)=colorflow(k0,0)+1
3056
178
c Store the label of the colour flow thus found
3057
179
                        colorflow(k0,colorflow(k0,0))=i
3058
180
                     elseif (spec_case)then
3059
181
c Special case: father and ipartners(k0) are both gluons, connected
3060
182
c by colour AND anticolour: the number of colour flows was overcounted
3061
183
c by one unit, so decrease it
3062
184
                         if( notagluon .or.
3063
185
     &                       ICOLUP(i1(1),j,i).eq.0 .or.
3064
186
     &                       ICOLUP(i1(2),j,i).eq.0 )then
3065
187
                            write(*,*)'Error #3 in set_matrices'
3066
188
                            write(*,*)i,j,l,k0,i1(1),i1(2)
3067
189
                            stop
3068
190
                         endif
3069
191
                         colorflow(k0,colorflow(k0,0))=i
3070
192
                         isspecial0=.true.
3071
193
                     endif
3072
194
                  endif
3073
195
               enddo
3074
196
            endif
3075
197
          enddo
3076
198
          if( ((nglu+nsngl).gt.(nexternal-2)) .or.
3077
199
     #       (isspecial0.and.(nglu+nsngl).ne.(nexternal-2)) )then
3078
200
           write(*,*)'Error #4 in set_matrices'
3079
201
           write(*,*)isspecial0,nglu,nsngl
3080
202
           stop
3081
203
          endif
3082
204
          isspecial=isspecial.or.isspecial0
3083
205
        enddo
3084
206
3085
207
      else if (split_type(qed_pos)) then
3086
208
        ! do nothing, the partner will be assigned at run-time 
3087
209
        ! (it is kinematics-dependent)
3088
210
        continue
3089
211
      endif
3090
212
      call check_mc_matrices
3091
213
      return
3092
214
      end
3093
215
3094
216
3095
217
3096
218
      subroutine check_mc_matrices
3097
219
      implicit none
3098
220
      include "nexternal.inc"
3099
221
      include "born_nhel.inc"
3100
222
c      include "fks.inc"
3101
223
      integer fks_j_from_i(nexternal,0:nexternal)
3102
224
     &     ,particle_type(nexternal),pdg_type(nexternal)
3103
225
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
3104
226
      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
3105
227
      common /MC_info/ ipartners,colorflow
3106
228
      integer i,j,ipart,iflow,ntot,ithere(1000)
3107
229
      integer fksfather
3108
230
      common/cfksfather/fksfather
3109
231
      logical notagluon
3110
232
      common/cnotagluon/notagluon
3111
233
      logical isspecial
3112
234
      common/cisspecial/isspecial
3113
235
3114
236
      include 'orders.inc'
3115
237
      logical split_type(nsplitorders) 
3116
238
      common /c_split_type/split_type
3117
239
3118
240
      logical is_leading_cflow(max_bcol)
3119
241
      integer num_leading_cflows
3120
242
      common/c_leading_cflows/is_leading_cflow,num_leading_cflows
3121
243
c
3122
244
      if(ipartners(0).gt.nexternal-1)then
3123
245
        write(*,*)'Error #1 in check_mc_matrices',ipartners(0)
3124
246
        stop
3125
247
      endif
3126
248
c
3127
249
      
3128
250
      if (split_type(QCD_pos)) then
3129
251
      ! these tests only apply for QCD-type splittings
3130
252
        do i=1,ipartners(0)
3131
253
          ipart=ipartners(i)
3132
254
          if( ipart.eq.fksfather .or.
3133
255
     #        ipart.le.0 .or. ipart.gt.nexternal-1 .or.
3134
256
     #        ( abs(particle_type(ipart)).ne.3 .and.
3135
257
     #          particle_type(ipart).ne.8 ) )then
3136
258
            write(*,*)'Error #2 in check_mc_matrices',i,ipart,
3137
259
     #  particle_type(ipart)
3138
260
            stop
3139
261
          endif
3140
262
        enddo
3141
263
c
3142
264
        do i=1,nexternal-1
3143
265
          ithere(i)=1
3144
266
        enddo
3145
267
        do i=1,ipartners(0)
3146
268
          ipart=ipartners(i)
3147
269
          ithere(ipart)=ithere(ipart)-1
3148
270
          if(ithere(ipart).lt.0)then
3149
271
            write(*,*)'Error #3 in check_mc_matrices',i,ipart
3150
272
            stop
3151
273
          endif
3152
274
        enddo
3153
275
c
3154
276
c ntot is the total number of colour plus anticolour partners of father
3155
277
        ntot=0
3156
278
        do i=1,ipartners(0)
3157
279
          ntot=ntot+colorflow(i,0)
3158
280
c
3159
281
          if( colorflow(i,0).le.0 .or.
3160
282
     #        colorflow(i,0).gt.max_bcol )then
3161
283
            write(*,*)'Error #4 in check_mc_matrices',i,colorflow(i,0)
3162
284
            stop
3163
285
          endif
3164
286
c
3165
287
          do j=1,max_bcol
3166
288
            ithere(j)=1
3167
289
          enddo
3168
290
          do j=1,colorflow(i,0)
3169
291
            iflow=colorflow(i,j)
3170
292
            ithere(iflow)=ithere(iflow)-1
3171
293
            if(ithere(iflow).lt.0)then
3172
294
              write(*,*)'Error #5 in check_mc_matrices',i,j,iflow
3173
295
              stop
3174
296
            endif
3175
297
          enddo
3176
298
c
3177
299
        enddo
3178
300
c
3179
301
        if( (notagluon.and.ntot.ne.num_leading_cflows) .or.
3180
302
     #    ( (.not.notagluon).and.
3181
303
     #      ( (.not.isspecial).and.ntot.ne.(2*num_leading_cflows) .or.
3182
304
     #        (isspecial.and.ntot.ne.num_leading_cflows) ) ) )then
3183
305
          write(*,*)'Error #6 in check_mc_matrices',
3184
306
     #     notagluon,ntot,num_leading_cflows,max_bcol
3185
307
          stop
3186
308
        endif
3187
309
c
3188
310
        if(num_leading_cflows.gt.max_bcol)then
3189
311
          write(*,*)'Error #7 in check_mc_matrices',
3190
312
     #     num_leading_cflows,max_bcol
3191
313
          stop
3192
314
        endif
3193
315
3194
316
      else if (split_type(QED_pos)) then
3195
317
        ! write here possible checks for QED-type splittings
3196
318
        continue
3197
319
      endif
3198
320
      return
3199
321
      end
3200
322
3201
323
3202
324
      subroutine find_ipartner_QED(pp, nofpartners)
3203
325
      implicit none
3204
326
      include 'nexternal.inc'
3205
327
      double precision pp(0:3, nexternal)
3206
328
      integer nofpartners
3207
329
3208
330
      integer fks_j_from_i(nexternal,0:nexternal)
3209
331
     &     ,particle_type(nexternal),pdg_type(nexternal)
3210
332
      common /c_fks_inc/fks_j_from_i,particle_type,pdg_type
3211
333
      double precision particle_charge(nexternal)
3212
334
      common /c_charges/particle_charge
3213
335
3214
336
      integer i_fks,j_fks,fksfather
3215
337
      common/fks_indices/i_fks,j_fks
3216
338
      double precision pmass(nexternal)
3217
339
      double precision zero
3218
340
      parameter (zero=0d0)
3219
341
3220
342
      include 'genps.inc'
3221
343
      include "born_nhel.inc"
3222
344
      integer idup(nexternal-1,maxproc)
3223
345
      integer mothup(2,nexternal-1,maxproc)
3224
346
      integer icolup(2,nexternal-1,max_bcol)
3225
347
      include 'born_leshouche.inc'
3226
348
      include 'coupl.inc'
3227
349
      integer ipartners(0:nexternal-1),colorflow(nexternal-1,0:max_bcol)
3228
350
      common /MC_info/ ipartners,colorflow
3229
351
c
3230
352
c     Shower MonteCarlo
3231
353
c     
3232
354
      character*10 shower_mc
3233
355
      common /cMonteCarloType/shower_mc
3234
356
3235
357
      logical found
3236
358
      logical same_state
3237
359
      double precision ppmin, ppnow
3238
360
      integer partner
3239
361
      integer i,j
3240
362
      double precision chargeprod
3241
363
      double precision dot
3242
364
3243
365
      include 'pmass.inc'
3244
366
      
3245
367
      found=.false.
3246
368
      ppmin=1d99
3247
369
      fksfather=min(i_fks,j_fks)
3248
370
3249
371
      if (shower_mc.eq.'PYTHIA8') then
3250
372
        ! this should follow what is done in TimeShower::setupQEDdip
3251
373
        ! first, look for the lowest-mass same- (opposite-)flavour pair of
3252
374
        ! particles in the opposite (same) state of the system
3253
375
        do j=1,nexternal
3254
376
          if (j.ne.fksfather.and.j.ne.i_fks) then
3255
377
            same_state = (j.gt.nincoming.and.fksfather.gt.nincoming).or.
3256
378
     $                   (j.le.nincoming.and.fksfather.le.nincoming)
3257
379
3258
380
            if ((pdg_type(j).eq.pdg_type(fksfather).and..not.same_state).or. 
3259
381
     $          (pdg_type(j).eq.-pdg_type(fksfather).and.same_state)) then
3260
382
3261
383
              ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j)
3262
384
              if (ppnow.lt.ppmin) then
3263
385
                found=.true.
3264
386
                partner=j
3265
387
              endif
3266
388
            endif
3267
389
          endif
3268
390
        enddo
3269
391
        
3270
392
        ! if no partner has been found, then look for the
3271
393
        ! lowest-mass/chargeprod pair
3272
394
        if (.not.found) then
3273
395
          do j=1,nexternal
3274
396
            if (j.ne.fksfather.and.j.ne.i_fks) then
3275
397
              if (particle_charge(fksfather).ne.0d0.and.particle_charge(j).ne.0d0) then
3276
398
                ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j) / 
3277
399
     $            (particle_charge(fksfather) * particle_charge(j))
3278
400
                if (ppnow.lt.ppmin) then
3279
401
                  found=.true.
3280
402
                  partner=j
3281
403
                endif
3282
404
              endif
3283
405
            endif
3284
406
          enddo
3285
407
        endif
3286
408
3287
409
        ! if no partner has been found, then look for the
3288
410
        ! lowest-mass pair
3289
411
        if (.not.found) then
3290
412
          do j=1,nexternal
3291
413
            if (j.ne.fksfather.and.j.ne.i_fks) then
3292
414
              ppnow=dot(pp(0,fksfather),pp(0,j)) - pmass(fksfather)*pmass(j) 
3293
415
              if (ppnow.lt.ppmin) then
3294
416
                found=.true.
3295
417
                partner=j
3296
418
              endif
Reviewer	Date Requested	Status
Olivier Mattelaer	2018-10-04	Approve on 2018-12-19
marco zaro	2018-10-04	Needs Fixing on 2018-11-21
Hua-Sheng Shao	2018-10-04	Pending
Valentin Hirschi	2018-10-04	Pending
davide.pagani.85	2018-10-04	Pending
Review via email: mp+356133@code.launchpad.net