MadGraph5_aMC@NLO

Merge lp:~maddevelopers/mg5amcnlo/LTS_dev into lp:mg5amcnlo/lts

  1. LTS_dev
  2. Merge into series2.0
Proposed by Olivier Mattelaer
Status: Merged
Merged at revision: 299
Proposed branch: lp:~maddevelopers/mg5amcnlo/LTS_dev
Merge into: lp:mg5amcnlo/lts
Diff against target: 4903 lines (+2027/-1902) (has conflicts)
57 files modified
MadSpin/madspin (+8/-6)
Template/LO/Source/run.inc (+2/-2)
Template/LO/Source/run_printout.f (+10/-3)
Template/LO/Source/rw_events.short.f (+0/-2)
Template/LO/SubProcesses/cuts.f (+1/-17)
Template/LO/SubProcesses/initcluster.f (+3/-2)
Template/LO/SubProcesses/reweight.f (+35/-34)
Template/LO/SubProcesses/setcuts.f (+2/-1)
Template/LO/SubProcesses/setscales.f (+9/-9)
Template/LO/bin/internal/Gridpack/gridrun (+1/-0)
Template/MadWeight/src/run.inc (+3/-2)
Template/MadWeight/src/setrun.f (+2/-1)
UpdateNotes.txt (+18/-0)
VERSION (+5/-0)
aloha/__init__.py (+1/-1)
bin/compile.py (+10/-2)
bin/create_aloha_release.py (+1/-1)
madgraph/core/base_objects.py (+13/-0)
madgraph/core/diagram_generation.py (+4/-2)
madgraph/interface/common_run_interface.py (+7/-4)
madgraph/interface/loop_interface.py (+1/-1)
madgraph/interface/madevent_interface.py (+9/-5)
madgraph/interface/madgraph_interface.py (+7/-7)
madgraph/interface/reweight_interface.py (+16/-12)
madgraph/iolibs/drawing_eps.py (+1/-1)
madgraph/iolibs/export_fks.py (+4/-4)
madgraph/iolibs/export_v4.py (+6/-6)
madgraph/iolibs/file_writers.py (+1/-1)
madgraph/iolibs/files.py (+1/-1)
madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc (+1/-1)
madgraph/loop/loop_exporters.py (+2/-2)
madgraph/madweight/write_MadWeight.py (+1/-1)
madgraph/various/banner.py (+45/-2)
madgraph/various/misc.py (+6/-0)
madgraph/various/process_checks.py (+3/-1)
mg5decay/decay_objects.py (+12/-9)
models/build_restriction_lib.py (+1/-1)
models/usermod.py (+1/-0)
tests/__init__.py (+1/-1)
tests/acceptance_tests/test_cmd_madevent.py (+2/-2)
tests/acceptance_tests/test_madspin.py (+1/-1)
tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f (+1/-1)
tests/input_files/IOTestsComparison/TestMadWeight/modification_to_cuts/cuts.f (+1723/-1739)
tests/input_files/full_sm/particles.py (+1/-1)
tests/input_files/mssm/particles.py (+1/-1)
tests/input_files/sm/particles.py (+1/-1)
tests/parallel_tests/decay_comparator.py (+1/-1)
tests/parallel_tests/me_comparator.py (+1/-1)
tests/parallel_tests/test_ML5MSSMQCD.py (+1/-1)
tests/test_manager.py (+2/-2)
tests/unit_tests/__init__.py (+1/-1)
tests/unit_tests/loop/test_loop_exporters.py (+1/-1)
tests/unit_tests/madweight/test_permutation.py (+1/-1)
tests/unit_tests/various/test_banner.py (+32/-0)
tests/unit_tests/various/test_cmd.py (+1/-1)
tests/unit_tests/various/test_misc.py (+1/-1)
tests/unit_tests/various/test_rambo.py (+1/-1)
Text conflict in UpdateNotes.txt
Text conflict in VERSION
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/LTS_dev
Reviewer Review Type Date Requested Status
MadTeam Pending
Review via email: mp+407493@code.launchpad.net

Description of the change

This is pure bug fixing for our long term stable release (2.9.x)
This contains some quite heavy change on the dynamical scale choice.
This will be released on Monday

To post a comment you must log in.

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'MadSpin/madspin'
2--- MadSpin/madspin 2018-05-10 15:14:33 +0000
3+++ MadSpin/madspin 2021-08-21 21:19:44 +0000
4@@ -17,6 +17,8 @@
5 """This is the main executable, a simple frontend to set up the PYTHONPATH
6 and call immediately the command line interface scripts"""
7
8+from __future__ import absolute_import
9+from __future__ import print_function
10 import sys
11 if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
12 sys.exit('MadSpin works only with python 2.6 or later (but not python 3.X).\n\
13@@ -68,7 +70,7 @@
14 try:
15 import pyreadline as readline
16 except:
17- print "For tab completion and history, install module readline."
18+ print("For tab completion and history, install module readline.")
19 else:
20 import rlcompleter
21
22@@ -99,7 +101,7 @@
23 pass
24
25 if __debug__:
26- print 'Running MG5 in debug mode'
27+ print('Running MG5 in debug mode')
28
29
30 # Set logging level according to the logging level given by options
31@@ -130,7 +132,7 @@
32 cmd_line.import_command_file(input_file)
33 cmd_line.run_cmd('quit')
34 else:
35- print "using input+file", input_file
36+ print("using input+file", input_file)
37 cmd_line = interface.MadSpinInterface()
38 cmd_line.use_rawinput = False
39 cmd_line.haspiping = False
40@@ -139,7 +141,7 @@
41 else:
42 # Interactive mode
43 if options.web:
44- if not os.environ.has_key('MADGRAPH_DATA'):
45+ if 'MADGRAPH_DATA' not in os.environ:
46 os.environ['MADGRAPH_DATA'] = root_path
47 os.environ['MADGRAPH_BASE'] = os.path.join(root_path,'input')
48 os.environ['REMOTE_USER'] = 'webmode'
49@@ -158,10 +160,10 @@
50 if not os.path.exists(os.path.join(os.environ['HOME'], '.mg5')):
51 os.mkdir(os.path.join(os.environ['HOME'], '.mg5'))
52 readline.write_history_file(history_file)
53- except Exception, error:
54+ except Exception as error:
55 pass
56 except KeyboardInterrupt:
57- print 'writting history and quit on KeyboardInterrupt'
58+ print('writting history and quit on KeyboardInterrupt')
59 pass
60
61
62
63=== modified file 'Template/LO/Source/run.inc'
64--- Template/LO/Source/run.inc 2020-10-01 14:23:34 +0000
65+++ Template/LO/Source/run.inc 2021-08-21 21:19:44 +0000
66@@ -5,9 +5,9 @@
67 c Scales
68 c
69 real*8 scale,scalefact,alpsfact
70- logical fixed_ren_scale,fixed_fac_scale,fixed_couplings,hmult
71+ logical fixed_ren_scale,fixed_fac_scale1,fixed_fac_scale2,fixed_couplings,hmult
72 integer ickkw,nhmult,asrwgtflavor, dynamical_scale_choice
73- common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale,
74+ common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale1,fixed_fac_scale2,
75 $ fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor,
76 $ dynamical_scale_choice
77 c
78
79=== modified file 'Template/LO/Source/run_printout.f'
80--- Template/LO/Source/run_printout.f 2021-02-12 15:28:58 +0000
81+++ Template/LO/Source/run_printout.f 2021-08-21 21:19:44 +0000
82@@ -55,11 +55,18 @@
83 else
84 write(6,*) 'Renormalization scale set on event-by-event basis'
85 endif
86- if(fixed_fac_scale) then
87+ if(fixed_fac_scale1.and.fixed_fac_scale2) then
88 write(6,*) 'Factorization scales fixed @ ',
89- & dsqrt(q2fact(1)),dsqrt(q2fact(2))
90- else
91+ & dsqrt(q2fact(1)),dsqrt(q2fact(2))
92+ else if(.not.fixed_fac_scale1.and..not.fixed_fac_scale2) then
93 write(6,*) 'Factorization scale set on event-by-event basis'
94+ else if(fixed_fac_scale1) then
95+ write(6,*) 'Factorization scales fixed for beam1 @ ',
96+ & dsqrt(q2fact(1)),dsqrt(q2fact(2))
97+ else
98+ write(6,*) 'Factorization scales fixed for beam2 @ ',
99+ & dsqrt(q2fact(1)),dsqrt(q2fact(2))
100+
101 endif
102
103 write(6,*)
104
105=== modified file 'Template/LO/Source/rw_events.short.f'
106--- Template/LO/Source/rw_events.short.f 2010-10-30 03:26:37 +0000
107+++ Template/LO/Source/rw_events.short.f 2021-08-21 21:19:44 +0000
108@@ -35,8 +35,6 @@
109 c
110 c include 'coupl.inc'
111 c real*8 scale
112-c logical fixed_ren_scale,fixed_fac_scale
113-c common/to_scale/scale,fixed_ren_scale,fixed_fac_scale
114
115 logical banner_open
116 integer lun_ban
117
118=== modified file 'Template/LO/SubProcesses/cuts.f'
119--- Template/LO/SubProcesses/cuts.f 2020-05-13 09:45:34 +0000
120+++ Template/LO/SubProcesses/cuts.f 2021-08-21 21:19:44 +0000
121@@ -256,22 +256,6 @@
122 c endif
123 c endif
124
125-c if(.not.fixed_fac_scale) then
126-c call set_fac_scale(P,q2fact)
127-c if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
128-c write(6,*)
129-c write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
130-c write(6,*) ' Dynamical renormalization scale choice '
131-c write(6,*) ' selected but user subroutine'
132-c write(6,*) ' set_fac_scale not edited in file:setpara.f'
133-c write(6,*) ' Switching to a fixed_fac_scale choice'
134-c write(6,*) ' with q2fact(i)=zmass**2'
135-c fixed_fac_scale=.true.
136-c q2fact(1)=91.2d0**2
137-c q2fact(2)=91.2d0**2
138-c write(6,*) 'scales=',q2fact(1),q2fact(2)
139-c endif
140-c endif
141
142 if(fixed_ren_scale) then
143 G = SQRT(4d0*PI*ALPHAS(scale))
144@@ -1237,7 +1221,7 @@
145 if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
146 endif
147
148- if(.not.fixed_fac_scale) then
149+ if(.not.fixed_fac_scale1.or..not.fixed_fac_scale2) then
150 call set_fac_scale(P,q2fact)
151 endif
152
153
154=== modified file 'Template/LO/SubProcesses/initcluster.f'
155--- Template/LO/SubProcesses/initcluster.f 2012-03-22 06:16:35 +0000
156+++ Template/LO/SubProcesses/initcluster.f 2021-08-21 21:19:44 +0000
157@@ -24,7 +24,7 @@
158 c
159
160 c if (ickkw.le.0) return
161- if (ickkw.le.0.and.xqcut.le.0d0.and.fixed_ren_scale.and.fixed_fac_scale) return
162+ if (ickkw.le.0.and.xqcut.le.0d0.and.fixed_ren_scale.and.fixed_fac_scale1.and.fixed_fac_scale2) return
163
164 c if(ickkw.eq.2.and.xqcut.le.0d0)then
165 c write(*,*)'Must set qcut > 0 for ickkw = 2'
166@@ -38,7 +38,8 @@
167 c q2fact(1) = scale**2 ! fact scale**2 for pdf1
168 c q2fact(2) = scale**2 ! fact scale**2 for pdf2
169 c fixed_ren_scale=.true.
170-c fixed_fac_scale=.true.
171+c fixed_fac_scale1=.true.
172+c fixed_fac_scale2=.true.
173 c endif
174 c
175 c initialize clustering map
176
177=== modified file 'Template/LO/SubProcesses/reweight.f'
178--- Template/LO/SubProcesses/reweight.f 2020-09-07 22:17:59 +0000
179+++ Template/LO/SubProcesses/reweight.f 2021-08-21 21:19:44 +0000
180@@ -628,7 +628,7 @@
181 external is_octet
182 setclscales=.true.
183
184- if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.scale.gt.0) then
185+ if(ickkw.le.0.and.xqcut.le.0d0.and.q2fact(1).gt.0.and.q2fact(2).gt.0.and.scale.gt.0) then
186 if(use_syst)then
187 s_scale=scale
188 n_qcd=nqcd(iconfig)
189@@ -689,10 +689,8 @@
190 q2bck(2)=q2fact(2)
191 first=.false.
192 else if(ickkw.gt.0) then
193- if(fixed_fac_scale) then
194- q2fact(1)=q2bck(1)
195- q2fact(2)=q2bck(2)
196- endif
197+ if(fixed_fac_scale1) q2fact(1)=q2bck(1)
198+ if (fixed_fac_scale2) q2fact(2)=q2bck(2)
199 endif
200
201 c Preparing graph particle information (ipart, needed to keep track of
202@@ -1081,8 +1079,10 @@
203 endif
204 endif
205
206- if(ickkw.eq.0.and.(fixed_fac_scale.or.q2fact(1).gt.0).and.
207- $ (fixed_ren_scale.or.scale.gt.0)) return
208+ if(ickkw.eq.0
209+ $ .and. (fixed_fac_scale1.or.q2fact(1).gt.0)
210+ $ .and. (fixed_fac_scale2.or.q2fact(2).gt.0)
211+ $ .and. (fixed_ren_scale.or.scale.gt.0)) return
212
213 c Ensure that last scales are at least as big as first scales
214 if(jlast(1).gt.0)
215@@ -1090,34 +1090,36 @@
216 if(jlast(2).gt.0)
217 $ pt2ijcl(jlast(2))=max(pt2ijcl(jlast(2)),pt2ijcl(jfirst(2)))
218
219- if(ickkw.gt.0.and.q2fact(1).gt.0) then
220+ if(ickkw.gt.0.and.q2fact(1).gt.0.and.q2fact(2).gt.0) then
221 c Use the fixed or previously set scale for central scale
222 if(jcentral(1).gt.0) pt2ijcl(jcentral(1))=q2fact(1)
223 if(jcentral(2).gt.0.and.jcentral(2).ne.jcentral(1))
224 $ pt2ijcl(jcentral(2))=q2fact(2)
225 endif
226
227- if(nexternal.eq.3.and.nincoming.eq.2.and.q2fact(1).eq.0) then
228- q2fact(1)=pt2ijcl(nexternal-2)
229- q2fact(2)=pt2ijcl(nexternal-2)
230+ if(nexternal.eq.3.and.nincoming.eq.2.and.(q2fact(1).eq.0.or.q2fact(2).eq.0)) then
231+ if(.not.fixed_fac_scale1) q2fact(1)=pt2ijcl(nexternal-2)
232+ if(.not.fixed_fac_scale2) q2fact(2)=pt2ijcl(nexternal-2)
233 endif
234
235- if(q2fact(1).eq.0d0) then
236+ if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
237 c Use the geom. average of central scale and first non-radiation vertex
238- if(jlast(1).gt.0) q2fact(1)=sqrt(pt2ijcl(jlast(1))*pt2ijcl(jcentral(1)))
239- if(jlast(2).gt.0) q2fact(2)=sqrt(pt2ijcl(jlast(2))*pt2ijcl(jcentral(2)))
240+ if(jlast(1).gt.0.and..not.fixed_fac_scale1) q2fact(1)=sqrt(pt2ijcl(jlast(1))*pt2ijcl(jcentral(1)))
241+ if(jlast(2).gt.0.and..not.fixed_fac_scale2) q2fact(2)=sqrt(pt2ijcl(jlast(2))*pt2ijcl(jcentral(2)))
242 if(jcentral(1).gt.0.and.jcentral(1).eq.jcentral(2))then
243 c We have a qcd line going through the whole event, use single scale
244- q2fact(1)=max(q2fact(1),q2fact(2))
245- q2fact(2)=q2fact(1)
246+ if(.not.fixed_fac_scale1.and..not.fixed_fac_scale2) then
247+ q2fact(1)=max(q2fact(1),q2fact(2))
248+ q2fact(2)=q2fact(1)
249+ endif
250 endif
251 endif
252- if(.not. fixed_fac_scale) then
253- q2fact(1)=scalefact**2*q2fact(1)
254- q2fact(2)=scalefact**2*q2fact(2)
255+ if(.not. fixed_fac_scale1.or. fixed_fac_scale2) then
256+ if(.not.fixed_fac_scale1) q2fact(1)=scalefact**2*q2fact(1)
257+ if(.not.fixed_fac_scale2) q2fact(2)=scalefact**2*q2fact(2)
258 if (.not.keepq2bck)then
259- q2bck(1)=q2fact(1)
260- q2bck(2)=q2fact(2)
261+ if(.not.fixed_fac_scale1) q2bck(1)=q2fact(1)
262+ if(.not.fixed_fac_scale2) q2bck(2)=q2fact(2)
263 endif
264 if (btest(mlevel,3))
265 $ write(*,*) 'Set central fact scales to ',sqrt(q2bck(1)),sqrt(q2bck(2))
266@@ -1155,30 +1157,33 @@
267
268 c Take care of case when jcentral are zero
269 if(jcentral(1).eq.0.and.jcentral(2).eq.0)then
270- if(q2fact(1).gt.0)then
271+ if(q2fact(1).gt.0.and..not.fixed_fac_scale1)then
272 pt2ijcl(nexternal-2)=q2fact(1)
273 if(nexternal.gt.3) pt2ijcl(nexternal-3)=q2fact(1)
274+ else if (q2fact(2).gt.0.and..not.fixed_fac_scale2)then
275+ pt2ijcl(nexternal-2)=q2fact(2)
276+ if(nexternal.gt.3) pt2ijcl(nexternal-3)=q2fact(2)
277 else
278- q2fact(1)=scalefact**2*pt2ijcl(nexternal-2)
279- q2fact(2)=scalefact**2*q2fact(1)
280+ if(.not.fixed_fac_scale1) q2fact(1)=scalefact**2*pt2ijcl(nexternal-2)
281+ if(.not.fixed_fac_scale2) q2fact(2)=scalefact**2*q2fact(1)
282 endif
283 elseif(jcentral(1).eq.0)then
284- q2fact(1) = scalefact**2*pt2ijcl(jfirst(1))
285+ if(.not.fixed_fac_scale1) q2fact(1) = scalefact**2*pt2ijcl(jfirst(1))
286 elseif(jcentral(2).eq.0)then
287- q2fact(2) = scalefact**2*pt2ijcl(jfirst(2))
288+ if(.not.fixed_fac_scale2) q2fact(2) = scalefact**2*pt2ijcl(jfirst(2))
289 elseif(ickkw.eq.2.or.(pdfwgt.and.ickkw.gt.0))then
290 c Total pdf weight is f1(x1,pt2E)*fj(x1*z,Q)/fj(x1*z,pt2E)
291 c f1(x1,pt2E) is given by DSIG, just need to set scale.
292 c Use the minimum scale found for fact scale in ME
293- if(jlast(1).gt.0.and.jfirst(1).le.jlast(1))
294+ if(jlast(1).gt.0.and.jfirst(1).le.jlast(1).and..not.fixed_fac_scale1)
295 $ q2fact(1)=scalefact**2*min(pt2ijcl(jfirst(1)),q2fact(1))
296- if(jlast(2).gt.0.and.jfirst(2).le.jlast(2))
297+ if(jlast(2).gt.0.and.jfirst(2).le.jlast(2).and..not.fixed_fac_scale2)
298 $ q2fact(2)=scalefact**2*min(pt2ijcl(jfirst(2)),q2fact(2))
299 endif
300
301 c Check that factorization scale is >= 2 GeV
302- if(lpp(1).ne.0.and.q2fact(1).lt.4d0.or.
303- $ lpp(2).ne.0.and.q2fact(2).lt.4d0)then
304+ if(lpp(1).ne.0.and.(q2fact(1).lt.4d0.and..not.fixed_fac_scale1).or.
305+ $ lpp(2).ne.0.and.(q2fact(2).lt.4d0.and..not.fixed_fac_scale2))then
306 if(nwarning.le.10) then
307 nwarning=nwarning+1
308 write(*,*) 'Warning: Too low fact scales: ',
309@@ -1743,10 +1748,6 @@
310 if (btest(mlevel,3))
311 $ write(*,*)' set fact scales for PS to ',
312 $ sqrt(q2fact(1)),sqrt(q2fact(2))
313- else if (abs(lpp(1)).ge.2.and.abs(lpp(1)).le.4) then
314- q2fact(1)=q2bck(1)
315- else if (abs(lpp(2)).ge.2.or.abs(lpp(2)).le.4) then
316- q2fact(2)=q2bck(2)
317 endif
318
319 if (btest(mlevel,3)) then
320
321=== modified file 'Template/LO/SubProcesses/setcuts.f'
322--- Template/LO/SubProcesses/setcuts.f 2020-09-07 22:17:59 +0000
323+++ Template/LO/SubProcesses/setcuts.f 2021-08-21 21:19:44 +0000
324@@ -142,7 +142,8 @@
325 scale=pmass(1)
326 fixed_ren_scale=.true.
327 endif
328- fixed_fac_scale=.true.
329+ fixed_fac_scale1=.true.
330+ fixed_fac_scale2=.true.
331 use_syst=.false.
332 endif
333
334
335=== modified file 'Template/LO/SubProcesses/setscales.f'
336--- Template/LO/SubProcesses/setscales.f 2021-03-16 12:01:38 +0000
337+++ Template/LO/SubProcesses/setscales.f 2021-08-21 21:19:44 +0000
338@@ -126,7 +126,7 @@
339 integer i
340 logical first
341 data first/.true./
342-
343+ double precision tempscale
344 c----------
345 c start
346 c----------
347@@ -134,8 +134,8 @@
348 if (dynamical_scale_choice.eq.-1) then
349 c Cluster external states until reducing the system to a 2->2 topology whose transverse mass is used for setting the scale.
350 c This is not done in this file due to the clustering.
351- q2factorization(1)=0d0 !factorization scale**2 for pdf1
352- q2factorization(2)=0d0 !factorization scale**2 for pdf2
353+ if(.not.fixed_fac_scale1) q2factorization(1)=0d0 !factorization scale**2 for pdf1
354+ if(.not.fixed_fac_scale2) q2factorization(2)=0d0 !factorization scale**2 for pdf2
355 elseif(dynamical_scale_choice.eq.0) then
356 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
357 cc USER DEFINE SCALE: ENTER YOUR CODE HERE cc
358@@ -143,9 +143,9 @@
359 cc dymamical_scale_choice to 0 in the run_card cc
360 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
361 c default: use the renormalization scale
362- call set_ren_scale(P,q2factorization(1))
363- q2factorization(1)=q2factorization(1)**2
364- q2factorization(2)=q2factorization(1) !factorization scale**2 for pdf2
365+ call set_ren_scale(P,tempscale)
366+ if(.not.fixed_fac_scale1) q2factorization(1)=tempscale**2
367+ if(.not.fixed_fac_scale2) q2factorization(2)=tempscale**2 !factorization scale**2 for pdf2
368
369 c
370 c-some examples of dynamical scales
371@@ -181,9 +181,9 @@
372 cc USER DEFINE SCALE: END of USER CODE cc
373 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
374 else
375- call set_ren_scale(P,q2factorization(1))
376- q2factorization(1)=q2factorization(1)**2
377- q2factorization(2)=q2factorization(1) !factorization scale**2 for pdf2
378+ call set_ren_scale(P,tempscale)
379+ if(.not.fixed_fac_scale1) q2factorization(1)=tempscale**2
380+ if(.not.fixed_fac_scale2) q2factorization(2)=tempscale**2 !factorization scale**2 for pdf2
381 endif
382
383
384
385=== modified file 'Template/LO/bin/internal/Gridpack/gridrun'
386--- Template/LO/bin/internal/Gridpack/gridrun 2021-05-31 07:36:21 +0000
387+++ Template/LO/bin/internal/Gridpack/gridrun 2021-08-21 21:19:44 +0000
388@@ -63,6 +63,7 @@
389
390 import logging
391 import logging.config
392+import internal.coloring_logging
393
394 try:
395 import psyco
396
397=== modified file 'Template/MadWeight/src/run.inc'
398--- Template/MadWeight/src/run.inc 2018-04-27 21:01:37 +0000
399+++ Template/MadWeight/src/run.inc 2021-08-21 21:19:44 +0000
400@@ -5,9 +5,10 @@
401 c Scales
402 c
403 real*8 scale,scalefact,alpsfact
404- logical fixed_ren_scale,fixed_fac_scale,fixed_couplings,hmult
405+ logical fixed_ren_scale,fixed_fac_scale1,fixed_fac_scale2
406+ logical fixed_couplings,hmult
407 integer ickkw,nhmult,asrwgtflavor,dynamical_scale_choice
408- common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale,
409+ common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale1,fixed_fac_scale2,
410 $ fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor,dynamical_scale_choice
411 c
412 c Collider
413
414=== modified file 'Template/MadWeight/src/setrun.f'
415--- Template/MadWeight/src/setrun.f 2018-04-27 21:01:37 +0000
416+++ Template/MadWeight/src/setrun.f 2021-08-21 21:19:44 +0000
417@@ -309,7 +309,8 @@
418 c user subroutines set_ren_scale and set_fac_scale in setpara.f
419
420 call get_logical(npara,param,value," fixed_ren_scale ",fixed_ren_scale,.true.)
421- call get_logical(npara,param,value," fixed_fac_scale ",fixed_fac_scale,.true.)
422+ call get_logical(npara,param,value," fixed_fac_scale1 ",fixed_fac_scale1,.true.)
423+ call get_logical(npara,param,value," fixed_fac_scale2 ",fixed_fac_scale2,.true.)
424 call get_real (npara,param,value," scale " ,scale,91.188d0)
425 call get_real (npara,param,value," dsqrt_q2fact1 " ,sf1 ,91.188d0)
426 call get_real (npara,param,value," dsqrt_q2fact2 " ,sf2 ,91.188d0)
427
428=== modified file 'UpdateNotes.txt'
429--- UpdateNotes.txt 2021-05-31 08:43:33 +0000
430+++ UpdateNotes.txt 2021-08-21 21:19:44 +0000
431@@ -1,6 +1,24 @@
432 Update notes for MadGraph5_aMC@NLO (in reverse time order)
433
434+<<<<<<< TREE
435 2.9.4(30/05/21)
436+=======
437+2.9.5 (22/08/21)
438+ OM+LM: [LO only] Fix the factorization scale dependence for lpp=2/3/4.
439+ This was claimed to be using fixed scale computation while in s\
440+ome case the scale was dynamical
441+ To have full flexibility we introduced two additional (hidden) \
442+parameter
443+ fixed_fac_scale1 and fixed_fac_scale2 that allow to choose fixe\
444+d scale for only one beam.
445+ OM: fix auto-width that was not following run_mode specification
446+ OM: fixing missing rwgt_info for reweighting with gridpack mode
447+ OM: fixing 'check' command with the skip_event mode
448+ OM: fixing auto-width computation for 3 body decay and identical particle which was sometimes leading to crash
449+ OM: Fix some potential infinite loop when running with python3
450+
451+2.9.4(28/05/21)
452+>>>>>>> MERGE-SOURCE
453 OM: Fix a python3 issue for madSpin when using in a gridpack mode (set ms_dir)
454 OM: Fix an issue for non positive definite matrix-element when using the
455 "set group_subprocesses False" mode of MG5aMc
456
457=== modified file 'VERSION'
458--- VERSION 2021-05-31 08:43:33 +0000
459+++ VERSION 2021-08-21 21:19:44 +0000
460@@ -1,5 +1,10 @@
461+<<<<<<< TREE
462 version = 2.9.4
463 date = 2021-05-30
464+=======
465+version = 2.9.5
466+date = 2021-08-22
467+>>>>>>> MERGE-SOURCE
468
469
470
471
472=== modified file 'aloha/__init__.py'
473--- aloha/__init__.py 2015-10-01 16:00:08 +0000
474+++ aloha/__init__.py 2021-08-21 21:19:44 +0000
475@@ -4,4 +4,4 @@
476 mp_precision = False # Tag for passing parameter in quadruple precision
477 aloha_prefix = 'mdl_'
478
479-class ALOHAERROR(Exception): pass
480\ No newline at end of file
481+class ALOHAERROR(Exception): pass
482
483=== modified file 'bin/compile.py'
484--- bin/compile.py 2021-01-26 14:28:08 +0000
485+++ bin/compile.py 2021-08-21 21:19:44 +0000
486@@ -18,6 +18,7 @@
487 import os
488 import sys
489 import logging
490+import logging.config
491 import time
492 import shutil
493 import subprocess
494@@ -33,11 +34,18 @@
495 import aloha.create_aloha as create_aloha
496 import madgraph.iolibs.files as files
497 import madgraph.various.misc as misc
498+import madgraph.interface.coloring_logging
499 import re
500
501 # Set logging level to error
502-logging.basicConfig(level=vars(logging)['INFO'],
503- format="%(message)s")
504+#logging.basicConfig(level=vars(logging)['INFO'],
505+# format="%(message)s")
506+level=50
507+logging.config.fileConfig(os.path.join(root_path, 'madgraph', 'interface', '.mg5_logging.conf'))
508+logging.root.setLevel(level)
509+logging.getLogger('madgraph').setLevel(level)
510+logging.getLogger('madevent').setLevel(level)
511+
512 pjoin = os.path.join
513
514 class Compile_MG5:
515
516=== modified file 'bin/create_aloha_release.py'
517--- bin/create_aloha_release.py 2019-04-17 14:39:47 +0000
518+++ bin/create_aloha_release.py 2021-08-21 21:19:44 +0000
519@@ -126,7 +126,7 @@
520 os.system('rm -rf %s ' % os.path.join(filepath,fname))
521
522 os.mkdir(filepath+'/vendor')
523-shutil.copytree(os.path.join(MG5DIR,'vendor','ply'),filepath+'/vendor/ply')
524+misc.copytree(os.path.join(MG5DIR,'vendor','ply'),filepath+'/vendor/ply')
525 files_routines.cp(MG5DIR +'/vendor/__init__.py', filepath+'/vendor/__init__.py')
526
527
528
529=== modified file 'madgraph/core/base_objects.py'
530--- madgraph/core/base_objects.py 2021-01-28 23:08:28 +0000
531+++ madgraph/core/base_objects.py 2021-08-21 21:19:44 +0000
532@@ -2391,6 +2391,8 @@
533 """Returns a nicely formatted string of the diagram content."""
534
535 pass_sanity = True
536+ removed_index = set()
537+
538 if self['vertices']:
539 mystr = '('
540 for vert in self['vertices']:
541@@ -2406,10 +2408,21 @@
542 if __debug__ and len(used_leg) != len(set(used_leg)):
543 pass_sanity = False
544 responsible = id(vert)
545+ if __debug__ and any(l['number'] in removed_index for l in vert['legs']):
546+ pass_sanity = False
547+ responsible = id(vert)
548
549 if self['vertices'].index(vert) < len(self['vertices']) - 1:
550 # Do not want ">" in the last vertex
551 mystr = mystr[:-1] + '>'
552+ if __debug__:
553+ if vert['legs'][-1]['number'] != min([l['number'] for l in vert['legs'][:-1]]):
554+ pass_sanity = False
555+ responsible = id(vert)
556+ for l in vert['legs']:
557+ removed_index.add(l['number'])
558+ removed_index.remove(vert['legs'][-1]['number'])
559+
560 lastleg = vert['legs'][-1]
561 if lastleg['polarization']:
562 mystr = mystr + str(lastleg['number']) + '(%s{%s})' % (str(lastleg['id']), lastleg['polarization']) + ','
563
564=== modified file 'madgraph/core/diagram_generation.py'
565--- madgraph/core/diagram_generation.py 2021-01-21 13:01:31 +0000
566+++ madgraph/core/diagram_generation.py 2021-08-21 21:19:44 +0000
567@@ -385,10 +385,12 @@
568 try:
569 return self.vertex_id[0] < other.vertex_id[0]
570 except TypeError as error:
571- if error.args == "'<' not supported between instances of 'tuple' and 'str'":
572+ if error.args == ("'<' not supported between instances of 'tuple' and 'str'",):
573 return False
574- else:
575+ elif error.args == ("'<' not supported between instances of 'str' and 'tuple'",):
576 return True
577+ else:
578+ raise Exception
579
580
581 for i, link in enumerate(self.links):
582
583=== modified file 'madgraph/interface/common_run_interface.py'
584--- madgraph/interface/common_run_interface.py 2021-05-26 12:28:50 +0000
585+++ madgraph/interface/common_run_interface.py 2021-08-21 21:19:44 +0000
586@@ -1799,9 +1799,12 @@
587 return
588
589 if self.options['run_mode'] ==2 and self.options['nb_core'] != 1:
590- nb_submit = min(self.options['nb_core'], nb_event//2500)
591+ nb_submit = min(int(self.options['nb_core']), nb_event//2500)
592 elif self.options['run_mode'] ==1:
593- nb_submit = min(self.options['cluster_size'], nb_event//25000)
594+ try:
595+ nb_submit = min(int(self.options['cluster_size']), nb_event//25000)
596+ except Exception:
597+ nb_submit =1
598 else:
599 nb_submit =1
600
601@@ -1934,7 +1937,7 @@
602
603 multicore = True
604 if self.options['run_mode'] in [0,1]:
605- multicore = False
606+ return False
607
608 lines = [l.strip() for l in open(card) if not l.strip().startswith('#')]
609 while lines and not lines[0].startswith('launch'):
610@@ -3238,7 +3241,7 @@
611 elif os.path.exists(pjoin(self.me_dir, args[1])):
612 self.options[args[0]] = pjoin(self.me_dir, args[1])
613 else:
614- raise self.InvalidCmd('Not a valid path: keep previous value: \'%s\'' % self.options[args[0]])
615+ raise self.InvalidCmd('Not a valid path: keep previous value: \'%s\' for %s instead of %s' % (self.options[args[0]], args[0], args[1]) )
616 else:
617 self.options[args[0]] = args[1]
618
619
620=== modified file 'madgraph/interface/loop_interface.py'
621--- madgraph/interface/loop_interface.py 2020-09-22 06:54:43 +0000
622+++ madgraph/interface/loop_interface.py 2021-08-21 21:19:44 +0000
623@@ -963,7 +963,7 @@
624 if os.path.exists(pjoin(install_dir1, 'collier')):
625 self.code['collier'] = pjoin(install_dir1, 'collier')
626 if os.path.exists(pjoin(install_dir2, 'golem95')):
627- self.code['glem'] = pjoin(install_dir2, 'golem95')
628+ self.code['golem'] = pjoin(install_dir2, 'golem95')
629 if os.path.exists(pjoin(install_dir1, 'ninja')):
630 self.code['ninja'] = pjoin(install_dir2, 'ninja','lib')
631
632
633=== modified file 'madgraph/interface/madevent_interface.py'
634--- madgraph/interface/madevent_interface.py 2021-05-25 20:35:55 +0000
635+++ madgraph/interface/madevent_interface.py 2021-08-21 21:19:44 +0000
636@@ -2992,8 +2992,10 @@
637 eradir = self.options['exrootanalysis_path']
638 if eradir and misc.is_executable(pjoin(eradir,'ExRootLHEFConverter')):
639 self.update_status("Create Root file", level='parton')
640- misc.gunzip('%s/%s/unweighted_events.lhe.gz' %
641- (pjoin(self.me_dir,'Events'), self.run_name))
642+ path = '%s/%s/unweighted_events.lhe.gz' % (pjoin(self.me_dir,'Events'), self.run_name)
643+
644+ if os.path.exists(path):
645+ misc.gunzip(path)
646
647 self.create_root_file('%s/unweighted_events.lhe' % self.run_name,
648 '%s/unweighted_events.root' % self.run_name)
649@@ -3108,7 +3110,8 @@
650 else:
651 run_card = self.run_card
652 self.run_card = run_card
653- self.cluster.modify_interface(self)
654+ if self.cluster:
655+ self.cluster.modify_interface(self)
656 if self.ninitial == 1:
657 run_card['lpp1'] = 0
658 run_card['lpp2'] = 0
659@@ -3129,7 +3132,7 @@
660 if not os.path.isfile(pjoin(run_card['bias_module'],mandatory_file)):
661 raise InvalidCmd("Could not find the mandatory file '%s' in bias module '%s'."%(
662 mandatory_file,run_card['bias_module']))
663- shutil.copytree(run_card['bias_module'], pjoin(self.me_dir,'Source','BIAS',
664+ misc.copytree(run_card['bias_module'], pjoin(self.me_dir,'Source','BIAS',
665 os.path.basename(run_card['bias_module'])))
666
667 #check expected parameters for the module.
668@@ -6572,7 +6575,8 @@
669 self.exec_cmd('systematics %s --from_card' % self.run_name,
670 postcmd=False,printcmd=False)
671 self.exec_cmd('decay_events -from_cards', postcmd=False)
672- elif self.run_card['use_syst']:
673+ elif self.run_card['use_syst'] and self.run_card['systematics_program'] == 'systematics':
674+ self.options['nb_core'] = 1
675 self.exec_cmd('systematics %s --from_card' %
676 pjoin('Events', self.run_name, 'unweighted_events.lhe.gz'),
677 postcmd=False,printcmd=False)
678
679=== modified file 'madgraph/interface/madgraph_interface.py'
680--- madgraph/interface/madgraph_interface.py 2021-05-28 09:43:47 +0000
681+++ madgraph/interface/madgraph_interface.py 2021-08-21 21:19:44 +0000
682@@ -5848,7 +5848,7 @@
683 add_options.remove('--local')
684 logger.warning('you are using a local installer. This is intended for debugging only!')
685 shutil.rmtree(pjoin(MG5DIR,'HEPTools','HEPToolsInstallers'))
686- shutil.copytree(os.path.abspath(pjoin(MG5DIR,os.path.pardir,
687+ misc.copytree(os.path.abspath(pjoin(MG5DIR,os.path.pardir,
688 'HEPToolsInstallers')),pjoin(MG5DIR,'HEPTools','HEPToolsInstallers'))
689
690 # Potential change in naming convention
691@@ -6639,13 +6639,13 @@
692 pattern = re.compile(r'''^=== renamed directory \'(?P<orig>[^\']*)\' => \'(?P<new>[^\']*)\'''')
693 #= = = renamed directory 'Template' => 'Template/LO'
694 for orig, new in pattern.findall(text):
695- shutil.copytree(pjoin(MG5DIR, orig), pjoin(MG5DIR, 'UPDATE_TMP'))
696+ misc.copytree(pjoin(MG5DIR, orig), pjoin(MG5DIR, 'UPDATE_TMP'))
697 full_path = os.path.dirname(pjoin(MG5DIR, new)).split('/')
698 for i, name in enumerate(full_path):
699 path = os.path.sep.join(full_path[:i+1])
700 if path and not os.path.isdir(path):
701 os.mkdir(path)
702- shutil.copytree(pjoin(MG5DIR, 'UPDATE_TMP'), pjoin(MG5DIR, new))
703+ misc.copytree(pjoin(MG5DIR, 'UPDATE_TMP'), pjoin(MG5DIR, new))
704 shutil.rmtree(pjoin(MG5DIR, 'UPDATE_TMP'))
705 # track rename since patch fail to apply those correctly.
706 pattern = re.compile(r'''=== renamed file \'(?P<orig>[^\']*)\' => \'(?P<new>[^\']*)\'''')
707@@ -8555,9 +8555,6 @@
708 logger.info('Pass to numerical integration for computing the widths:')
709 else:
710 logger.info('No need for N body-decay (N>2). Results are in %s' % opts['output'])
711-
712-
713-
714 return decay_info
715
716 # Do the MadEvent integration!!
717@@ -8586,7 +8583,10 @@
718 me_cmd = madevent_interface.MadEventCmd(decay_dir)
719 for name, val in self.options.items():
720 if name in me_cmd.options and me_cmd.options[name] != val:
721- self.exec_cmd('set %s %s --no_save' % (name, val))
722+ try:
723+ me_cmd.exec_cmd('set %s %s --no_save' % (name, val))
724+ except madgraph.InvalidCmd:
725+ continue
726 #me_cmd.options.update(self.options)
727 #me_cmd.configure_run_mode(self.options['run_mode'])
728 #self.define_child_cmd_interface(me_cmd, interface=False)
729
730=== modified file 'madgraph/interface/reweight_interface.py'
731--- madgraph/interface/reweight_interface.py 2021-02-23 08:59:44 +0000
732+++ madgraph/interface/reweight_interface.py 2021-08-21 21:19:44 +0000
733@@ -715,7 +715,6 @@
734 else:
735 rw_dir = pjoin(path_me, 'rw_me')
736
737-
738 if not '--keep_card' in args:
739 if self.has_nlo and self.rwgt_mode != "LO":
740 rwdir_virt = rw_dir.replace('rw_me', 'rw_mevirt')
741@@ -1641,12 +1640,12 @@
742
743 for i in range(2):
744 pdir = pjoin(path_me,data['paths'][i])
745- if os.path.exists(pdir):
746- try:
747- shutil.rmtree(pjoin(path_me,data['paths'][0]))
748- except Exception as error:
749- misc.sprint('fail to rm rwgt dir:', error)
750- pass
751+ if os.path.exists(pdir):
752+ try:
753+ shutil.rmtree(pdir)
754+ except Exception as error:
755+ misc.sprint('fail to rm rwgt dir:', error)
756+ pass
757
758 # 1. prepare the interface----------------------------------------------
759 mgcmd = self.mg5cmd
760@@ -1805,7 +1804,6 @@
761 path_me = self.me_dir
762 else:
763 path_me = self.rwgt_dir
764-
765 self.id_to_path = {}
766 self.id_to_path_second = {}
767 rwgt_dir_possibility = ['rw_me','rw_me_%s' % self.nb_library,'rw_mevirt','rw_mevirt_%s' % self.nb_library]
768@@ -1850,7 +1848,6 @@
769 mymod.set_madloop_path(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources'))
770 if (self.second_model or self.second_process or self.dedicated_path):
771 break
772-
773 data = self.id_to_path
774 if onedir not in ["rw_me", "rw_mevirt"]:
775 data = self.id_to_path_second
776@@ -1861,7 +1858,6 @@
777 all_prefix = [''.join([i.decode() for i in j]).strip().lower() for j in mymod.get_prefix()]
778 prefix_set = set(all_prefix)
779
780-
781 hel_dict={}
782 for prefix in prefix_set:
783 if hasattr(mymod,'%sprocess_nhel' % prefix):
784@@ -1965,6 +1961,7 @@
785 to_save['rwgt_mode'] = self.rwgt_mode
786 to_save['rwgt_name'] = self.options['rwgt_name']
787 to_save['allow_missing_finalstate'] = self.options['allow_missing_finalstate']
788+ to_save['nb_library'] = self.nb_library
789
790 name = pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl')
791 save_load_object.save_to_file(name, to_save)
792@@ -1976,10 +1973,16 @@
793 obj = save_load_object.load_from_file( pjoin(self.rwgt_dir, 'rw_me', 'rwgt.pkl'))
794
795 self.has_standalone_dir = True
796- self.options = {'curr_dir': os.path.realpath(os.getcwd()),
797- 'rwgt_name': None}
798+ if 'rwgt_info' in self.options:
799+ self.options = {'rwgt_info': self.options['rwgt_info']}
800+ else:
801+ self.options = {}
802+ self.options.update({'curr_dir': os.path.realpath(os.getcwd()),
803+ 'rwgt_name': None})
804+
805 if keep_name:
806 self.options['rwgt_name'] = obj['rwgt_name']
807+
808 self.options['allow_missing_finalstate'] = obj['allow_missing_finalstate']
809 old_rwgt = obj['rwgt_dir']
810
811@@ -2000,6 +2003,7 @@
812 self.second_process = obj['second_process']
813 self.second_model = obj['second_model']
814 self.has_nlo = obj['has_nlo']
815+ self.nb_library = obj['nb_library']
816 if not self.rwgt_mode:
817 self.rwgt_mode = obj['rwgt_mode']
818 logger.info("mode set to %s" % self.rwgt_mode)
819
820=== modified file 'madgraph/iolibs/drawing_eps.py'
821--- madgraph/iolibs/drawing_eps.py 2020-02-12 08:38:24 +0000
822+++ madgraph/iolibs/drawing_eps.py 2021-08-21 21:19:44 +0000
823@@ -284,7 +284,7 @@
824
825 #add the code in the correct format
826 x, y = self.rescale(line.begin.pos_x, line.begin.pos_y)
827- self.text += '%s %s moveto'%(x, y)
828+ self.text += '%s %s moveto \n'%(x, y)
829 self.text += self.line_format(line.begin.pos_x, line.begin.pos_y,
830 line.end.pos_x+0.01*direction[0], line.end.pos_y+0.01*direction[1],
831 '%s Fhiggsl' % (curvature*7))
832
833=== modified file 'madgraph/iolibs/export_fks.py'
834--- madgraph/iolibs/export_fks.py 2021-02-14 22:27:00 +0000
835+++ madgraph/iolibs/export_fks.py 2021-08-21 21:19:44 +0000
836@@ -98,7 +98,7 @@
837 raise MadGraph5Error("No valid MG_ME path given for MG4 run directory creation.")
838 logger.info('initialize a new directory: %s' % \
839 os.path.basename(dir_path))
840- shutil.copytree(os.path.join(mgme_dir, 'Template', 'NLO'), dir_path, True)
841+ misc.copytree(os.path.join(mgme_dir, 'Template', 'NLO'), dir_path, True)
842 # misc.copytree since dir_path already exists
843 misc.copytree(pjoin(self.mgme_dir, 'Template', 'Common'),dir_path)
844 # Copy plot_card
845@@ -307,7 +307,7 @@
846 except OSError as error:
847 pass
848 model_path = model.get('modelpath')
849- shutil.copytree(model_path,
850+ misc.copytree(model_path,
851 pjoin(self.dir_path,'bin','internal','ufomodel'),
852 ignore=shutil.ignore_patterns(*IGNORE_PATTERNS))
853 if hasattr(model, 'restrict_card'):
854@@ -881,7 +881,7 @@
855
856 elif output_dependencies == 'internal':
857 StdHEP_internal_path = pjoin(self.dir_path,'Source','StdHEP')
858- shutil.copytree(StdHep_path,StdHEP_internal_path, symlinks=True)
859+ misc.copytree(StdHep_path,StdHEP_internal_path, symlinks=True)
860 # Create the links to the lib folder
861 linkfiles = ['libstdhep.a', 'libFmcfio.a']
862 for file in linkfiles:
863@@ -3412,7 +3412,7 @@
864 raise MadGraph5Error("No valid MG_ME path given for MG4 run directory creation.")
865 logger.info('initialize a new directory: %s' % \
866 os.path.basename(dir_path))
867- shutil.copytree(os.path.join(mgme_dir, 'Template', 'NLO'), dir_path, True)
868+ misc.copytree(os.path.join(mgme_dir, 'Template', 'NLO'), dir_path, True)
869 # misc.copytree since dir_path already exists
870 misc.copytree(pjoin(self.mgme_dir, 'Template', 'Common'),
871 dir_path)
872
873=== modified file 'madgraph/iolibs/export_v4.py'
874--- madgraph/iolibs/export_v4.py 2021-05-20 21:17:50 +0000
875+++ madgraph/iolibs/export_v4.py 2021-08-21 21:19:44 +0000
876@@ -255,7 +255,7 @@
877 "No valid MG_ME path given for MG4 run directory creation."
878 logger.info('initialize a new directory: %s' % \
879 os.path.basename(self.dir_path))
880- shutil.copytree(pjoin(self.mgme_dir, 'Template/LO'),
881+ misc.copytree(pjoin(self.mgme_dir, 'Template/LO'),
882 self.dir_path, True)
883 # misc.copytree since dir_path already exists
884 misc.copytree(pjoin(self.mgme_dir, 'Template/Common'),
885@@ -280,7 +280,7 @@
886 # if os.path.isfile(filename):
887 # files.cp(filename, pjoin(self.dir_path,name))
888 # elif os.path.isdir(filename):
889-# shutil.copytree(filename, pjoin(self.dir_path,name), True)
890+# misc.copytree(filename, pjoin(self.dir_path,name), True)
891 # misc.copytree since dir_path already exists
892 misc.copytree(pjoin(self.mgme_dir, 'Template/Common'),
893 self.dir_path)
894@@ -2950,9 +2950,9 @@
895 super(ProcessExporterFortranMW, self).copy_template(model)
896
897 # Add the MW specific file
898- shutil.copytree(pjoin(MG5DIR,'Template','MadWeight'),
899+ misc.copytree(pjoin(MG5DIR,'Template','MadWeight'),
900 pjoin(self.dir_path, 'Source','MadWeight'), True)
901- shutil.copytree(pjoin(MG5DIR,'madgraph','madweight'),
902+ misc.copytree(pjoin(MG5DIR,'madgraph','madweight'),
903 pjoin(self.dir_path, 'bin','internal','madweight'), True)
904 files.mv(pjoin(self.dir_path, 'Source','MadWeight','src','setrun.f'),
905 pjoin(self.dir_path, 'Source','setrun.f'))
906@@ -2999,7 +2999,7 @@
907 pass
908 model_path = model.get('modelpath')
909 # This is not safe if there is a '##' or '-' in the path.
910- shutil.copytree(model_path,
911+ misc.copytree(model_path,
912 pjoin(self.dir_path,'bin','internal','ufomodel'),
913 ignore=shutil.ignore_patterns(*IGNORE_PATTERNS))
914 if hasattr(model, 'restrict_card'):
915@@ -3846,7 +3846,7 @@
916 pass
917 model_path = model.get('modelpath')
918 # This is not safe if there is a '##' or '-' in the path.
919- shutil.copytree(model_path,
920+ misc.copytree(model_path,
921 pjoin(self.dir_path,'bin','internal','ufomodel'),
922 ignore=shutil.ignore_patterns(*IGNORE_PATTERNS))
923 if hasattr(model, 'restrict_card'):
924
925=== modified file 'madgraph/iolibs/file_writers.py'
926--- madgraph/iolibs/file_writers.py 2020-12-06 20:55:21 +0000
927+++ madgraph/iolibs/file_writers.py 2021-08-21 21:19:44 +0000
928@@ -456,7 +456,7 @@
929 else:
930 if not line.endswith('\n'):
931 line = '%s\n' % line
932- file.writelines(self, line)
933+ super(FileWriter,self).writelines(line)
934 else:
935 removed.append(line)
936
937
938=== modified file 'madgraph/iolibs/files.py'
939--- madgraph/iolibs/files.py 2019-07-23 09:56:18 +0000
940+++ madgraph/iolibs/files.py 2021-08-21 21:19:44 +0000
941@@ -152,7 +152,7 @@
942 try:
943 if os.path.exists(path2):
944 path2 = os.path.join(path2, os.path.split(path1)[1])
945- shutil.copytree(path1, path2)
946+ misc.copytree(path1, path2)
947 except IOError as why:
948 if error:
949 raise
950
951=== modified file 'madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc'
952--- madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc 2021-02-09 14:08:25 +0000
953+++ madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc 2021-08-21 21:19:44 +0000
954@@ -407,7 +407,7 @@
955 IMPLICIT REAL*8(A-H,O-Z)
956 INTEGER NEXTERNAL, NINCOMING
957 PARAMETER (NEXTERNAL=%(nexternal)d,NINCOMING=%(nincoming)d)
958- DIMENSION XM(NEXTERNAL-NINCOMING),P(4,NEXTERNAL-NINCOMING)
959+ DIMENSION XM(*),P(4,*)
960 DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4),B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING),E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5)
961 SAVE ACC,ITMAX,IBEGIN,IWARN
962 DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
963
964=== modified file 'madgraph/loop/loop_exporters.py'
965--- madgraph/loop/loop_exporters.py 2021-05-20 21:17:50 +0000
966+++ madgraph/loop/loop_exporters.py 2021-08-21 21:19:44 +0000
967@@ -136,7 +136,7 @@
968
969 if self.dependencies=='internal':
970 new_CT_path = pjoin(targetPath,'Source','CutTools')
971- shutil.copytree(self.cuttools_dir, new_CT_path, symlinks=True)
972+ misc.copytree(self.cuttools_dir, new_CT_path, symlinks=True)
973
974 current = misc.detect_current_compiler(os.path.join(new_CT_path,
975 'makefile'))
976@@ -2029,7 +2029,7 @@
977 elif tir_name == "iregi":
978 # This is the right paths for IREGI
979 new_iregi_path = pjoin(targetPath,os.path.pardir,'Source','IREGI')
980- shutil.copytree(pjoin(libpath,os.path.pardir), new_iregi_path,
981+ misc.copytree(pjoin(libpath,os.path.pardir), new_iregi_path,
982 symlinks=True)
983
984 current = misc.detect_current_compiler(
985
986=== modified file 'madgraph/madweight/write_MadWeight.py'
987--- madgraph/madweight/write_MadWeight.py 2020-03-26 22:18:37 +0000
988+++ madgraph/madweight/write_MadWeight.py 2021-08-21 21:19:44 +0000
989@@ -1001,7 +1001,7 @@
990 if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
991 call UPDATE_AS_PARAM()
992 endif
993- if(.not.fixed_fac_scale) then
994+ if(.not.fixed_fac_scale1.or..not.fixed_fac_scale2) then
995 call set_fac_scale(momenta(0,1),q2fact)
996 endif
997 """
998
999=== modified file 'madgraph/various/banner.py'
1000--- madgraph/various/banner.py 2021-03-17 17:12:05 +0000
1001+++ madgraph/various/banner.py 2021-08-21 21:19:44 +0000
1002@@ -2995,7 +2995,9 @@
1003 self.add_param("pdlabel", "nn23lo1", allowed=['lhapdf', 'cteq6_m','cteq6_l', 'cteq6l1','nn23lo', 'nn23lo1', 'nn23nlo']),
1004 self.add_param("lhaid", 230000, hidden=True)
1005 self.add_param("fixed_ren_scale", False)
1006- self.add_param("fixed_fac_scale", False)
1007+ self.add_param("fixed_fac_scale", False, hidden=True, include=False, comment="define if the factorization scale is fixed or not. You can define instead fixed_fac_scale1 and fixed_fac_scale2 if you want to make that choice per beam")
1008+ self.add_param("fixed_fac_scale1", False, hidden=True)
1009+ self.add_param("fixed_fac_scale2", False, hidden=True)
1010 self.add_param("scale", 91.1880)
1011 self.add_param("dsqrt_q2fact1", 91.1880, fortran_name="sf1")
1012 self.add_param("dsqrt_q2fact2", 91.1880, fortran_name="sf2")
1013@@ -3300,10 +3302,51 @@
1014 if abs(self['lpp%s' % i ]) == 2 and self['dsqrt_q2fact%s'%i] == 91.188:
1015 logger.warning("Since 2.7.1 Photon from proton are using fixed scale value of muf [dsqrt_q2fact%s] as the cut of the Improved Weizsaecker-Williams formula. Please edit it accordingly." % i)
1016 time.sleep(5)
1017+
1018+ # check that fixed_fac_scale(1/2) is setting as expected
1019+ # if lpp=2/3/4 -> default is that beam in fixed scale
1020+ # check that fixed_fac_scale is not setup if fixed_fac_scale1/2 are
1021+ # check that both fixed_fac_scale1/2 are defined together
1022+ # ensure that fixed_fac_scale1 and fixed_fac_scale2 are setup as needed
1023+ if 'fixed_fac_scale1' in self.user_set:
1024+ if 'fixed_fac_scale2' in self.user_set:
1025+ if 'fixed_fac_scale' in self.user_set:
1026+ if not (self['fixed_fac_scale'] == self['fixed_fac_scale2'] == self['fixed_fac_scale2']):
1027+ logger.warning('Both fixed_fac_scale, fixed_fac_scale1 and fixed_fac_scale2 are defined. The value of fixed_fac_scale is ignored')
1028+ elif 'fixed_fac_scale' in self.user_set:
1029+ logger.warning('Both fixed_fac_scale, fixed_fac_scale1 are defined but not fixed_fac_scale2. The value of fixed_fac_scale2 will be set to the one of fixed_fac_scale')
1030+ self['fixed_fac_scale2'] = self['fixed_fac_scale']
1031+ elif self['lpp2'] !=0:
1032+ raise Exception('fixed_fac_scale2 not defined while fixed_fac_scale1 is. Please fix your run_card.')
1033+ elif 'fixed_fac_scale2' in self.user_set:
1034+ if 'fixed_fac_scale' in self.user_set:
1035+ logger.warning('Both fixed_fac_scale, fixed_fac_scale2 are defined but not fixed_fac_scale1. The value of fixed_fac_scale1 will be set to the one of fixed_fac_scale')
1036+ self['fixed_fac_scale1'] = self['fixed_fac_scale']
1037+ elif self['lpp1'] !=0:
1038+ raise Exception('fixed_fac_scale1 not defined while fixed_fac_scale2 is. Please fix your run_card.')
1039+ else:
1040+ if 'fixed_fac_scale' in self.user_set:
1041+ if self['lpp1'] in [2,3,4]:
1042+ logger.warning('fixed factorization scale is used for beam1. You can prevent this by setting fixed_fac_scale1 to False')
1043+ self['fixed_fac_scale1'] = True
1044+ self['fixed_fac_scale2'] = self['fixed_fac_scale']
1045+ elif self['lpp2'] in [2,3,4]:
1046+ logger.warning('fixed factorization scale is used for beam2. You can prevent this by setting fixed_fac_scale2 to False')
1047+ self['fixed_fac_scale1'] = self['fixed_fac_scale']
1048+ self['fixed_fac_scale2'] = True
1049+ else:
1050+ self['fixed_fac_scale1'] = self['fixed_fac_scale']
1051+ self['fixed_fac_scale2'] = self['fixed_fac_scale']
1052+ elif self['lpp1'] !=0 or self['lpp2']!=0:
1053+ raise Exception('fixed_fac_scale not defined whithin your run_card. Plase fix this.')
1054+
1055+
1056+
1057
1058 # if both lpp1/2 are on PA mode -> force fixed factorization scale
1059 if abs(self['lpp1']) in [2, 3,4] and abs(self['lpp2']) in [2, 3,4] and not self['fixed_fac_scale']:
1060- raise InvalidRunCard("Having both beam in elastic photon mode requires fixed_fac_scale to be on True [since this is use as cutoff]")
1061+ if 'fixed_fac_scale1' not in self.user_set or 'fixed_fac_scale2' not in self.user_set:
1062+ raise InvalidRunCard("Having both beam in elastic photon mode requires fixed_fac_scale to be on True [since this is use as cutoff]. If you really want a running scale here, please define fixed_fac_scale1 on False and fixed_fac_scale2 on False")
1063
1064 if six.PY2 and self['hel_recycling']:
1065 self['hel_recycling'] = False
1066
1067=== modified file 'madgraph/various/misc.py'
1068--- madgraph/various/misc.py 2021-03-05 20:18:16 +0000
1069+++ madgraph/various/misc.py 2021-08-21 21:19:44 +0000
1070@@ -461,6 +461,12 @@
1071 name = "%s[%s-%s]%s" % (base, first[len(base):], last[len(base):],end)
1072 return name
1073
1074+def copytree(*args, **opts):
1075+
1076+ if 'copy_function' not in opts:
1077+ opts['copy_function'] = shutil.copy
1078+ return misc.copytree(*args, **opts)
1079+
1080 #===============================================================================
1081 # Compiler which returns smart output error in case of trouble
1082 #===============================================================================
1083
1084=== modified file 'madgraph/various/process_checks.py'
1085--- madgraph/various/process_checks.py 2020-11-25 20:43:03 +0000
1086+++ madgraph/various/process_checks.py 2021-08-21 21:19:44 +0000
1087@@ -400,7 +400,8 @@
1088 if events:
1089 ids = [l.get('id') for l in sorted_legs]
1090 import MadSpin.decay as madspin
1091- if not hasattr(self, 'event_file'):
1092+ if not hasattr(self, 'event_file') or self.event_file.inputfile.closed:
1093+ print( "reset")
1094 fsock = open(events)
1095 self.event_file = madspin.Event(fsock)
1096
1097@@ -419,6 +420,7 @@
1098 for part in event.values():
1099 m = part['momentum']
1100 p.append([m.E, m.px, m.py, m.pz])
1101+ fsock.close()
1102 return p, 1
1103
1104 nincoming = len([leg for leg in sorted_legs if leg.get('state') == False])
1105
1106=== modified file 'mg5decay/decay_objects.py'
1107--- mg5decay/decay_objects.py 2020-09-16 07:01:30 +0000
1108+++ mg5decay/decay_objects.py 2021-08-21 21:19:44 +0000
1109@@ -3418,7 +3418,9 @@
1110 numbers = [l.get('number') for l in self['final_legs'] if l.get('id') == handling]
1111
1112 for new_numbers in itertools.permutations(numbers):
1113+
1114 mapping_id = dict([(o,n) for o,n in zip(numbers, new_numbers) if o!=n])
1115+
1116 if not mapping_id:
1117 out.append(self)
1118 continue
1119@@ -3432,6 +3434,7 @@
1120 channel['helas_number'] = None
1121 # diagram written by IdentifyHelasTag
1122 channel['std_diagram'] = None
1123+
1124 for l,vertex in enumerate(self['vertices']):
1125 new_vertex = copy.copy(vertex)
1126 new_vertex['legs'] = base_objects.LegList()
1127@@ -3442,16 +3445,16 @@
1128 new_leg.set('number', mapping_id[leg['number']])
1129 new_vertex['legs'].append(new_leg)
1130 else:
1131- new_vertex['legs'].append(leg)
1132- min_id = min(min_id, leg['number'])
1133-
1134- if min_id != new_vertex['legs'][-1]['number']:
1135- if l != len(self['vertices']) -1:
1136- mapping_id[new_vertex['legs'][-1]['number']] = min_id
1137- new_vertex['legs'][-1]['number'] = min_id
1138+ new_leg = copy.copy(leg)
1139+ new_vertex['legs'].append(new_leg)
1140+ min_id = min(min_id, new_leg['number'])
1141+
1142+ if l != len(self['vertices']) -1:
1143+ new_vertex['legs'][-1]['number'] = min_id
1144+ mapping_id[vertex['legs'][-1]['number']] = min_id
1145+
1146 channel['vertices'].append(new_vertex)
1147- out.append(channel)
1148-
1149+ out.append(channel)
1150
1151 # do the recursion
1152 if len(remain_id) > 1:
1153
1154=== modified file 'models/build_restriction_lib.py'
1155--- models/build_restriction_lib.py 2015-10-01 16:00:08 +0000
1156+++ models/build_restriction_lib.py 2021-08-21 21:19:44 +0000
1157@@ -50,4 +50,4 @@
1158
1159
1160
1161-
1162\ No newline at end of file
1163+
1164
1165=== modified file 'models/usermod.py'
1166--- models/usermod.py 2019-07-27 21:18:13 +0000
1167+++ models/usermod.py 2021-08-21 21:19:44 +0000
1168@@ -588,6 +588,7 @@
1169 logger.warning('The particle name \'%s\' is present in both model with different pdg code' % name)
1170 logger.warning('The particle coming from the plug-in model will be rename to \'%s%s\'' % (name, self.addon))
1171 particle.name = '%s%s' % (name, self.addon)
1172+ particle.antiname = '%s%s' % (particle.antiname, self.addon)
1173 self.particles.append(particle)
1174 return
1175 elif identify:
1176
1177=== modified file 'tests/__init__.py'
1178--- tests/__init__.py 2020-01-22 15:16:50 +0000
1179+++ tests/__init__.py 2021-08-21 21:19:44 +0000
1180@@ -1,1 +1,1 @@
1181-to_bypass = []
1182\ No newline at end of file
1183+to_bypass = []
1184
1185=== modified file 'tests/acceptance_tests/test_cmd_madevent.py'
1186--- tests/acceptance_tests/test_cmd_madevent.py 2021-01-21 13:01:31 +0000
1187+++ tests/acceptance_tests/test_cmd_madevent.py 2021-08-21 21:19:44 +0000
1188@@ -188,11 +188,11 @@
1189 if self.debugging:
1190 if os.path.isdir(pjoin(MG5DIR,'BackUp_tmp_test')):
1191 shutil.rmtree(pjoin(MG5DIR,'BackUp_tmp_test'))
1192- shutil.copytree(pjoin(MG5DIR,'tmp_test'),
1193+ misc.copytree(pjoin(MG5DIR,'tmp_test'),
1194 pjoin(MG5DIR,'BackUp_tmp_test'))
1195 else:
1196 shutil.rmtree(pjoin(MG5DIR,'tmp_test'))
1197- shutil.copytree(pjoin(MG5DIR,'BackUp_tmp_test'),pjoin(MG5DIR,'tmp_test'))
1198+ misc.copytree(pjoin(MG5DIR,'BackUp_tmp_test'),pjoin(MG5DIR,'tmp_test'))
1199
1200 biased_events = lhe_parser.EventFile(pjoin(self.out_dir, 'Events','run_01','unweighted_events.lhe.gz'))
1201 unbiased_events = lhe_parser.EventFile(pjoin(self.out_dir, 'Events','run_02','unweighted_events.lhe.gz'))
1202
1203=== modified file 'tests/acceptance_tests/test_madspin.py'
1204--- tests/acceptance_tests/test_madspin.py 2019-04-17 14:39:47 +0000
1205+++ tests/acceptance_tests/test_madspin.py 2021-08-21 21:19:44 +0000
1206@@ -164,4 +164,4 @@
1207 muon_in +=1
1208 self.assertEqual(muon_in, 1)
1209 self.assertEqual(nb_dec, 189)
1210- self.assertEqual(nb_muon, 100)
1211\ No newline at end of file
1212+ self.assertEqual(nb_muon, 100)
1213
1214=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f'
1215--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f 2021-03-25 07:55:31 +0000
1216+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f 2021-08-21 21:19:44 +0000
1217@@ -193,7 +193,7 @@
1218 R=XRAN1(IDUM)*ANS
1219 SUMHEL=0D0
1220 DO I=1,NCOMB
1221- SUMHEL=SUMHEL+TS(I)
1222+ SUMHEL=SUMHEL+DABS(TS(I))
1223 IF(R.LT.SUMHEL)THEN
1224 WRITE(HEL_BUFF,'(20i5)')(NHEL(II,I),II=1,NEXTERNAL)
1225 ANS=DSIGN(ANS,TS(I))
1226
1227=== modified file 'tests/input_files/IOTestsComparison/TestMadWeight/modification_to_cuts/cuts.f'
1228--- tests/input_files/IOTestsComparison/TestMadWeight/modification_to_cuts/cuts.f 2020-05-13 10:56:30 +0000
1229+++ tests/input_files/IOTestsComparison/TestMadWeight/modification_to_cuts/cuts.f 2021-08-21 21:19:44 +0000
1230@@ -1,1742 +1,1726 @@
1231- logical function pass_point(p)
1232-c************************************************************************
1233-c This function is called from sample to see if it needs to
1234-c bother calculating the weight from all the different conficurations
1235-c You can either just return true, or have it call passcuts
1236-c************************************************************************
1237- implicit none
1238-c
1239-c Arguments
1240-c
1241- double precision p
1242-c
1243-c External
1244-c
1245- logical passcuts
1246- external passcuts
1247-c-----
1248-c Begin Code
1249-c-----
1250- pass_point = .true.
1251-c pass_point = passcuts(p)
1252- end
1253-C
1254- LOGICAL FUNCTION PASSCUTS(P)
1255-C**************************************************************************
1256-C INPUT:
1257-C P(0:3,1) MOMENTUM OF INCOMING PARTON
1258-C P(0:3,2) MOMENTUM OF INCOMING PARTON
1259-C P(0:3,3) MOMENTUM OF ...
1260-C ALL MOMENTA ARE IN THE REST FRAME!!
1261-C COMMON/JETCUTS/ CUTS ON JETS
1262-C OUTPUT:
1263-C TRUE IF EVENTS PASSES ALL CUTS LISTED
1264-C**************************************************************************
1265- IMPLICIT NONE
1266-c
1267-c Constants
1268-c
1269- include 'maxparticles.inc'
1270- include 'nexternal.inc'
1271-C
1272-C ARGUMENTS
1273-C
1274- REAL*8 P(0:3,nexternal)
1275-
1276-C
1277-C LOCAL
1278-C
1279- LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
1280- LOGICAL DEBUG
1281- integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
1282- REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
1283- real*8 ptemp(0:3), ptemp2(0:3)
1284- character*20 formstr
1285-C
1286-C PARAMETERS
1287-C
1288- real*8 PI
1289- parameter( PI = 3.14159265358979323846d0 )
1290-C
1291-C EXTERNAL
1292-C
1293- REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
1294- logical cut_bw,setclscales,dummy_cuts
1295- external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
1296- external dummy_cuts
1297-C
1298-C GLOBAL
1299-C
1300- include 'run.inc'
1301- include 'cuts.inc'
1302-
1303- double precision ptjet(nexternal)
1304- double precision ptheavyjet(nexternal)
1305- double precision ptlepton(nexternal)
1306- double precision temp
1307-
1308-C VARIABLES TO SPECIFY JETS
1309- DOUBLE PRECISION PJET(NEXTERNAL,0:3)
1310- DOUBLE PRECISION PTMIN
1311- DOUBLE PRECISION PT1,PT2
1312-
1313-C INTEGERS FOR COUNTING.
1314- INTEGER K,L,J1,J2
1315-
1316-C VARIABLES FOR KT CUT
1317- DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
1318- DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
1319-
1320- double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
1321- double precision emin(nincoming+1:nexternal)
1322- double precision r2min(nincoming+1:nexternal,nincoming+1:nexternal)
1323- double precision s_min(nexternal,nexternal)
1324- double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
1325- double precision emax(nincoming+1:nexternal)
1326- double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
1327- double precision s_max(nexternal,nexternal)
1328- double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
1329- double precision inclHtmin,inclHtmax
1330- common/to_cuts/ etmin, emin, etamax, r2min, s_min,
1331- $ etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
1332-
1333- double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
1334- logical jetor
1335- common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
1336-
1337- double precision ptlmin4(4),ptlmax4(4)
1338- common/to_lepton_cuts/ ptlmin4,ptlmax4
1339-
1340-c
1341-c Special cuts
1342-c
1343-
1344- integer lbw(0:nexternal) !Use of B.W.
1345- common /to_BW/ lbw
1346-C
1347-C SPECIAL CUTS
1348-C
1349- LOGICAL IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
1350- LOGICAL IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
1351- LOGICAL IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
1352- logical do_cuts(nexternal)
1353- COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
1354- . IS_A_ONIUM, do_cuts
1355-
1356-C
1357-C MERGING SCALE CUT
1358-C
1359-C Retrieve which external particles undergo the ktdurham and ptlund cuts.
1360- LOGICAL is_pdg_for_merging_cut(NEXTERNAL)
1361- logical from_decay(-(nexternal+3):nexternal)
1362- COMMON /TO_MERGE_CUTS/is_pdg_for_merging_cut, from_decay
1363-
1364-C
1365-C ADDITIONAL VARIABLES FOR PTLUND CUT
1366-C
1367- INTEGER NMASSLESS
1368- DOUBLE PRECISION PINC(NINCOMING,0:3)
1369- DOUBLE PRECISION PRADTEMP(0:3), PRECTEMP(0:3), PEMTTEMP(0:3)
1370- DOUBLE PRECISION PTMINSAVE, RHOPYTHIA
1371- EXTERNAL RHOPYTHIA
1372-C
1373-C FLAVOUR INFORMATION NECESSARY TO RECONSTRUCT PTLUND
1374-C
1375- INTEGER JETFLAVOUR(NEXTERNAL), INCFLAVOUR(NINCOMING)
1376- include 'maxamps.inc'
1377- integer idup(nexternal,maxproc,maxsproc)
1378- integer mothup(2,nexternal)
1379- integer icolup(2,nexternal,maxflow,maxsproc)
1380- include 'leshouche.inc'
1381-
1382-C
1383-C Keep track of whether cuts already calculated for this event
1384-C
1385- LOGICAL CUTSDONE,CUTSPASSED
1386- COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
1387- DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
1388-
1389-C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
1390-
1391- double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
1392- common/to_xqcuts/xqcutij,xqcuti
1393-
1394-c jet cluster algorithm
1395- integer nQCD !,NJET,JET(nexternal)
1396-c double precision plab(0:3, nexternal)
1397- double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
1398-c double precision rfj,sycut,palg,fastjetdmerge
1399-c integer njet_eta
1400-c Photon isolation
1401- integer nph,nem,nin
1402- double precision ptg,chi_gamma_iso,iso_getdrv40
1403- double precision Etsum(0:nexternal)
1404- real drlist(nexternal)
1405- double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
1406- logical alliso
1407-C Sort array of results: ismode>0 for real, isway=0 for ascending order
1408- integer ismode,isway,izero,isorted(nexternal)
1409- parameter (ismode=1)
1410- parameter (isway=0)
1411- parameter (izero=0)
1412-
1413- include 'coupl.inc'
1414-
1415-C
1416-C
1417-c
1418- DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
1419-
1420-c put momenta in common block for couplings.f
1421- double precision pp(0:3,max_particles)
1422- common /momenta_pp/pp
1423-
1424- DATA DEBUG/.FALSE./
1425-
1426-C-----
1427-C BEGIN CODE
1428-C-----
1429-
1430-
1431-
1432- PASSCUTS=.TRUE. !EVENT IS OK UNLESS OTHERWISE CHANGED
1433- IF (FIRSTTIME) THEN
1434- FIRSTTIME=.FALSE.
1435-c Preparation for reweighting by setting up clustering by diagrams
1436- ! Remove for MW!call initcluster()
1437-c
1438-c
1439- write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
1440- write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
1441- write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
1442- write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
1443- write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
1444- write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
1445- write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
1446- write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
1447- do j=nincoming+1,nexternal-1
1448- write(*,formstr) 'd R #',j,' >',(-0.0,i=nincoming+1,j),
1449- & (r2min(i,j),i=j+1,nexternal)
1450- do i=j+1,nexternal
1451- r2min(i,j)=r2min(i,j)*dabs(r2min(i,j)) !Since r2 returns distance squared
1452- r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
1453- enddo
1454- enddo
1455- do j=nincoming+1,nexternal-1
1456- write(*,formstr) 's min #',j,'>',
1457- & (s_min(i,j),i=nincoming+1,nexternal)
1458- enddo
1459- do j=nincoming+1,nexternal-1
1460- write(*,formstr) 'xqcutij #',j,'>',
1461- & (xqcutij(i,j),i=nincoming+1,nexternal)
1462- enddo
1463-
1464-cc
1465-cc Set the strong coupling
1466-cc
1467-c call set_ren_scale(P,scale)
1468-c
1469-cc Check that the user funtions for setting the scales
1470-cc have been edited if the choice of an event-by-event
1471-cc scale choice has been made
1472-c
1473-c if(.not.fixed_ren_scale) then
1474-c if(scale.eq.0d0) then
1475-c write(6,*)
1476-c write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
1477-c write(6,*) ' Dynamical renormalization scale choice '
1478-c write(6,*) ' selected but user subroutine'
1479-c write(6,*) ' set_ren_scale not edited in file:setpara.f'
1480-c write(6,*) ' Switching to a fixed_ren_scale choice'
1481-c write(6,*) ' with scale=zmass'
1482-c scale=91.2d0
1483-c write(6,*) 'scale=',scale
1484-c fixed_ren_scale=.true.
1485-c call set_ren_scale(P,scale)
1486-c endif
1487-c endif
1488-
1489-c if(.not.fixed_fac_scale) then
1490-c call set_fac_scale(P,q2fact)
1491-c if(q2fact(1).eq.0d0.or.q2fact(2).eq.0d0) then
1492-c write(6,*)
1493-c write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
1494-c write(6,*) ' Dynamical renormalization scale choice '
1495-c write(6,*) ' selected but user subroutine'
1496-c write(6,*) ' set_fac_scale not edited in file:setpara.f'
1497-c write(6,*) ' Switching to a fixed_fac_scale choice'
1498-c write(6,*) ' with q2fact(i)=zmass**2'
1499-c fixed_fac_scale=.true.
1500-c q2fact(1)=91.2d0**2
1501-c q2fact(2)=91.2d0**2
1502-c write(6,*) 'scales=',q2fact(1),q2fact(2)
1503-c endif
1504-c endif
1505-
1506- if(fixed_ren_scale) then
1507- G = SQRT(4d0*PI*ALPHAS(scale))
1508- call update_as_param()
1509- endif
1510-
1511-c Put momenta in the common block to zero to start
1512- do i=0,3
1513- do j=1,max_particles
1514- pp(i,j) = 0d0
1515- enddo
1516- enddo
1517-
1518- ENDIF ! IF FIRSTTIME
1519-
1520- IF (CUTSDONE) THEN
1521- PASSCUTS=CUTSPASSED
1522- RETURN
1523- ENDIF
1524- CUTSDONE=.TRUE.
1525-c CUTSPASSED=.FALSE.
1526-
1527-c
1528-c Make sure have reasonable 4-momenta
1529-c
1530- if (p(0,1) .le. 0d0) then
1531- passcuts=.false.
1532- return
1533- endif
1534-
1535-c Also make sure there's no INF or NAN
1536- do i=1,nexternal
1537- do j=0,3
1538- if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
1539- passcuts=.false.
1540- return
1541- endif
1542- enddo
1543- enddo
1544-
1545-c
1546-c Limit S_hat
1547-c
1548-c if (x1*x2*stot .gt. 500**2) then
1549-c passcuts=.false.
1550-c return
1551-c endif
1552-
1553-C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
1554-
1555- if(debug) write (*,*) '============================='
1556- if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED '
1557- if(debug) write (*,*) '============================='
1558-c
1559-c p_t min & max cuts
1560-c
1561- do i=nincoming+1,nexternal
1562- if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),' ',etmin(i),
1563- $ ':',etmax(i)
1564- notgood=(pt(p(0,i)) .lt. etmin(i)).or.
1565- & (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
1566- if (notgood) then
1567- if(debug) write (*,*) i,' -> fails'
1568- passcuts=.false.
1569- return
1570- endif
1571- enddo
1572-c
1573-c missing ET min & max cut + Invariant mass of leptons and neutrino
1574-c nb: missing Et defined as the vector sum over the neutrino's pt
1575-c
1576-c-- reset ptemp(0:3)
1577- do j=0,3
1578- ptemp(j)=0 ! for the neutrino
1579- ptemp2(j)=0 ! for the leptons
1580- enddo
1581-c- sum over the momenta
1582- do i=nincoming+1,nexternal
1583- if(is_a_nu(i)) then
1584- if(debug) write (*,*) i,' -> neutrino '
1585- do j=0,3
1586- ptemp(j)=ptemp(j)+p(j,i)
1587- enddo
1588- elseif(is_a_l(i)) then
1589- if(debug) write (*,*) i,' -> lepton '
1590- do j=0,3
1591- ptemp2(j)=ptemp2(j)+p(j,i)
1592- enddo
1593- endif
1594-
1595- enddo
1596-c- check the et
1597- if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
1598- if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),' ',misset,':',missetmax
1599- if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
1600- if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
1601- if(ptemp(0).gt.0d0) then
1602- notgood=(pt(ptemp(0)) .lt. misset).or.
1603- & (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
1604- if (notgood) then
1605- if(debug) write (*,*) ' missing et cut -> fails'
1606- passcuts=.false.
1607- return
1608- endif
1609- endif
1610- if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
1611- if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
1612- if(debug) write (*,*) 'lepton invariant mass -> fails'
1613- passcuts=.false.
1614- return
1615- endif
1616- endif
1617-c
1618-c pt cut on heavy particles
1619-c gives min(pt) for (at least) one heavy particle
1620-c
1621- if(ptheavy.gt.0d0)then
1622- passcuts=.false.
1623- foundheavy=.false.
1624- do i=nincoming+1,nexternal
1625- if(is_heavy(i)) then
1626- if(debug) write (*,*) i,' -> heavy '
1627- foundheavy=.true.
1628- if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
1629- endif
1630- enddo
1631-
1632- if(.not.passcuts.and.foundheavy)then
1633- if(debug) write (*,*) ' heavy particle cut -> fails'
1634- return
1635- else
1636- passcuts=.true.
1637- endif
1638- endif
1639-c
1640-c E min & max cuts
1641-c
1642- do i=nincoming+1,nexternal
1643- if(debug) write (*,*) 'p(0,',i,')=',p(0,i),' ',emin(i),':',emax(i)
1644- notgood=(p(0,i) .le. emin(i)).or.
1645- & (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
1646- if (notgood) then
1647- if(debug) write (*,*) i,' -> fails'
1648- passcuts=.false.
1649- return
1650- endif
1651- enddo
1652-c
1653-c Rapidity min & max cuts
1654-c
1655- do i=nincoming+1,nexternal
1656- if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),' ',etamin(i),':',etamax(i)
1657- notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
1658- & (abs(rap(p(0,i))) .lt. etamin(i))
1659- if (notgood) then
1660- if(debug) write (*,*) i,' -> fails'
1661- passcuts=.false.
1662- return
1663- endif
1664- enddo
1665-c
1666-c DeltaR min & max cuts
1667-c
1668- do i=nincoming+1,nexternal-1
1669- do j=i+1,nexternal
1670- if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
1671- if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
1672- if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
1673- tmp=r2(p(0,i),p(0,j))
1674- notgood=(tmp .lt. r2min(j,i)).or.
1675- $ (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
1676- if (notgood) then
1677- if(debug) write (*,*) i,j,' -> fails'
1678- passcuts=.false.
1679- return
1680- endif
1681- endif
1682- enddo
1683- enddo
1684-
1685-
1686-c s-channel min & max pt of sum of 4-momenta
1687-c
1688- do i=nincoming+1,nexternal-1
1689- do j=i+1,nexternal
1690- if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
1691- if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
1692- if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
1693- tmp=PtDot(p(0,i),p(0,j))
1694- notgood=(tmp .lt. ptll_min(j,i).or.
1695- $ ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
1696- if (notgood) then
1697- if(debug) write (*,*) i,j,' -> fails'
1698- passcuts=.false.
1699- return
1700- endif
1701- endif
1702- enddo
1703- enddo
1704-
1705-
1706-
1707-
1708-c
1709-c s-channel min & max invariant mass cuts
1710-c
1711- do i=nincoming+1,nexternal-1
1712- do j=i+1,nexternal
1713- if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
1714- if(debug) write (*,*) s_min(j,i),s_max(j,i)
1715- if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
1716- tmp=SumDot(p(0,i),p(0,j),+1d0)
1717- if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
1718- notgood=(tmp .lt. s_min(j,i).or.
1719- $ s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i))
1720- if (notgood) then
1721- if(debug) write (*,*) i,j,' -> fails'
1722- passcuts=.false.
1723- return
1724- endif
1725- else
1726- notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i))
1727- if (notgood) then
1728- if(debug) write (*,*) i,j,' -> fails'
1729- passcuts=.false.
1730- return
1731- endif
1732- endif
1733- endif
1734- enddo
1735- enddo
1736-C $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
1737-c
1738-c B.W. phase space cuts
1739-c
1740- pass_bw=cut_bw(p)
1741-c JA 4/8/11 always check pass_bw
1742- if ( pass_bw.and..not.CUTSPASSED) then
1743- passcuts=.false.
1744- if(debug) write (*,*) ' pass_bw -> fails'
1745- return
1746- endif
1747-C $E$DESACTIVATE_BW_CUT$E$ This is a Tag for MadWeight
1748- CUTSPASSED=.FALSE.
1749-C
1750-C maximal and minimal pt of the jets sorted by pt
1751-c
1752- njets=0
1753- nheavyjets=0
1754-
1755-c- fill ptjet with the pt's of the jets.
1756- do i=nincoming+1,nexternal
1757- if(is_a_j(i)) then
1758- njets=njets+1
1759- ptjet(njets)=pt(p(0,i))
1760- endif
1761- if(is_a_b(i)) then
1762- nheavyjets=nheavyjets+1
1763- ptheavyjet(nheavyjets)=pt(p(0,i))
1764- endif
1765-
1766- enddo
1767- if(debug) write (*,*) 'not yet ordered ',njets,' ',ptjet
1768-
1769-C----------------------------------------------------------------------------
1770-C DURHAM_KT CUT
1771-C----------------------------------------------------------------------------
1772-
1773- IF ( KT_DURHAM .GT. 0D0) THEN
1774-
1775-C RESET JET MOMENTA
1776- njets=0
1777- DO I=1,NEXTERNAL
1778- DO J=0,3
1779- PJET(I,J) = 0D0
1780- ENDDO
1781- ENDDO
1782-
1783- do i=nincoming+1,nexternal
1784- if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
1785- njets=njets+1
1786- DO J=0,3
1787- PJET(NJETS,J) = P(J,I)
1788- ENDDO
1789- endif
1790- enddo
1791-
1792-C COUNT NUMBER OF MASSLESS OUTGOING PARTICLES, SINCE WE DO NOT WANT
1793-C TO APPLY A CUT FOR JUST A SINGLE MASSIVE PARTICLE IN THE FINAL STATE.
1794- NMASSLESS = 0
1795- DO I=NINCOMING+1,NEXTERNAL
1796- if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and.
1797- & is_a_j(i).or.is_a_b(i)) then
1798- NMASSLESS = NMASSLESS + 1
1799- ENDIF
1800- ENDDO
1801-
1802-C DURHAM KT SEPARATION CUT
1803- IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
1804-
1805- PTMIN = EBEAM(1) + EBEAM(2)
1806-
1807- DO I=1,NJETS
1808-
1809-C PT WITH RESPECT TO Z AXIS FOR HADRONIC COLLISIONS
1810- IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
1811- PT1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2)
1812- PTMIN = MIN( PTMIN, PT1 )
1813- ENDIF
1814-
1815- DO J=I+1,NJETS
1816-C GET ANGLE BETWEEN JETS
1817- PABS1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2 + PJET(I,3)**2)
1818- PABS2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2 + PJET(J,3)**2)
1819-C CHECK IF 3-MOMENTA DO NOT VANISH
1820- IF(PABS1*PABS2 .NE. 0D0) THEN
1821- COSTH = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) + PJET(I,3)*PJET(J,3) )/(PABS1*PABS2)
1822- ELSE
1823-C IF 3-MOMENTA VANISH, MAKE JET COSTH = 1D0 SO THAT JET MEASURE VANISHES
1824- COSTH = 1D0
1825- ENDIF
1826-C GET PT AND ETA OF JETS
1827- PT2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2)
1828- ETA1 = 0.5D0*LOG( (PJET(I,0) + PJET(I,3)) / (PJET(I,0) - PJET(I,3)) )
1829- ETA2 = 0.5D0*LOG( (PJET(J,0) + PJET(J,3)) / (PJET(J,0) - PJET(J,3)) )
1830-C GET COSH OF DELTA ETA, COS OF DELTA PHI
1831- COSH_DETA = DCOSH( ETA1 - ETA2 )
1832- COS_DPHI = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) ) / (PT1*PT2)
1833- DPHI = DACOS( COS_DPHI )
1834- IF ( (LPP(1).EQ.0) .AND. (LPP(2).EQ.0)) THEN
1835-C KT FOR E+E- COLLISION
1836- PTNOW = DSQRT( 2D0*MIN(PJET(I,0)**2,PJET(J,0)**2)*( 1D0-COSTH ) )
1837- ELSE
1838-C HADRONIC KT, FASTJET DEFINITION
1839- PTNOW = DSQRT( MIN(PT1**2,PT2**2)*( (ETA1 - ETA2 )**2 + DPHI**2 )/(D_PARAMETER**2) )
1840- ENDIF
1841-
1842- PTMIN = MIN( PTMIN, PTNOW )
1843-
1844- ENDDO ! LOOP OVER NJET
1845-
1846- ENDDO ! LOOP OVER NJET
1847-
1848-C CHECK COMPATIBILITY WITH CUT
1849- IF( (PTMIN .LT. KT_DURHAM)) THEN
1850- PASSCUTS = .FALSE.
1851- RETURN
1852- ENDIF
1853-
1854- ENDIF ! IF NJETS.GT. 0
1855-
1856- ENDIF ! KT_DURHAM .GT. 0D0
1857-
1858-C----------------------------------------------------------------------------
1859-C PTLUND CUT
1860-C----------------------------------------------------------------------------
1861-
1862- IF(PT_LUND .GT. 0D0 ) THEN
1863-
1864-C Reset jet momenta
1865- NJETS=0
1866- DO I=1,NEXTERNAL
1867- JETFLAVOUR(I) = 0
1868- DO J=0,3
1869- PJET(I,J) = 0D0
1870- ENDDO
1871- ENDDO
1872-
1873-C Fill incoming particle momenta
1874- DO I=1,NINCOMING
1875- INCFLAVOUR(I) = IDUP(I,1,1)
1876- DO J=0,3
1877- PINC(I,J) = P(J,I)
1878- ENDDO
1879- ENDDO
1880-
1881-C Fill final jet momenta
1882- DO I=NINCOMING+1,NEXTERNAL
1883- if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
1884- NJETS=NJETS+1
1885- JETFLAVOUR(NJETS) = IDUP(I,1,1)
1886- DO J=0,3
1887- PJET(NJETS,J) = P(J,I)
1888- ENDDO
1889- ENDIF
1890- ENDDO
1891-
1892-C PROCESS WITH EXACTLY TWO MASSLESS OUTGOING PARTICLES IS SPECIAL
1893-C BECAUSE AN ENERGY-SHARING VARIABLE LIKE "Z" DOES NOT MAKE SENSE.
1894-C IN THIS CASE, ONLY APPLY MINIMAL pT W.R.T BEAM CUT.
1895-C THIS CUT WILL ONLY APPLY TO THE TWO-MASSLESS PARTICLE STATE.
1896- NMASSLESS = 0
1897- DO I=NINCOMING+1,NEXTERNAL
1898- if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and.
1899- & is_a_j(i).or.is_a_b(i)) THEN
1900- NMASSLESS = NMASSLESS + 1
1901- ENDIF
1902- ENDDO
1903- IF (NMASSLESS .EQ. 2 .AND. NJETS .EQ. 2 .AND.
1904- & NEXTERNAL-NINCOMING .EQ. 2) THEN
1905- PTMINSAVE = EBEAM(1) + EBEAM(2)
1906- DO I=NINCOMING+1,NEXTERNAL
1907- if( .not. from_decay(I) ) then
1908- PTMINSAVE = MIN(PTMINSAVE, PT(p(0,i)))
1909- ENDIF
1910- ENDDO
1911-C CHECK COMPATIBILITY WITH CUT
1912- IF ( ((LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND.
1913- & PTMINSAVE .LT. PT_LUND) THEN
1914- PASSCUTS = .FALSE.
1915- RETURN
1916- ENDIF
1917-C RESET NJETS TO AVOID FURTHER MERGING SCALE CUT.
1918- NJETS=0
1919- ENDIF
1920-
1921-C PYTHIA PT SEPARATION CUT
1922- IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
1923-
1924- PTMINSAVE = EBEAM(1) + EBEAM(2)
1925-
1926- DO I=1,NJETS
1927-
1928- PTMIN = EBEAM(1) + EBEAM(2)
1929- PTMINSAVE = MIN(PTMIN,PTMINSAVE)
1930-
1931-C Compute pythia ISR separation between i-jet and incoming.
1932-C Only SM-like emissions off the beam are possible.
1933- IF ( ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND.
1934- & ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
1935-C Check separation to first incoming particle
1936- DO L=0,3
1937- PRADTEMP(L) = PINC(1,L)
1938- PEMTTEMP(L) = PJET(I,L)
1939- PRECTEMP(L) = PINC(2,L)
1940- ENDDO
1941- PT1 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(1),
1942- & JETFLAVOUR(I), -1, -1)
1943- PTMIN = MIN( PTMIN, PT1 )
1944-C Check separation to second incoming particle
1945- DO L=0,3
1946- PRADTEMP(L) = PINC(2,L)
1947- PEMTTEMP(L) = PJET(I,L)
1948- PRECTEMP(L) = PINC(1,L)
1949- ENDDO
1950- PT2 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(2),
1951- & JETFLAVOUR(I), -1, -1)
1952- PTMIN = MIN( PTMIN, PT2 )
1953- ENDIF
1954-
1955-C Compute pythia FSR separation between two jets,
1956-C without any knowledge of colour connections
1957- DO J=1,NJETS
1958- DO K=1,NJETS
1959- IF ( I .NE. J .AND. I .NE. K .AND. J .NE. K ) THEN
1960-
1961-C Check separation between final partons i and j, with k as spectator
1962- DO L=0,3
1963- PRADTEMP(L) = PJET(J,L)
1964- PEMTTEMP(L) = PJET(I,L)
1965- PRECTEMP(L) = PJET(K,L)
1966- ENDDO
1967-
1968- TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
1969- & JETFLAVOUR(J), JETFLAVOUR(I), 1, 1)
1970-C Only SM-like emissions off the beam are possible, no additional
1971-C BSM particles will be produced as as shower emissions.
1972- IF ( ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
1973- PTMIN = MIN(PTMIN, TEMP);
1974- ENDIF
1975-
1976- TEMP = RHOPYTHIA( PEMTTEMP, PRADTEMP, PRECTEMP,
1977- & JETFLAVOUR(I), JETFLAVOUR(J), 1, 1)
1978-C Only SM-like emissions off the beam are possible, no additional
1979-C BSM particles will be produced as as shower emissions.
1980- IF ( ABS(JETFLAVOUR(J)) .LT. 30 ) THEN
1981- PTMIN = MIN(PTMIN, TEMP);
1982- ENDIF
1983-
1984- ENDIF
1985-
1986- ENDDO ! LOOP OVER NJET
1987- ENDDO ! LOOP OVER NJET
1988-
1989-C Compute pythia FSR separation between two jets, with initial spectator
1990- IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
1991- DO J=1,NJETS
1992-
1993-C BSM particles can only be radiators, and will not be produced
1994-C as shower emissions.
1995- IF ( ABS(JETFLAVOUR(I)) .GT. 1000000 ) THEN
1996- EXIT
1997- ENDIF
1998-
1999-C Allow both initial partons as recoiler
2000- IF ( I .NE. J ) THEN
2001-
2002-C Check with first initial as recoiler
2003- DO L=0,3
2004- PRADTEMP(L) = PJET(J,L)
2005- PEMTTEMP(L) = PJET(I,L)
2006- PRECTEMP(L) = PINC(1,L)
2007- ENDDO
2008- TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
2009- & JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
2010- IF ( LPP(1) .NE. 0 ) THEN
2011- PTMIN = MIN(PTMIN, TEMP);
2012- ENDIF
2013- DO L=0,3
2014- PRADTEMP(L) = PJET(J,L)
2015- PEMTTEMP(L) = PJET(I,L)
2016- PRECTEMP(L) = PINC(2,L)
2017- ENDDO
2018- TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
2019- & JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
2020- IF ( LPP(2) .NE. 0 ) THEN
2021- PTMIN = MIN(PTMIN, TEMP);
2022- ENDIF
2023- ENDIF
2024- ENDDO ! LOOP OVER NJET
2025- ENDIF
2026-
2027- PTMINSAVE = MIN(PTMIN,PTMINSAVE)
2028-
2029- ENDDO ! LOOP OVER NJET
2030-
2031-C CHECK COMPATIBILITY WITH CUT
2032- IF (PTMINSAVE .LT. PT_LUND) THEN
2033- PASSCUTS = .FALSE.
2034- RETURN
2035- ENDIF
2036-
2037- ENDIF ! IF NJETS.GT. 0
2038-
2039- ENDIF ! PT_LUND .GT. 0D0
2040-
2041-C----------------------------------------------------------------------------
2042-C----------------------------------------------------------------------------
2043-
2044-
2045-
2046-
2047-c- check existance of jets if jet cuts are on
2048- if(njets.lt.1.and.(htjmin.gt.0.or.ptj1min.gt.0).or.
2049- $ njets.lt.2.and.ptj2min.gt.0.or.
2050- $ njets.lt.3.and.ptj3min.gt.0.or.
2051- $ njets.lt.4.and.ptj4min.gt.0)then
2052- if(debug) write (*,*) i, ' too few jets -> fails'
2053- passcuts=.false.
2054- return
2055- endif
2056-
2057-c - sort jet pts
2058- do i=1,njets-1
2059- do j=i+1,njets
2060- if(ptjet(j).gt.ptjet(i)) then
2061- temp=ptjet(i)
2062- ptjet(i)=ptjet(j)
2063- ptjet(j)=temp
2064- endif
2065- enddo
2066- enddo
2067- if(debug) write (*,*) 'ordered ',njets,' ',ptjet
2068-c
2069-c Use "and" or "or" prescriptions
2070-c
2071- inclht=0
2072-
2073- if(njets.gt.0) then
2074-
2075- notgood=.not.jetor
2076- if(debug) write (*,*) 'jetor :',jetor
2077- if(debug) write (*,*) '0',notgood
2078-
2079- do i=1,min(njets,4)
2080- if(debug) write (*,*) i,ptjet(i), ' ',ptjmin4(i),
2081- $ ':',ptjmax4(i)
2082- if(jetor) then
2083-c--- if one of the jets does not pass, the event is rejected
2084- notgood=notgood.or.(ptjmax4(i).ge.0d0 .and.
2085- $ ptjet(i).gt.ptjmax4(i)) .or.
2086- $ (ptjet(i).lt.ptjmin4(i))
2087- if(debug) write (*,*) i,' notgood total:', notgood
2088- else
2089-c--- all cuts must fail to reject the event
2090- notgood=notgood.and.(ptjmax4(i).ge.0d0 .and.
2091- $ ptjet(i).gt.ptjmax4(i) .or.
2092- $ (ptjet(i).lt.ptjmin4(i)))
2093- if(debug) write (*,*) i,' notgood total:', notgood
2094- endif
2095- enddo
2096-
2097-
2098- if (notgood) then
2099- if(debug) write (*,*) i, ' multiple pt -> fails'
2100- passcuts=.false.
2101- return
2102- endif
2103-
2104-c---------------------------
2105-c Ht cuts
2106-C---------------------------
2107- htj=ptjet(1)
2108-
2109- do i=2,njets
2110- htj=htj+ptjet(i)
2111- if(debug) write (*,*) i, 'htj ',htj
2112- if(debug.and.i.le.4) write (*,*) 'htmin ',i,' ', htjmin4(i),':',htjmax4(i)
2113- if(i.le.4)then
2114- if(htj.lt.htjmin4(i) .or.
2115- $ htjmax4(i).ge.0d0.and.htj.gt.htjmax4(i)) then
2116- if(debug) write (*,*) i, ' ht -> fails'
2117- passcuts=.false.
2118- return
2119- endif
2120- endif
2121- enddo
2122-
2123- if(htj.lt.htjmin.or.htjmax.ge.0d0.and.htj.gt.htjmax)then
2124- if(debug) write (*,*) i, ' htj -> fails'
2125- passcuts=.false.
2126- return
2127- endif
2128-
2129- inclht=htj
2130-
2131- endif !if there are jets
2132-
2133- if(nheavyjets.gt.0) then
2134- do i=1,nheavyjets
2135- inclht=inclht+ptheavyjet(i)
2136- enddo
2137- endif !if there are heavyjets
2138-
2139- if(inclht.lt.inclHtmin.or.
2140- $ inclHtmax.ge.0d0.and.inclht.gt.inclHtmax)then
2141- if(debug) write (*,*) ' inclhtmin=',inclHtmin,' -> fails'
2142- passcuts=.false.
2143- return
2144- endif
2145-
2146-C
2147-C maximal and minimal pt of the leptons sorted by pt
2148-c
2149- nleptons=0
2150-
2151- if(ptl1min.gt.0.or.ptl2min.gt.0.or.ptl3min.gt.0.or.ptl4min.gt.0.or.
2152- $ ptl1max.ge.0d0.or.ptl2max.ge.0d0.or.
2153- $ ptl3max.ge.0d0.or.ptl4max.ge.0d0) then
2154-
2155-c - fill ptlepton with the pt's of the leptons.
2156- do i=nincoming+1,nexternal
2157- if(is_a_l(i)) then
2158- nleptons=nleptons+1
2159- ptlepton(nleptons)=pt(p(0,i))
2160- endif
2161- enddo
2162- if(debug) write (*,*) 'not yet ordered ',njets,' ',ptjet
2163-
2164-c - check existance of leptons if lepton cuts are on
2165- if(nleptons.lt.1.and.ptl1min.gt.0.or.
2166- $ nleptons.lt.2.and.ptl2min.gt.0.or.
2167- $ nleptons.lt.3.and.ptl3min.gt.0.or.
2168- $ nleptons.lt.4.and.ptl4min.gt.0)then
2169- if(debug) write (*,*) i, ' too few leptons -> fails'
2170- passcuts=.false.
2171- return
2172- endif
2173-
2174-c - sort lepton pts
2175- do i=1,nleptons-1
2176- do j=i+1,nleptons
2177- if(ptlepton(j).gt.ptlepton(i)) then
2178- temp=ptlepton(i)
2179- ptlepton(i)=ptlepton(j)
2180- ptlepton(j)=temp
2181- endif
2182- enddo
2183- enddo
2184- if(debug) write (*,*) 'ordered ',nleptons,' ',ptlepton
2185-
2186- if(nleptons.gt.0) then
2187-
2188- notgood = .false.
2189- do i=1,min(nleptons,4)
2190- if(debug) write (*,*) i,ptlepton(i), ' ',ptlmin4(i),':',ptlmax4(i)
2191-c--- if one of the leptons does not pass, the event is rejected
2192- notgood=notgood.or.
2193- $ (ptlmax4(i).ge.0d0.and.ptlepton(i).gt.ptlmax4(i)).or.
2194- $ (ptlepton(i).lt.ptlmin4(i))
2195- if(debug) write (*,*) i,' notgood total:', notgood
2196- enddo
2197-
2198-
2199- if (notgood) then
2200- if(debug) write (*,*) i, ' multiple pt -> fails'
2201- passcuts=.false.
2202- return
2203- endif
2204- endif
2205- endif
2206-C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
2207-C SPECIAL CUTS
2208-C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2209-
2210-C REQUIRE AT LEAST ONE JET WITH PT>XPTJ
2211-
2212- IF(xptj.gt.0d0) THEN
2213- xvar=0
2214- do i=nincoming+1,nexternal
2215- if(is_a_j(i)) xvar=max(xvar,pt(p(0,i)))
2216- enddo
2217- if (xvar .lt. xptj) then
2218- passcuts=.false.
2219- return
2220- endif
2221- ENDIF
2222-
2223-C REQUIRE AT LEAST ONE PHOTON WITH PT>XPTA
2224-
2225- IF(xpta.gt.0d0) THEN
2226- xvar=0
2227- do i=nincoming+1,nexternal
2228- if(is_a_a(i)) xvar=max(xvar,pt(p(0,i)))
2229- enddo
2230- if (xvar .lt. xpta) then
2231- passcuts=.false.
2232- return
2233- endif
2234- ENDIF
2235-
2236-C REQUIRE AT LEAST ONE B WITH PT>XPTB
2237-
2238- IF(xptb.gt.0d0) THEN
2239- xvar=0
2240- do i=nincoming+1,nexternal
2241- if(is_a_b(i)) xvar=max(xvar,pt(p(0,i)))
2242- enddo
2243- if (xvar .lt. xptb) then
2244- passcuts=.false.
2245- return
2246- endif
2247- ENDIF
2248-
2249-C REQUIRE AT LEAST ONE LEPTON WITH PT>XPTL
2250-
2251- IF(xptl.gt.0d0) THEN
2252- xvar=0
2253- do i=nincoming+1,nexternal
2254- if(is_a_l(i)) xvar=max(xvar,pt(p(0,i)))
2255- enddo
2256- if (xvar .lt. xptl) then
2257- passcuts=.false.
2258- if(debug) write (*,*) ' xptl -> fails'
2259- return
2260- endif
2261- ENDIF
2262-C
2263-C WBF CUTS: TWO TYPES
2264-C
2265-C FIRST TYPE: implemented by FM
2266-C
2267-C 1. FIND THE 2 HARDEST JETS
2268-C 2. REQUIRE |RAP(J)|>XETAMIN
2269-C 3. REQUIRE RAP(J1)*ETA(J2)<0
2270-C
2271-C SECOND TYPE : added by Simon de Visscher 1-08-2007
2272-C
2273-C 1. FIND THE 2 HARDEST JETS
2274-C 2. REQUIRE |RAP(J1)-RAP(J2)|>DELTAETA
2275-C 3. REQUIRE RAP(J1)*RAP(J2)<0
2276-C
2277-C
2278- hardj1=0
2279- hardj2=0
2280- ptmax1=0d0
2281- ptmax2=0d0
2282-
2283-C-- START IF AT LEAST ONE OF THE CUTS IS ACTIVATED
2284-
2285- IF(XETAMIN.GT.0D0.OR.DELTAETA.GT.0D0) THEN
2286-
2287-C-- FIND THE HARDEST JETS
2288-
2289- do i=nincoming+1,nexternal
2290- if(is_a_j(i)) then
2291-c write (*,*) i,pt(p(0,i))
2292- if(pt(p(0,i)).gt.ptmax1) then
2293- hardj2=hardj1
2294- ptmax2=ptmax1
2295- hardj1=i
2296- ptmax1=pt(p(0,i))
2297- elseif(pt(p(0,i)).gt.ptmax2) then
2298- hardj2=i
2299- ptmax2=pt(p(0,i))
2300- endif
2301-c write (*,*) hardj1,hardj2,ptmax1,ptmax2
2302- endif
2303- enddo
2304-
2305- if (hardj2.eq.0) goto 21 ! bypass vbf cut since not enough jets
2306-
2307-C-- NOW APPLY THE CUT I
2308-
2309- if (abs(rap(p(0,hardj1))) .lt. xetamin
2310- & .or.abs(rap(p(0,hardj2))) .lt. xetamin
2311- & .or.rap(p(0,hardj1))*rap(p(0,hardj2)) .gt.0d0) then
2312- passcuts=.false.
2313- return
2314- endif
2315-
2316-
2317-C-- NOW APPLY THE CUT II
2318-
2319- if (abs(rap(p(0,hardj1))-rap(p(0,hardj2))) .lt. deltaeta) then
2320- passcuts=.false.
2321- return
2322- endif
2323-
2324-c write (*,*) hardj1,hardj2,rap(p(0,hardj1)),rap(p(0,hardj2))
2325-
2326- ENDIF
2327-
2328-c Begin photon isolation
2329-c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
2330-c Use is made of parton cm frame momenta. If this must be
2331-c changed, pQCD used below must be redefined
2332-c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
2333-c If we do not require a mimimum jet energy, there's no need to apply
2334-c jet clustering and all that.
2335- if (ptgmin.ne.0d0) then
2336-
2337-c Put all (light) QCD partons in momentum array for jet clustering.
2338-c From the run_card.dat, maxjetflavor defines if b quark should be
2339-c considered here (via the logical variable 'is_a_jet'). nQCD becomes
2340-c the number of (light) QCD partons at the real-emission level (i.e. one
2341-c more than the Born).
2342-
2343- nQCD=0
2344- do j=nincoming+1,nexternal
2345- if (is_a_j(j)) then
2346- nQCD=nQCD+1
2347- do i=0,3
2348- pQCD(i,nQCD)=p(i,j) ! Use C.o.M. frame momenta
2349- enddo
2350- endif
2351- enddo
2352-
2353- nph=0
2354- do j=nincoming+1,nexternal
2355- if (is_a_a(j)) then
2356- nph=nph+1
2357- do i=0,3
2358- pgamma(i,nph)=p(i,j) ! Use C.o.M. frame momenta
2359- enddo
2360- endif
2361- enddo
2362- if(nph.eq.0) goto 444
2363-
2364- if(isoEM)then
2365- nem=nph
2366- do k=1,nem
2367- do i=0,3
2368- pem(i,k)=pgamma(i,k)
2369- enddo
2370- enddo
2371- do j=nincoming+1,nexternal
2372- if (is_a_l(j)) then
2373- nem=nem+1
2374- do i=0,3
2375- pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
2376- enddo
2377- endif
2378- enddo
2379- endif
2380-
2381- alliso=.true.
2382-
2383- j=0
2384- dowhile(j.lt.nph.and.alliso)
2385-c Loop over all photons
2386- j=j+1
2387-
2388- ptg=pt(pgamma(0,j))
2389- if(ptg.lt.ptgmin)then
2390- passcuts=.false.
2391- return
2392- endif
2393-
2394-c Isolate from hadronic energy
2395- do i=1,nQCD
2396- drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
2397- enddo
2398- call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
2399- Etsum(0)=0.d0
2400- nin=0
2401- do i=1,nQCD
2402- if(dble(drlist(isorted(i))).le.R0gamma)then
2403- nin=nin+1
2404- Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
2405- endif
2406- enddo
2407- do i=1,nin
2408- alliso=alliso .and.
2409- # Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
2410- # R0gamma,xn,epsgamma,ptg)
2411- enddo
2412-
2413-c Isolate from EM energy
2414- if(isoEM.and.nem.gt.1)then
2415- do i=1,nem
2416- drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
2417- enddo
2418- call sortzv(drlist,isorted,nem,ismode,isway,izero)
2419-c First of list must be the photon: check this, and drop it
2420- if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
2421- write(*,*)'Error #1 in photon isolation'
2422- write(*,*)j,isorted(1),drlist(isorted(1))
2423- stop
2424- endif
2425- Etsum(0)=0.d0
2426- nin=0
2427- do i=2,nem
2428- if(dble(drlist(isorted(i))).le.R0gamma)then
2429- nin=nin+1
2430- Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
2431- endif
2432- enddo
2433- do i=1,nin
2434- alliso=alliso .and.
2435- # Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
2436- # R0gamma,xn,epsgamma,ptg)
2437- enddo
2438-
2439- endif
2440-
2441-c End of loop over photons
2442- enddo
2443-
2444- if(.not.alliso)then
2445- passcuts=.false.
2446- return
2447- endif
2448- endif
2449-
2450- 444 continue
2451-c End photon isolation
2452-
2453-c
2454-c call the dummy_cuts function to check plugin/user defined cuts
2455-c
2456-
2457- if(.not.dummy_cuts(P))then
2458- passcuts=.false.
2459- return
2460- endif
2461-
2462-
2463-C...Set couplings if event passed cuts
2464-
2465- 21 if(.not.fixed_ren_scale) then
2466- call set_ren_scale(P,scale)
2467- if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
2468- endif
2469-
2470- if(.not.fixed_fac_scale) then
2471- call set_fac_scale(P,q2fact)
2472- endif
2473-
2474-c
2475-c Here we cluster event and reset factorization and renormalization
2476-c scales on an event-by-event basis, as well as check xqcut for jets
2477-c
2478-c Note the following condition is the first line of setclscales
2479-
2480- IF (FIRSTTIME2) THEN
2481- FIRSTTIME2=.FALSE.
2482- write(6,*) 'alpha_s for scale ',scale,' is ', G**2/(16d0*atan(1d0))
2483- ENDIF
2484-
2485- if(debug) write (*,*) '============================='
2486- if(debug) write (*,*) ' EVENT PASSED THE CUTS '
2487- if(debug) write (*,*) '============================='
2488-
2489- CUTSPASSED=.TRUE.
2490-
2491- RETURN
2492- END
2493-
2494-
2495-C
2496-C FUNCTION FOR ISOLATION
2497-C
2498-
2499- function iso_getdrv40(p1,p2)
2500- implicit none
2501- real*8 iso_getdrv40,p1(0:3),p2(0:3)
2502- real*8 iso_getdr
2503-c
2504- iso_getdrv40=iso_getdr(p1(0),p1(1),p1(2),p1(3),
2505- # p2(0),p2(1),p2(2),p2(3))
2506- return
2507- end
2508-
2509-
2510- function iso_getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
2511- implicit none
2512- real*8 iso_getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
2513- # iso_getpseudorap,iso_getdelphi
2514-c
2515- deta=iso_getpseudorap(en1,ptx1,pty1,pl1)-
2516- # iso_getpseudorap(en2,ptx2,pty2,pl2)
2517- dphi=iso_getdelphi(ptx1,pty1,ptx2,pty2)
2518- iso_getdr=sqrt(dphi**2+deta**2)
2519- return
2520- end
2521-
2522-
2523- function iso_getpseudorap(en,ptx,pty,pl)
2524- implicit none
2525- real*8 iso_getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
2526- parameter (tiny=1.d-5)
2527-c
2528- pt=sqrt(ptx**2+pty**2)
2529- if(pt.lt.tiny.and.abs(pl).lt.tiny)then
2530- eta=sign(1.d0,pl)*1.d8
2531- else
2532- th=atan2(pt,pl)
2533- eta=-log(tan(th/2.d0))
2534- endif
2535- iso_getpseudorap=eta
2536- return
2537- end
2538-
2539-
2540- function iso_getdelphi(ptx1,pty1,ptx2,pty2)
2541- implicit none
2542- real*8 iso_getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
2543- parameter (tiny=1.d-5)
2544-c
2545- pt1=sqrt(ptx1**2+pty1**2)
2546- pt2=sqrt(ptx2**2+pty2**2)
2547- if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
2548- tmp=ptx1*ptx2+pty1*pty2
2549- tmp=tmp/(pt1*pt2)
2550- if(abs(tmp).gt.1.d0+tiny)then
2551- write(*,*)'Cosine larger than 1'
2552- stop
2553- elseif(abs(tmp).ge.1.d0)then
2554- tmp=sign(1.d0,tmp)
2555- endif
2556- tmp=acos(tmp)
2557- else
2558- tmp=1.d8
2559- endif
2560- iso_getdelphi=tmp
2561- return
2562- end
2563-
2564- function chi_gamma_iso(dr,R0,xn,epsgamma,pTgamma)
2565-c Eq.(3.4) of Phys.Lett. B429 (1998) 369-374 [hep-ph/9801442]
2566- implicit none
2567- real*8 chi_gamma_iso,dr,R0,xn,epsgamma,pTgamma
2568- real*8 tmp,axn
2569-c
2570- axn=abs(xn)
2571- tmp=epsgamma*pTgamma
2572- if(axn.ne.0.d0)then
2573- tmp=tmp*( (1-cos(dr))/(1-cos(R0)) )**axn
2574- endif
2575- chi_gamma_iso=tmp
2576- return
2577- end
2578-
2579-
2580-*
2581-* $Id: sortzv.F,v 1.1.1.1 1996/02/15 17:49:50 mclareni Exp $
2582-*
2583-* $Log: sortzv.F,v $
2584-* Revision 1.1.1.1 1996/02/15 17:49:50 mclareni
2585-* Kernlib
2586-*
2587-*
2588-c$$$#include "kerngen/pilot.h"
2589- SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
2590-C
2591-C CERN PROGLIB# M101 SORTZV .VERSION KERNFOR 3.15 820113
2592-C ORIG. 02/10/75
2593-C
2594- DIMENSION A(N1),INDEX(N1)
2595-C
2596-C
2597- N = N1
2598- IF (N.LE.0) RETURN
2599- IF (NSORT.NE.0) GO TO 2
2600- DO 1 I=1,N
2601- 1 INDEX(I)=I
2602-C
2603- 2 IF (N.EQ.1) RETURN
2604- IF (MODE) 10,20,30
2605- 10 STOP 5 ! CALL SORTTI (A,INDEX,N)
2606- GO TO 40
2607-C
2608- 20 STOP 5 ! CALL SSORTTC(A,INDEX,N)
2609- GO TO 40
2610-C
2611- 30 CALL SORTTF (A,INDEX,N)
2612-C
2613- 40 IF (NWAY.EQ.0) GO TO 50
2614- N2 = N/2
2615- DO 41 I=1,N2
2616- ISWAP = INDEX(I)
2617- K = N+1-I
2618- INDEX(I) = INDEX(K)
2619- 41 INDEX(K) = ISWAP
2620- 50 RETURN
2621- END
2622-* ========================================
2623- SUBROUTINE SORTTF (A,INDEX,N1)
2624-C
2625- DIMENSION A(N1),INDEX(N1)
2626-C
2627- N = N1
2628- DO 3 I1=2,N
2629- I3 = I1
2630- I33 = INDEX(I3)
2631- AI = A(I33)
2632- 1 I2 = I3/2
2633- IF (I2) 3,3,2
2634- 2 I22 = INDEX(I2)
2635- IF (AI.LE.A (I22)) GO TO 3
2636- INDEX (I3) = I22
2637- I3 = I2
2638- GO TO 1
2639- 3 INDEX (I3) = I33
2640- 4 I3 = INDEX (N)
2641- INDEX (N) = INDEX (1)
2642- AI = A(I3)
2643- N = N-1
2644- IF (N-1) 12,12,5
2645- 5 I1 = 1
2646- 6 I2 = I1 + I1
2647- IF (I2.LE.N) I22= INDEX(I2)
2648- IF (I2-N) 7,9,11
2649- 7 I222 = INDEX (I2+1)
2650- IF (A(I22)-A(I222)) 8,9,9
2651- 8 I2 = I2+1
2652- I22 = I222
2653- 9 IF (AI-A(I22)) 10,11,11
2654- 10 INDEX(I1) = I22
2655- I1 = I2
2656- GO TO 6
2657- 11 INDEX (I1) = I3
2658- GO TO 4
2659- 12 INDEX (1) = I3
2660- RETURN
2661- END
2662-* ========================================
2663- SUBROUTINE SORTTI (A,INDEX,N1)
2664-C
2665- INTEGER A,AI
2666- DIMENSION A(N1),INDEX(N1)
2667-C
2668- N = N1
2669- DO 3 I1=2,N
2670- I3 = I1
2671- I33 = INDEX(I3)
2672- AI = A(I33)
2673- 1 I2 = I3/2
2674- IF (I2) 3,3,2
2675- 2 I22 = INDEX(I2)
2676- IF (AI.LE.A (I22)) GO TO 3
2677- INDEX (I3) = I22
2678- I3 = I2
2679- GO TO 1
2680- 3 INDEX (I3) = I33
2681- 4 I3 = INDEX (N)
2682- INDEX (N) = INDEX (1)
2683- AI = A(I3)
2684- N = N-1
2685- IF (N-1) 12,12,5
2686- 5 I1 = 1
2687- 6 I2 = I1 + I1
2688- IF (I2.LE.N) I22= INDEX(I2)
2689- IF (I2-N) 7,9,11
2690- 7 I222 = INDEX (I2+1)
2691- IF (A(I22)-A(I222)) 8,9,9
2692- 8 I2 = I2+1
2693- I22 = I222
2694- 9 IF (AI-A(I22)) 10,11,11
2695- 10 INDEX(I1) = I22
2696- I1 = I2
2697- GO TO 6
2698- 11 INDEX (I1) = I3
2699- GO TO 4
2700- 12 INDEX (1) = I3
2701- RETURN
2702- END
2703-* ========================================
2704- SUBROUTINE SORTTC (A,INDEX,N1)
2705-C
2706- INTEGER A,AI
2707- DIMENSION A(N1),INDEX(N1)
2708-C
2709- N = N1
2710- DO 3 I1=2,N
2711- I3 = I1
2712- I33 = INDEX(I3)
2713- AI = A(I33)
2714- 1 I2 = I3/2
2715- IF (I2) 3,3,2
2716- 2 I22 = INDEX(I2)
2717- IF(ICMPCH(AI,A(I22)))3,3,21
2718- 21 INDEX (I3) = I22
2719- I3 = I2
2720- GO TO 1
2721- 3 INDEX (I3) = I33
2722- 4 I3 = INDEX (N)
2723- INDEX (N) = INDEX (1)
2724- AI = A(I3)
2725- N = N-1
2726- IF (N-1) 12,12,5
2727- 5 I1 = 1
2728- 6 I2 = I1 + I1
2729- IF (I2.LE.N) I22= INDEX(I2)
2730- IF (I2-N) 7,9,11
2731- 7 I222 = INDEX (I2+1)
2732- IF (ICMPCH(A(I22),A(I222))) 8,9,9
2733- 8 I2 = I2+1
2734- I22 = I222
2735- 9 IF (ICMPCH(AI,A(I22))) 10,11,11
2736- 10 INDEX(I1) = I22
2737- I1 = I2
2738- GO TO 6
2739- 11 INDEX (I1) = I3
2740- GO TO 4
2741- 12 INDEX (1) = I3
2742- RETURN
2743- END
2744-* ========================================
2745- FUNCTION ICMPCH(IC1,IC2)
2746-C FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
2747-C ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
2748-C ICMPCH=0 IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
2749-C ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
2750- I1=IC1
2751- I2=IC2
2752- IF(I1.GE.0.AND.I2.GE.0)GOTO 40
2753- IF(I1.GE.0)GOTO 60
2754- IF(I2.GE.0)GOTO 80
2755- I1=-I1
2756- I2=-I2
2757- IF(I1-I2)80,70,60
2758- 40 IF(I1-I2)60,70,80
2759- 60 ICMPCH=-1
2760- RETURN
2761- 70 ICMPCH=0
2762- RETURN
2763- 80 ICMPCH=1
2764- RETURN
2765- END
2766-
2767-c************************************************************************
2768-c Returns pTLund (i.e. the Pythia8 evolution variable) between two
2769-c particles with momenta prad and pemt with momentum prec as spectator
2770-c************************************************************************
2771-
2772- DOUBLE PRECISION FUNCTION RHOPYTHIA(PRAD, PEMT, PREC, FLAVRAD,
2773- & FLAVEMT, RADTYPE, RECTYPE)
2774-
2775- IMPLICIT NONE
2776-c
2777-c Arguments
2778-c
2779- DOUBLE PRECISION PRAD(0:3),PEMT(0:3), PREC(0:3)
2780- INTEGER FLAVRAD, FLAVEMT, RADTYPE,RECTYPE
2781-c
2782-c Local
2783-c
2784- DOUBLE PRECISION Q(0:3),SUM(0:3), qBR(0:3), qAR(0:3)
2785- DOUBLE PRECISION Q2, m2Rad, m2Emt, m2Dip, qBR2, qAR2, x1, x2, z
2786- DOUBLE PRECISION m2RadAft, m2EmtAft, m2Rec, m2RadBef, m2ar, rescale
2787- DOUBLE PRECISION TEMP, lambda13, k1, k3, m2Final
2788- DOUBLE PRECISION m0u, m0d, m0c, m0s, m0t, m0b, m0w, m0z, m0x
2789- DOUBLE PRECISION PRECAFT(0:3)
2790- INTEGER emtsign
2791- INTEGER idRadBef
2792- LOGICAL allowed
2793-
2794-c-----
2795-c Begin Code
2796-c-----
2797-
2798-C Set masses. Currently no way of getting those?
2799- m0u = 0.0
2800- m0d = 0.0
2801- m0c = 1.5
2802- m0s = 0.0
2803- m0b = 4.7
2804- m0t = 172.5
2805- m0w = 80.4
2806- m0z = 91.188
2807- m0x = 400.0
2808-
2809-C Store recoiler momentum (since FI splittings require recoiler
2810-C rescaling)
2811- PRECAFT(0) = PREC(0)
2812- PRECAFT(1) = PREC(1)
2813- PRECAFT(2) = PREC(2)
2814- PRECAFT(3) = PREC(3)
2815-C Get sign of emitted momentum
2816- emtsign = 1
2817- if(radtype .eq. -1) emtsign = -1
2818-
2819-C Get virtuality
2820- Q(0) = pRad(0) + emtsign*pEmt(0)
2821- Q(1) = pRad(1) + emtsign*pEmt(1)
2822- Q(2) = pRad(2) + emtsign*pEmt(2)
2823- Q(3) = pRad(3) + emtsign*pEmt(3)
2824- Q2 = emtsign * ( Q(0)**2 - Q(1)**2 - Q(2)**2 - Q(3)**2 );
2825-
2826-C Reset allowed
2827- allowed = .true.
2828-
2829-C Splitting not possible for negative virtuality.
2830- if ( Q2 .lt. 0.0 ) allowed = .false.
2831-
2832-C Try to reconstruct flavour of radiator before emission.
2833- idRadBef = 0
2834-C gluon radiation: idBef = idAft
2835- if (abs(flavEmt) .eq. 21 .or. abs(flavEmt) .eq. 22 ) idRadBef=flavRad
2836-C final state gluon splitting: idBef = 21
2837- if (radtype .eq. 1 .and. flavEmt .eq. -flavRad) idRadBef=21
2838-C final state quark -> gluon conversion
2839- if (radtype .eq. 1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=flavEmt
2840-C initial state gluon splitting: idBef = -idEmt
2841- if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=-flavEmt
2842-C initial state gluon -> quark conversion
2843- if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. flavEmt) idRadBef=21
2844-C W-boson radiation
2845- if (flavEmt .eq. 24) idRadBef = flavRad+1
2846- if (flavEmt .eq. -24) idRadBef = flavRad-1
2847-
2848-C Set particle masses.
2849- m2RadAft = 0.0
2850- m2EmtAft = 0.0
2851- m2Rec = 0.0
2852- m2RadBef = 0.0
2853-
2854- m2RadAft = pRad(0)*pRad(0)-pRad(1)*pRad(1)-pRad(2)*pRad(2)-pRad(3)*pRad(3)
2855- m2EmtAft = pEmt(0)*pEmt(0)-pEmt(1)*pEmt(1)-pEmt(2)*pEmt(2)-pEmt(3)*pEmt(3)
2856- m2Rec = pRec(0)*pRec(0)-pRec(1)*pRec(1)-pRec(2)*pRec(2)-pRec(3)*pRec(3)
2857-
2858- if (m2RadAft .lt. 1d-4) m2RadAft = 0.0
2859- if (m2EmtAft .lt. 1d-4) m2EmtAft = 0.0
2860- if (m2Rec .lt. 1d-4) m2Rec = 0.0
2861-
2862- if (abs(flavRad) .ne. 21 .and. abs(flavRad) .ne. 22 .and.
2863- & abs(flavEmt) .ne. 24 .and.
2864- & abs(flavRad) .ne. abs(flavEmt)) then
2865- m2RadBef = m2RadAft
2866- else if (abs(flavEmt) .eq. 24) then
2867- if (idRadBef .ne. 0) then
2868- if( abs(idRadBef) .eq. 4 ) m2RadBef = m0c**2
2869- if( abs(idRadBef) .eq. 5 ) m2RadBef = m0b**2
2870- if( abs(idRadBef) .eq. 6 ) m2RadBef = m0t**2
2871- if( abs(idRadBef) .eq. 9000001 ) m2RadBef = m0x**2
2872- endif
2873- else if (radtype .eq. -1) then
2874- if (abs(flavRad) .eq. 21 .and. abs(flavEmt) .eq. 21) m2RadBef = m2EmtAft
2875- endif
2876-
2877-C Calculate dipole mass for final-state radiation.
2878- m2Final = 0.0
2879- m2Final = m2Final + (pRad(0) + pRec(0) + pEmt(0))**2
2880- m2Final = m2Final - (pRad(1) + pRec(1) + pEmt(1))**2
2881- m2Final = m2Final - (pRad(2) + pRec(2) + pEmt(2))**2
2882- m2Final = m2Final - (pRad(3) + pRec(3) + pEmt(3))**2
2883-
2884-C Final state splitting not possible for negative dipole mass.
2885- if ( radtype .eq. 1 .and. m2Final .lt. 0.0 ) allowed = .false.
2886-
2887-C Rescale recoiler for final-intial splittings.
2888- rescale = 1.0
2889- if (radtype .eq. 1 .and. rectype .eq. -1) then
2890- m2ar = m2Final - 2.0*Q2 + 2.0*m2RadBef
2891- rescale = (1.0 - (Q2 - m2RadBef) / (m2ar-m2RadBef))
2892- & /(1.0 + (Q2 - m2RadBef) / (m2ar-m2RadBef))
2893- pRecAft(0) = pRecAft(0)*rescale
2894- pRecAft(1) = pRecAft(1)*rescale
2895- pRecAft(2) = pRecAft(2)*rescale
2896- pRecAft(3) = pRecAft(3)*rescale
2897- endif
2898-
2899-C Final-initial splitting not possible for negative rescaling.
2900- if ( rescale .lt. 0.0 ) allowed = .false.
2901-
2902-C Construct dipole momentum for FSR.
2903- sum(0) = pRad(0) + pRecAft(0) + pEmt(0)
2904- sum(1) = pRad(1) + pRecAft(1) + pEmt(1)
2905- sum(2) = pRad(2) + pRecAft(2) + pEmt(2)
2906- sum(3) = pRad(3) + pRecAft(3) + pEmt(3)
2907- m2Dip = sum(0)**2 - sum(1)**2 - sum(2)**2 - sum(3)**2
2908-
2909-C Construct 2->3 variables for FSR
2910- x1 = 2. * ( sum(0)*pRad(0) - sum(1)*pRad(1)
2911- & - sum(2)*pRad(2) - sum(3)*pRad(3) ) / m2Dip
2912- x2 = 2. * ( sum(0)*pRecAft(0) - sum(1)*pRecAft(1)
2913- & - sum(2)*pRecAft(2) - sum(3)*pRecAft(3) ) / m2Dip
2914-
2915-C Final state splitting not possible for ill-defined
2916-C 3-body-variables.
2917- if ( radtype .eq. 1 .and. x1 .lt. 0.0 ) allowed = .false.
2918- if ( radtype .eq. 1 .and. x1 .gt. 1.0 ) allowed = .false.
2919- if ( radtype .eq. 1 .and. x2 .lt. 0.0 ) allowed = .false.
2920- if ( radtype .eq. 1 .and. x2 .gt. 1.0 ) allowed = .false.
2921-
2922-C Auxiliary variables for massive FSR
2923- lambda13 = DSQRT( (Q2 - m2RadAft - m2EmtAft )**2 - 4.0 * m2RadAft*m2EmtAft)
2924- k1 = ( Q2 - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
2925- k3 = ( Q2 - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
2926-
2927-C Construct momenta of dipole before/after splitting for ISR
2928- qBR(0) = pRad(0) + pRec(0) - pEmt(0)
2929- qBR(1) = pRad(1) + pRec(1) - pEmt(1)
2930- qBR(2) = pRad(2) + pRec(2) - pEmt(2)
2931- qBR(3) = pRad(3) + pRec(3) - pEmt(3)
2932- qBR2 = qBR(0)**2 - qBR(1)**2 - qBR(2)**2 - qBR(3)**2
2933-
2934- qAR(0) = pRad(0) + pRec(0)
2935- qAR(1) = pRad(1) + pRec(1)
2936- qAR(2) = pRad(2) + pRec(2)
2937- qAR(3) = pRad(3) + pRec(3)
2938- qAR2 = qAR(0)**2 - qAR(1)**2 - qAR(2)**2 - qAR(3)**2
2939-
2940-C Calculate z of splitting, different for FSR and ISR
2941- z = 1.0 / (1.0 - k1 -k3) * ( x1 / (2.0-x2) - k3)
2942- if(radtype .eq. -1 ) z = qBR2 / qAR2;
2943-
2944-C Splitting not possible for ill-defined energy sharing.
2945- if ( z .lt. 0.0 .or. z .gt. 1.0 ) allowed = .false.
2946-
2947-C pT^2 = separation * virtuality (corrected with mass for FSR)
2948- if (radtype .eq. 1) temp = z*(1-z)*(Q2 - m2RadBef)
2949- if (radtype .eq. -1) temp = (1-z)*Q2
2950-
2951-C Check threshold in ISR
2952- if (radtype .ne. 1) then
2953- if ((abs(flavRad) .eq. 4 .or. abs(flavEmt) .eq. 4)
2954- & .and. dsqrt(temp) .le. 2.0*m0c**2 ) temp = (1.-z)*(Q2+m0c**2)
2955- if ((abs(flavRad) .eq. 5 .or. abs(flavEmt) .eq. 5)
2956- & .and. dsqrt(temp) .le. 2.0*m0b**2 ) temp = (1.-z)*(Q2+m0b**2)
2957- endif
2958-
2959-C Kinematically impossible splittings should not be included in the
2960-C pT definition!
2961- if( .not. allowed) temp = 1d15
2962-
2963- if(temp .lt. 0.0) temp = 0.0
2964-
2965-C Return pT
2966- rhoPythia = dsqrt(temp);
2967-
2968- RETURN
2969- END
2970+ logical function pass_point(p)
2971+c************************************************************************
2972+c This function is called from sample to see if it needs to
2973+c bother calculating the weight from all the different conficurations
2974+c You can either just return true, or have it call passcuts
2975+c************************************************************************
2976+ implicit none
2977+c
2978+c Arguments
2979+c
2980+ double precision p
2981+c
2982+c External
2983+c
2984+ logical passcuts
2985+ external passcuts
2986+c-----
2987+c Begin Code
2988+c-----
2989+ pass_point = .true.
2990+c pass_point = passcuts(p)
2991+ end
2992+C
2993+ LOGICAL FUNCTION PASSCUTS(P)
2994+C**************************************************************************
2995+C INPUT:
2996+C P(0:3,1) MOMENTUM OF INCOMING PARTON
2997+C P(0:3,2) MOMENTUM OF INCOMING PARTON
2998+C P(0:3,3) MOMENTUM OF ...
2999+C ALL MOMENTA ARE IN THE REST FRAME!!
3000+C COMMON/JETCUTS/ CUTS ON JETS
3001+C OUTPUT:
3002+C TRUE IF EVENTS PASSES ALL CUTS LISTED
3003+C**************************************************************************
3004+ IMPLICIT NONE
3005+c
3006+c Constants
3007+c
3008+ include 'maxparticles.inc'
3009+ include 'nexternal.inc'
3010+C
3011+C ARGUMENTS
3012+C
3013+ REAL*8 P(0:3,nexternal)
3014+
3015+C
3016+C LOCAL
3017+C
3018+ LOGICAL FIRSTTIME,FIRSTTIME2,pass_bw,notgood,good,foundheavy
3019+ LOGICAL DEBUG
3020+ integer i,j,njets,nheavyjets,nleptons,hardj1,hardj2
3021+ REAL*8 XVAR,ptmax1,ptmax2,htj,tmp,inclht
3022+ real*8 ptemp(0:3), ptemp2(0:3)
3023+ character*20 formstr
3024+C
3025+C PARAMETERS
3026+C
3027+ real*8 PI
3028+ parameter( PI = 3.14159265358979323846d0 )
3029+C
3030+C EXTERNAL
3031+C
3032+ REAL*8 R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,PtDot
3033+ logical cut_bw,setclscales,dummy_cuts
3034+ external R2,DOT,ET,RAP,DJ,SumDot,pt,ALPHAS,cut_bw,setclscales,PtDot
3035+ external dummy_cuts
3036+C
3037+C GLOBAL
3038+C
3039+ include 'run.inc'
3040+ include 'cuts.inc'
3041+
3042+ double precision ptjet(nexternal)
3043+ double precision ptheavyjet(nexternal)
3044+ double precision ptlepton(nexternal)
3045+ double precision temp
3046+
3047+C VARIABLES TO SPECIFY JETS
3048+ DOUBLE PRECISION PJET(NEXTERNAL,0:3)
3049+ DOUBLE PRECISION PTMIN
3050+ DOUBLE PRECISION PT1,PT2
3051+
3052+C INTEGERS FOR COUNTING.
3053+ INTEGER K,L,J1,J2
3054+
3055+C VARIABLES FOR KT CUT
3056+ DOUBLE PRECISION PTNOW,COSTH,PABS1,PABS2
3057+ DOUBLE PRECISION ETA1,ETA2,COSH_DETA,COS_DPHI,KT1SQ,KT2SQ, DPHI
3058+
3059+ double precision etmin(nincoming+1:nexternal),etamax(nincoming+1:nexternal)
3060+ double precision emin(nincoming+1:nexternal)
3061+ double precision r2min(nincoming+1:nexternal,nincoming+1:nexternal)
3062+ double precision s_min(nexternal,nexternal)
3063+ double precision etmax(nincoming+1:nexternal),etamin(nincoming+1:nexternal)
3064+ double precision emax(nincoming+1:nexternal)
3065+ double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
3066+ double precision s_max(nexternal,nexternal)
3067+ double precision ptll_min(nexternal,nexternal),ptll_max(nexternal,nexternal)
3068+ double precision inclHtmin,inclHtmax
3069+ common/to_cuts/ etmin, emin, etamax, r2min, s_min,
3070+ $ etmax, emax, etamin, r2max, s_max, ptll_min, ptll_max, inclHtmin,inclHtmax
3071+
3072+ double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
3073+ logical jetor
3074+ common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
3075+
3076+ double precision ptlmin4(4),ptlmax4(4)
3077+ common/to_lepton_cuts/ ptlmin4,ptlmax4
3078+
3079+c
3080+c Special cuts
3081+c
3082+
3083+ integer lbw(0:nexternal) !Use of B.W.
3084+ common /to_BW/ lbw
3085+C
3086+C SPECIAL CUTS
3087+C
3088+ LOGICAL IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
3089+ LOGICAL IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
3090+ LOGICAL IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
3091+ logical do_cuts(nexternal)
3092+ COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,
3093+ . IS_A_ONIUM, do_cuts
3094+
3095+C
3096+C MERGING SCALE CUT
3097+C
3098+C Retrieve which external particles undergo the ktdurham and ptlund cuts.
3099+ LOGICAL is_pdg_for_merging_cut(NEXTERNAL)
3100+ logical from_decay(-(nexternal+3):nexternal)
3101+ COMMON /TO_MERGE_CUTS/is_pdg_for_merging_cut, from_decay
3102+
3103+C
3104+C ADDITIONAL VARIABLES FOR PTLUND CUT
3105+C
3106+ INTEGER NMASSLESS
3107+ DOUBLE PRECISION PINC(NINCOMING,0:3)
3108+ DOUBLE PRECISION PRADTEMP(0:3), PRECTEMP(0:3), PEMTTEMP(0:3)
3109+ DOUBLE PRECISION PTMINSAVE, RHOPYTHIA
3110+ EXTERNAL RHOPYTHIA
3111+C
3112+C FLAVOUR INFORMATION NECESSARY TO RECONSTRUCT PTLUND
3113+C
3114+ INTEGER JETFLAVOUR(NEXTERNAL), INCFLAVOUR(NINCOMING)
3115+ include 'maxamps.inc'
3116+ integer idup(nexternal,maxproc,maxsproc)
3117+ integer mothup(2,nexternal)
3118+ integer icolup(2,nexternal,maxflow,maxsproc)
3119+ include 'leshouche.inc'
3120+
3121+C
3122+C Keep track of whether cuts already calculated for this event
3123+C
3124+ LOGICAL CUTSDONE,CUTSPASSED
3125+ COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
3126+ DATA CUTSDONE,CUTSPASSED/.FALSE.,.FALSE./
3127+
3128+C $B$ MW_NEW_DEF $E$ !this is a tag for MadWeight
3129+
3130+ double precision xqcutij(nexternal,nexternal),xqcuti(nexternal)
3131+ common/to_xqcuts/xqcutij,xqcuti
3132+
3133+c jet cluster algorithm
3134+ integer nQCD !,NJET,JET(nexternal)
3135+c double precision plab(0:3, nexternal)
3136+ double precision pQCD(0:3,nexternal)!,PJET(0:3,nexternal)
3137+c double precision rfj,sycut,palg,fastjetdmerge
3138+c integer njet_eta
3139+c Photon isolation
3140+ integer nph,nem,nin
3141+ double precision ptg,chi_gamma_iso,iso_getdrv40
3142+ double precision Etsum(0:nexternal)
3143+ real drlist(nexternal)
3144+ double precision pgamma(0:3,nexternal),pem(0:3,nexternal)
3145+ logical alliso
3146+C Sort array of results: ismode>0 for real, isway=0 for ascending order
3147+ integer ismode,isway,izero,isorted(nexternal)
3148+ parameter (ismode=1)
3149+ parameter (isway=0)
3150+ parameter (izero=0)
3151+
3152+ include 'coupl.inc'
3153+
3154+C
3155+C
3156+c
3157+ DATA FIRSTTIME,FIRSTTIME2/.TRUE.,.TRUE./
3158+
3159+c put momenta in common block for couplings.f
3160+ double precision pp(0:3,max_particles)
3161+ common /momenta_pp/pp
3162+
3163+ DATA DEBUG/.FALSE./
3164+
3165+C-----
3166+C BEGIN CODE
3167+C-----
3168+
3169+
3170+
3171+ PASSCUTS=.TRUE. !EVENT IS OK UNLESS OTHERWISE CHANGED
3172+ IF (FIRSTTIME) THEN
3173+ FIRSTTIME=.FALSE.
3174+c Preparation for reweighting by setting up clustering by diagrams
3175+ ! Remove for MW!call initcluster()
3176+c
3177+c
3178+ write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'i8)'
3179+ write(*,formstr) 'Particle',(i,i=nincoming+1,nexternal)
3180+ write(formstr,'(a,i2.2,a)')'(a10,',nexternal,'f8.1)'
3181+ write(*,formstr) 'Et >',(etmin(i),i=nincoming+1,nexternal)
3182+ write(*,formstr) 'E >',(emin(i),i=nincoming+1,nexternal)
3183+ write(*,formstr) 'Eta <',(etamax(i),i=nincoming+1,nexternal)
3184+ write(*,formstr) 'xqcut: ',(xqcuti(i),i=nincoming+1,nexternal)
3185+ write(formstr,'(a,i2.2,a)')'(a,i2,a,',nexternal,'f8.1)'
3186+ do j=nincoming+1,nexternal-1
3187+ write(*,formstr) 'd R #',j,' >',(-0.0,i=nincoming+1,j),
3188+ & (r2min(i,j),i=j+1,nexternal)
3189+ do i=j+1,nexternal
3190+ r2min(i,j)=r2min(i,j)*dabs(r2min(i,j)) !Since r2 returns distance squared
3191+ r2max(i,j)=r2max(i,j)*dabs(r2max(i,j))
3192+ enddo
3193+ enddo
3194+ do j=nincoming+1,nexternal-1
3195+ write(*,formstr) 's min #',j,'>',
3196+ & (s_min(i,j),i=nincoming+1,nexternal)
3197+ enddo
3198+ do j=nincoming+1,nexternal-1
3199+ write(*,formstr) 'xqcutij #',j,'>',
3200+ & (xqcutij(i,j),i=nincoming+1,nexternal)
3201+ enddo
3202+
3203+cc
3204+cc Set the strong coupling
3205+cc
3206+c call set_ren_scale(P,scale)
3207+c
3208+cc Check that the user funtions for setting the scales
3209+cc have been edited if the choice of an event-by-event
3210+cc scale choice has been made
3211+c
3212+c if(.not.fixed_ren_scale) then
3213+c if(scale.eq.0d0) then
3214+c write(6,*)
3215+c write(6,*) '* >>>>>>>>>ERROR<<<<<<<<<<<<<<<<<<<<<<<*'
3216+c write(6,*) ' Dynamical renormalization scale choice '
3217+c write(6,*) ' selected but user subroutine'
3218+c write(6,*) ' set_ren_scale not edited in file:setpara.f'
3219+c write(6,*) ' Switching to a fixed_ren_scale choice'
3220+c write(6,*) ' with scale=zmass'
3221+c scale=91.2d0
3222+c write(6,*) 'scale=',scale
3223+c fixed_ren_scale=.true.
3224+c call set_ren_scale(P,scale)
3225+c endif
3226+c endif
3227+
3228+
3229+ if(fixed_ren_scale) then
3230+ G = SQRT(4d0*PI*ALPHAS(scale))
3231+ call update_as_param()
3232+ endif
3233+
3234+c Put momenta in the common block to zero to start
3235+ do i=0,3
3236+ do j=1,max_particles
3237+ pp(i,j) = 0d0
3238+ enddo
3239+ enddo
3240+
3241+ ENDIF ! IF FIRSTTIME
3242+
3243+ IF (CUTSDONE) THEN
3244+ PASSCUTS=CUTSPASSED
3245+ RETURN
3246+ ENDIF
3247+ CUTSDONE=.TRUE.
3248+c CUTSPASSED=.FALSE.
3249+
3250+c
3251+c Make sure have reasonable 4-momenta
3252+c
3253+ if (p(0,1) .le. 0d0) then
3254+ passcuts=.false.
3255+ return
3256+ endif
3257+
3258+c Also make sure there's no INF or NAN
3259+ do i=1,nexternal
3260+ do j=0,3
3261+ if(p(j,i).gt.1d32.or.p(j,i).ne.p(j,i))then
3262+ passcuts=.false.
3263+ return
3264+ endif
3265+ enddo
3266+ enddo
3267+
3268+c
3269+c Limit S_hat
3270+c
3271+c if (x1*x2*stot .gt. 500**2) then
3272+c passcuts=.false.
3273+c return
3274+c endif
3275+
3276+C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
3277+
3278+ if(debug) write (*,*) '============================='
3279+ if(debug) write (*,*) ' EVENT STARTS TO BE CHECKED '
3280+ if(debug) write (*,*) '============================='
3281+c
3282+c p_t min & max cuts
3283+c
3284+ do i=nincoming+1,nexternal
3285+ if(debug) write (*,*) 'pt(',i,')=',pt(p(0,i)),' ',etmin(i),
3286+ $ ':',etmax(i)
3287+ notgood=(pt(p(0,i)) .lt. etmin(i)).or.
3288+ & (etmax(i).ge.0d0.and.pt(p(0,i)) .gt. etmax(i))
3289+ if (notgood) then
3290+ if(debug) write (*,*) i,' -> fails'
3291+ passcuts=.false.
3292+ return
3293+ endif
3294+ enddo
3295+c
3296+c missing ET min & max cut + Invariant mass of leptons and neutrino
3297+c nb: missing Et defined as the vector sum over the neutrino's pt
3298+c
3299+c-- reset ptemp(0:3)
3300+ do j=0,3
3301+ ptemp(j)=0 ! for the neutrino
3302+ ptemp2(j)=0 ! for the leptons
3303+ enddo
3304+c- sum over the momenta
3305+ do i=nincoming+1,nexternal
3306+ if(is_a_nu(i)) then
3307+ if(debug) write (*,*) i,' -> neutrino '
3308+ do j=0,3
3309+ ptemp(j)=ptemp(j)+p(j,i)
3310+ enddo
3311+ elseif(is_a_l(i)) then
3312+ if(debug) write (*,*) i,' -> lepton '
3313+ do j=0,3
3314+ ptemp2(j)=ptemp2(j)+p(j,i)
3315+ enddo
3316+ endif
3317+
3318+ enddo
3319+c- check the et
3320+ if(debug.and.ptemp(0).eq.0d0) write (*,*) 'No et miss in event'
3321+ if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Et miss =',pt(ptemp(0)),' ',misset,':',missetmax
3322+ if(debug.and.ptemp2(0).eq.0d0) write (*,*) 'No leptons in event'
3323+ if(debug.and.ptemp(0).gt.0d0) write (*,*) 'Energy of leptons =',pt(ptemp2(0))
3324+ if(ptemp(0).gt.0d0) then
3325+ notgood=(pt(ptemp(0)) .lt. misset).or.
3326+ & (missetmax.ge.0d0.and.pt(ptemp(0)) .gt. missetmax)
3327+ if (notgood) then
3328+ if(debug) write (*,*) ' missing et cut -> fails'
3329+ passcuts=.false.
3330+ return
3331+ endif
3332+ endif
3333+ if (mmnl.gt.0d0.or.mmnlmax.ge.0d0)then
3334+ if(dsqrt(SumDot(ptemp,ptemp2,1d0)).lt.mmnl.or.mmnlmax.ge.0d0.and.dsqrt(SumDot(ptemp, ptemp2,1d0)).gt.mmnlmax) then
3335+ if(debug) write (*,*) 'lepton invariant mass -> fails'
3336+ passcuts=.false.
3337+ return
3338+ endif
3339+ endif
3340+c
3341+c pt cut on heavy particles
3342+c gives min(pt) for (at least) one heavy particle
3343+c
3344+ if(ptheavy.gt.0d0)then
3345+ passcuts=.false.
3346+ foundheavy=.false.
3347+ do i=nincoming+1,nexternal
3348+ if(is_heavy(i)) then
3349+ if(debug) write (*,*) i,' -> heavy '
3350+ foundheavy=.true.
3351+ if(pt(p(0,i)).gt.ptheavy) passcuts=.true.
3352+ endif
3353+ enddo
3354+
3355+ if(.not.passcuts.and.foundheavy)then
3356+ if(debug) write (*,*) ' heavy particle cut -> fails'
3357+ return
3358+ else
3359+ passcuts=.true.
3360+ endif
3361+ endif
3362+c
3363+c E min & max cuts
3364+c
3365+ do i=nincoming+1,nexternal
3366+ if(debug) write (*,*) 'p(0,',i,')=',p(0,i),' ',emin(i),':',emax(i)
3367+ notgood=(p(0,i) .le. emin(i)).or.
3368+ & (emax(i).ge.0d0 .and. p(0,i) .gt. emax(i))
3369+ if (notgood) then
3370+ if(debug) write (*,*) i,' -> fails'
3371+ passcuts=.false.
3372+ return
3373+ endif
3374+ enddo
3375+c
3376+c Rapidity min & max cuts
3377+c
3378+ do i=nincoming+1,nexternal
3379+ if(debug) write (*,*) 'abs(rap(',i,'))=',abs(rap(p(0,i))),' ',etamin(i),':',etamax(i)
3380+ notgood=(etamax(i).ge.0.and.abs(rap(p(0,i))) .gt. etamax(i)).or.
3381+ & (abs(rap(p(0,i))) .lt. etamin(i))
3382+ if (notgood) then
3383+ if(debug) write (*,*) i,' -> fails'
3384+ passcuts=.false.
3385+ return
3386+ endif
3387+ enddo
3388+c
3389+c DeltaR min & max cuts
3390+c
3391+ do i=nincoming+1,nexternal-1
3392+ do j=i+1,nexternal
3393+ if(debug) write (*,*) 'r2(',i, ',' ,j,')=',dsqrt(r2(p(0,i),p(0,j)))
3394+ if(debug) write (*,*) dsqrt(r2min(j,i)),dsqrt(r2max(j,i))
3395+ if(r2min(j,i).gt.0.or.r2max(j,i).ge.0d0) then
3396+ tmp=r2(p(0,i),p(0,j))
3397+ notgood=(tmp .lt. r2min(j,i)).or.
3398+ $ (r2max(j,i).ge.0d0 .and. tmp .gt. r2max(j,i))
3399+ if (notgood) then
3400+ if(debug) write (*,*) i,j,' -> fails'
3401+ passcuts=.false.
3402+ return
3403+ endif
3404+ endif
3405+ enddo
3406+ enddo
3407+
3408+
3409+c s-channel min & max pt of sum of 4-momenta
3410+c
3411+ do i=nincoming+1,nexternal-1
3412+ do j=i+1,nexternal
3413+ if(debug)write (*,*) 'ptll(',i,',',j,')=',PtDot(p(0,i),p(0,j))
3414+ if(debug)write (*,*) ptll_min(j,i),ptll_max(j,i)
3415+ if(ptll_min(j,i).gt.0.or.ptll_max(j,i).ge.0d0) then
3416+ tmp=PtDot(p(0,i),p(0,j))
3417+ notgood=(tmp .lt. ptll_min(j,i).or.
3418+ $ ptll_max(j,i).ge.0d0 .and. tmp.gt.ptll_max(j,i))
3419+ if (notgood) then
3420+ if(debug) write (*,*) i,j,' -> fails'
3421+ passcuts=.false.
3422+ return
3423+ endif
3424+ endif
3425+ enddo
3426+ enddo
3427+
3428+
3429+
3430+
3431+c
3432+c s-channel min & max invariant mass cuts
3433+c
3434+ do i=nincoming+1,nexternal-1
3435+ do j=i+1,nexternal
3436+ if(debug) write (*,*) 's(',i,',',j,')=',Sumdot(p(0,i),p(0,j),+1d0)
3437+ if(debug) write (*,*) s_min(j,i),s_max(j,i)
3438+ if(s_min(j,i).gt.0.or.s_max(j,i).ge.0d0) then
3439+ tmp=SumDot(p(0,i),p(0,j),+1d0)
3440+ if(s_min(j,i).le.s_max(j,i) .or. s_max(j,i).lt.0d0)then
3441+ notgood=(tmp .lt. s_min(j,i).or.
3442+ $ s_max(j,i).ge.0d0 .and. tmp .gt. s_max(j,i))
3443+ if (notgood) then
3444+ if(debug) write (*,*) i,j,' -> fails'
3445+ passcuts=.false.
3446+ return
3447+ endif
3448+ else
3449+ notgood=(tmp .lt. s_min(j,i).and.tmp .gt. s_max(j,i))
3450+ if (notgood) then
3451+ if(debug) write (*,*) i,j,' -> fails'
3452+ passcuts=.false.
3453+ return
3454+ endif
3455+ endif
3456+ endif
3457+ enddo
3458+ enddo
3459+C $B$DESACTIVATE_BW_CUT$B$ This is a Tag for MadWeight
3460+c
3461+c B.W. phase space cuts
3462+c
3463+ pass_bw=cut_bw(p)
3464+c JA 4/8/11 always check pass_bw
3465+ if ( pass_bw.and..not.CUTSPASSED) then
3466+ passcuts=.false.
3467+ if(debug) write (*,*) ' pass_bw -> fails'
3468+ return
3469+ endif
3470+C $E$DESACTIVATE_BW_CUT$E$ This is a Tag for MadWeight
3471+ CUTSPASSED=.FALSE.
3472+C
3473+C maximal and minimal pt of the jets sorted by pt
3474+c
3475+ njets=0
3476+ nheavyjets=0
3477+
3478+c- fill ptjet with the pt's of the jets.
3479+ do i=nincoming+1,nexternal
3480+ if(is_a_j(i)) then
3481+ njets=njets+1
3482+ ptjet(njets)=pt(p(0,i))
3483+ endif
3484+ if(is_a_b(i)) then
3485+ nheavyjets=nheavyjets+1
3486+ ptheavyjet(nheavyjets)=pt(p(0,i))
3487+ endif
3488+
3489+ enddo
3490+ if(debug) write (*,*) 'not yet ordered ',njets,' ',ptjet
3491+
3492+C----------------------------------------------------------------------------
3493+C DURHAM_KT CUT
3494+C----------------------------------------------------------------------------
3495+
3496+ IF ( KT_DURHAM .GT. 0D0) THEN
3497+
3498+C RESET JET MOMENTA
3499+ njets=0
3500+ DO I=1,NEXTERNAL
3501+ DO J=0,3
3502+ PJET(I,J) = 0D0
3503+ ENDDO
3504+ ENDDO
3505+
3506+ do i=nincoming+1,nexternal
3507+ if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
3508+ njets=njets+1
3509+ DO J=0,3
3510+ PJET(NJETS,J) = P(J,I)
3511+ ENDDO
3512+ endif
3513+ enddo
3514+
3515+C COUNT NUMBER OF MASSLESS OUTGOING PARTICLES, SINCE WE DO NOT WANT
3516+C TO APPLY A CUT FOR JUST A SINGLE MASSIVE PARTICLE IN THE FINAL STATE.
3517+ NMASSLESS = 0
3518+ DO I=NINCOMING+1,NEXTERNAL
3519+ if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and.
3520+ & is_a_j(i).or.is_a_b(i)) then
3521+ NMASSLESS = NMASSLESS + 1
3522+ ENDIF
3523+ ENDDO
3524+
3525+C DURHAM KT SEPARATION CUT
3526+ IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
3527+
3528+ PTMIN = EBEAM(1) + EBEAM(2)
3529+
3530+ DO I=1,NJETS
3531+
3532+C PT WITH RESPECT TO Z AXIS FOR HADRONIC COLLISIONS
3533+ IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
3534+ PT1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2)
3535+ PTMIN = MIN( PTMIN, PT1 )
3536+ ENDIF
3537+
3538+ DO J=I+1,NJETS
3539+C GET ANGLE BETWEEN JETS
3540+ PABS1 = DSQRT(PJET(I,1)**2 + PJET(I,2)**2 + PJET(I,3)**2)
3541+ PABS2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2 + PJET(J,3)**2)
3542+C CHECK IF 3-MOMENTA DO NOT VANISH
3543+ IF(PABS1*PABS2 .NE. 0D0) THEN
3544+ COSTH = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) + PJET(I,3)*PJET(J,3) )/(PABS1*PABS2)
3545+ ELSE
3546+C IF 3-MOMENTA VANISH, MAKE JET COSTH = 1D0 SO THAT JET MEASURE VANISHES
3547+ COSTH = 1D0
3548+ ENDIF
3549+C GET PT AND ETA OF JETS
3550+ PT2 = DSQRT(PJET(J,1)**2 + PJET(J,2)**2)
3551+ ETA1 = 0.5D0*LOG( (PJET(I,0) + PJET(I,3)) / (PJET(I,0) - PJET(I,3)) )
3552+ ETA2 = 0.5D0*LOG( (PJET(J,0) + PJET(J,3)) / (PJET(J,0) - PJET(J,3)) )
3553+C GET COSH OF DELTA ETA, COS OF DELTA PHI
3554+ COSH_DETA = DCOSH( ETA1 - ETA2 )
3555+ COS_DPHI = ( PJET(I,1)*PJET(J,1) + PJET(I,2)*PJET(J,2) ) / (PT1*PT2)
3556+ DPHI = DACOS( COS_DPHI )
3557+ IF ( (LPP(1).EQ.0) .AND. (LPP(2).EQ.0)) THEN
3558+C KT FOR E+E- COLLISION
3559+ PTNOW = DSQRT( 2D0*MIN(PJET(I,0)**2,PJET(J,0)**2)*( 1D0-COSTH ) )
3560+ ELSE
3561+C HADRONIC KT, FASTJET DEFINITION
3562+ PTNOW = DSQRT( MIN(PT1**2,PT2**2)*( (ETA1 - ETA2 )**2 + DPHI**2 )/(D_PARAMETER**2) )
3563+ ENDIF
3564+
3565+ PTMIN = MIN( PTMIN, PTNOW )
3566+
3567+ ENDDO ! LOOP OVER NJET
3568+
3569+ ENDDO ! LOOP OVER NJET
3570+
3571+C CHECK COMPATIBILITY WITH CUT
3572+ IF( (PTMIN .LT. KT_DURHAM)) THEN
3573+ PASSCUTS = .FALSE.
3574+ RETURN
3575+ ENDIF
3576+
3577+ ENDIF ! IF NJETS.GT. 0
3578+
3579+ ENDIF ! KT_DURHAM .GT. 0D0
3580+
3581+C----------------------------------------------------------------------------
3582+C PTLUND CUT
3583+C----------------------------------------------------------------------------
3584+
3585+ IF(PT_LUND .GT. 0D0 ) THEN
3586+
3587+C Reset jet momenta
3588+ NJETS=0
3589+ DO I=1,NEXTERNAL
3590+ JETFLAVOUR(I) = 0
3591+ DO J=0,3
3592+ PJET(I,J) = 0D0
3593+ ENDDO
3594+ ENDDO
3595+
3596+C Fill incoming particle momenta
3597+ DO I=1,NINCOMING
3598+ INCFLAVOUR(I) = IDUP(I,1,1)
3599+ DO J=0,3
3600+ PINC(I,J) = P(J,I)
3601+ ENDDO
3602+ ENDDO
3603+
3604+C Fill final jet momenta
3605+ DO I=NINCOMING+1,NEXTERNAL
3606+ if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) ) then
3607+ NJETS=NJETS+1
3608+ JETFLAVOUR(NJETS) = IDUP(I,1,1)
3609+ DO J=0,3
3610+ PJET(NJETS,J) = P(J,I)
3611+ ENDDO
3612+ ENDIF
3613+ ENDDO
3614+
3615+C PROCESS WITH EXACTLY TWO MASSLESS OUTGOING PARTICLES IS SPECIAL
3616+C BECAUSE AN ENERGY-SHARING VARIABLE LIKE "Z" DOES NOT MAKE SENSE.
3617+C IN THIS CASE, ONLY APPLY MINIMAL pT W.R.T BEAM CUT.
3618+C THIS CUT WILL ONLY APPLY TO THE TWO-MASSLESS PARTICLE STATE.
3619+ NMASSLESS = 0
3620+ DO I=NINCOMING+1,NEXTERNAL
3621+ if(is_pdg_for_merging_cut(i) .and. .not. from_decay(I) .and.
3622+ & is_a_j(i).or.is_a_b(i)) THEN
3623+ NMASSLESS = NMASSLESS + 1
3624+ ENDIF
3625+ ENDDO
3626+ IF (NMASSLESS .EQ. 2 .AND. NJETS .EQ. 2 .AND.
3627+ & NEXTERNAL-NINCOMING .EQ. 2) THEN
3628+ PTMINSAVE = EBEAM(1) + EBEAM(2)
3629+ DO I=NINCOMING+1,NEXTERNAL
3630+ if( .not. from_decay(I) ) then
3631+ PTMINSAVE = MIN(PTMINSAVE, PT(p(0,i)))
3632+ ENDIF
3633+ ENDDO
3634+C CHECK COMPATIBILITY WITH CUT
3635+ IF ( ((LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND.
3636+ & PTMINSAVE .LT. PT_LUND) THEN
3637+ PASSCUTS = .FALSE.
3638+ RETURN
3639+ ENDIF
3640+C RESET NJETS TO AVOID FURTHER MERGING SCALE CUT.
3641+ NJETS=0
3642+ ENDIF
3643+
3644+C PYTHIA PT SEPARATION CUT
3645+ IF(NJETS.GT.0 .AND. NMASSLESS .GT. 0) THEN
3646+
3647+ PTMINSAVE = EBEAM(1) + EBEAM(2)
3648+
3649+ DO I=1,NJETS
3650+
3651+ PTMIN = EBEAM(1) + EBEAM(2)
3652+ PTMINSAVE = MIN(PTMIN,PTMINSAVE)
3653+
3654+C Compute pythia ISR separation between i-jet and incoming.
3655+C Only SM-like emissions off the beam are possible.
3656+ IF ( ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) .AND.
3657+ & ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
3658+C Check separation to first incoming particle
3659+ DO L=0,3
3660+ PRADTEMP(L) = PINC(1,L)
3661+ PEMTTEMP(L) = PJET(I,L)
3662+ PRECTEMP(L) = PINC(2,L)
3663+ ENDDO
3664+ PT1 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(1),
3665+ & JETFLAVOUR(I), -1, -1)
3666+ PTMIN = MIN( PTMIN, PT1 )
3667+C Check separation to second incoming particle
3668+ DO L=0,3
3669+ PRADTEMP(L) = PINC(2,L)
3670+ PEMTTEMP(L) = PJET(I,L)
3671+ PRECTEMP(L) = PINC(1,L)
3672+ ENDDO
3673+ PT2 = RHOPYTHIA(PRADTEMP, PEMTTEMP, PRECTEMP, INCFLAVOUR(2),
3674+ & JETFLAVOUR(I), -1, -1)
3675+ PTMIN = MIN( PTMIN, PT2 )
3676+ ENDIF
3677+
3678+C Compute pythia FSR separation between two jets,
3679+C without any knowledge of colour connections
3680+ DO J=1,NJETS
3681+ DO K=1,NJETS
3682+ IF ( I .NE. J .AND. I .NE. K .AND. J .NE. K ) THEN
3683+
3684+C Check separation between final partons i and j, with k as spectator
3685+ DO L=0,3
3686+ PRADTEMP(L) = PJET(J,L)
3687+ PEMTTEMP(L) = PJET(I,L)
3688+ PRECTEMP(L) = PJET(K,L)
3689+ ENDDO
3690+
3691+ TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
3692+ & JETFLAVOUR(J), JETFLAVOUR(I), 1, 1)
3693+C Only SM-like emissions off the beam are possible, no additional
3694+C BSM particles will be produced as as shower emissions.
3695+ IF ( ABS(JETFLAVOUR(I)) .LT. 30 ) THEN
3696+ PTMIN = MIN(PTMIN, TEMP);
3697+ ENDIF
3698+
3699+ TEMP = RHOPYTHIA( PEMTTEMP, PRADTEMP, PRECTEMP,
3700+ & JETFLAVOUR(I), JETFLAVOUR(J), 1, 1)
3701+C Only SM-like emissions off the beam are possible, no additional
3702+C BSM particles will be produced as as shower emissions.
3703+ IF ( ABS(JETFLAVOUR(J)) .LT. 30 ) THEN
3704+ PTMIN = MIN(PTMIN, TEMP);
3705+ ENDIF
3706+
3707+ ENDIF
3708+
3709+ ENDDO ! LOOP OVER NJET
3710+ ENDDO ! LOOP OVER NJET
3711+
3712+C Compute pythia FSR separation between two jets, with initial spectator
3713+ IF ( (LPP(1).NE.0) .OR. (LPP(2).NE.0)) THEN
3714+ DO J=1,NJETS
3715+
3716+C BSM particles can only be radiators, and will not be produced
3717+C as shower emissions.
3718+ IF ( ABS(JETFLAVOUR(I)) .GT. 1000000 ) THEN
3719+ EXIT
3720+ ENDIF
3721+
3722+C Allow both initial partons as recoiler
3723+ IF ( I .NE. J ) THEN
3724+
3725+C Check with first initial as recoiler
3726+ DO L=0,3
3727+ PRADTEMP(L) = PJET(J,L)
3728+ PEMTTEMP(L) = PJET(I,L)
3729+ PRECTEMP(L) = PINC(1,L)
3730+ ENDDO
3731+ TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
3732+ & JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
3733+ IF ( LPP(1) .NE. 0 ) THEN
3734+ PTMIN = MIN(PTMIN, TEMP);
3735+ ENDIF
3736+ DO L=0,3
3737+ PRADTEMP(L) = PJET(J,L)
3738+ PEMTTEMP(L) = PJET(I,L)
3739+ PRECTEMP(L) = PINC(2,L)
3740+ ENDDO
3741+ TEMP = RHOPYTHIA( PRADTEMP, PEMTTEMP, PRECTEMP,
3742+ & JETFLAVOUR(J), JETFLAVOUR(I), 1, -1);
3743+ IF ( LPP(2) .NE. 0 ) THEN
3744+ PTMIN = MIN(PTMIN, TEMP);
3745+ ENDIF
3746+ ENDIF
3747+ ENDDO ! LOOP OVER NJET
3748+ ENDIF
3749+
3750+ PTMINSAVE = MIN(PTMIN,PTMINSAVE)
3751+
3752+ ENDDO ! LOOP OVER NJET
3753+
3754+C CHECK COMPATIBILITY WITH CUT
3755+ IF (PTMINSAVE .LT. PT_LUND) THEN
3756+ PASSCUTS = .FALSE.
3757+ RETURN
3758+ ENDIF
3759+
3760+ ENDIF ! IF NJETS.GT. 0
3761+
3762+ ENDIF ! PT_LUND .GT. 0D0
3763+
3764+C----------------------------------------------------------------------------
3765+C----------------------------------------------------------------------------
3766+
3767+
3768+
3769+
3770+c- check existance of jets if jet cuts are on
3771+ if(njets.lt.1.and.(htjmin.gt.0.or.ptj1min.gt.0).or.
3772+ $ njets.lt.2.and.ptj2min.gt.0.or.
3773+ $ njets.lt.3.and.ptj3min.gt.0.or.
3774+ $ njets.lt.4.and.ptj4min.gt.0)then
3775+ if(debug) write (*,*) i, ' too few jets -> fails'
3776+ passcuts=.false.
3777+ return
3778+ endif
3779+
3780+c - sort jet pts
3781+ do i=1,njets-1
3782+ do j=i+1,njets
3783+ if(ptjet(j).gt.ptjet(i)) then
3784+ temp=ptjet(i)
3785+ ptjet(i)=ptjet(j)
3786+ ptjet(j)=temp
3787+ endif
3788+ enddo
3789+ enddo
3790+ if(debug) write (*,*) 'ordered ',njets,' ',ptjet
3791+c
3792+c Use "and" or "or" prescriptions
3793+c
3794+ inclht=0
3795+
3796+ if(njets.gt.0) then
3797+
3798+ notgood=.not.jetor
3799+ if(debug) write (*,*) 'jetor :',jetor
3800+ if(debug) write (*,*) '0',notgood
3801+
3802+ do i=1,min(njets,4)
3803+ if(debug) write (*,*) i,ptjet(i), ' ',ptjmin4(i),
3804+ $ ':',ptjmax4(i)
3805+ if(jetor) then
3806+c--- if one of the jets does not pass, the event is rejected
3807+ notgood=notgood.or.(ptjmax4(i).ge.0d0 .and.
3808+ $ ptjet(i).gt.ptjmax4(i)) .or.
3809+ $ (ptjet(i).lt.ptjmin4(i))
3810+ if(debug) write (*,*) i,' notgood total:', notgood
3811+ else
3812+c--- all cuts must fail to reject the event
3813+ notgood=notgood.and.(ptjmax4(i).ge.0d0 .and.
3814+ $ ptjet(i).gt.ptjmax4(i) .or.
3815+ $ (ptjet(i).lt.ptjmin4(i)))
3816+ if(debug) write (*,*) i,' notgood total:', notgood
3817+ endif
3818+ enddo
3819+
3820+
3821+ if (notgood) then
3822+ if(debug) write (*,*) i, ' multiple pt -> fails'
3823+ passcuts=.false.
3824+ return
3825+ endif
3826+
3827+c---------------------------
3828+c Ht cuts
3829+C---------------------------
3830+ htj=ptjet(1)
3831+
3832+ do i=2,njets
3833+ htj=htj+ptjet(i)
3834+ if(debug) write (*,*) i, 'htj ',htj
3835+ if(debug.and.i.le.4) write (*,*) 'htmin ',i,' ', htjmin4(i),':',htjmax4(i)
3836+ if(i.le.4)then
3837+ if(htj.lt.htjmin4(i) .or.
3838+ $ htjmax4(i).ge.0d0.and.htj.gt.htjmax4(i)) then
3839+ if(debug) write (*,*) i, ' ht -> fails'
3840+ passcuts=.false.
3841+ return
3842+ endif
3843+ endif
3844+ enddo
3845+
3846+ if(htj.lt.htjmin.or.htjmax.ge.0d0.and.htj.gt.htjmax)then
3847+ if(debug) write (*,*) i, ' htj -> fails'
3848+ passcuts=.false.
3849+ return
3850+ endif
3851+
3852+ inclht=htj
3853+
3854+ endif !if there are jets
3855+
3856+ if(nheavyjets.gt.0) then
3857+ do i=1,nheavyjets
3858+ inclht=inclht+ptheavyjet(i)
3859+ enddo
3860+ endif !if there are heavyjets
3861+
3862+ if(inclht.lt.inclHtmin.or.
3863+ $ inclHtmax.ge.0d0.and.inclht.gt.inclHtmax)then
3864+ if(debug) write (*,*) ' inclhtmin=',inclHtmin,' -> fails'
3865+ passcuts=.false.
3866+ return
3867+ endif
3868+
3869+C
3870+C maximal and minimal pt of the leptons sorted by pt
3871+c
3872+ nleptons=0
3873+
3874+ if(ptl1min.gt.0.or.ptl2min.gt.0.or.ptl3min.gt.0.or.ptl4min.gt.0.or.
3875+ $ ptl1max.ge.0d0.or.ptl2max.ge.0d0.or.
3876+ $ ptl3max.ge.0d0.or.ptl4max.ge.0d0) then
3877+
3878+c - fill ptlepton with the pt's of the leptons.
3879+ do i=nincoming+1,nexternal
3880+ if(is_a_l(i)) then
3881+ nleptons=nleptons+1
3882+ ptlepton(nleptons)=pt(p(0,i))
3883+ endif
3884+ enddo
3885+ if(debug) write (*,*) 'not yet ordered ',njets,' ',ptjet
3886+
3887+c - check existance of leptons if lepton cuts are on
3888+ if(nleptons.lt.1.and.ptl1min.gt.0.or.
3889+ $ nleptons.lt.2.and.ptl2min.gt.0.or.
3890+ $ nleptons.lt.3.and.ptl3min.gt.0.or.
3891+ $ nleptons.lt.4.and.ptl4min.gt.0)then
3892+ if(debug) write (*,*) i, ' too few leptons -> fails'
3893+ passcuts=.false.
3894+ return
3895+ endif
3896+
3897+c - sort lepton pts
3898+ do i=1,nleptons-1
3899+ do j=i+1,nleptons
3900+ if(ptlepton(j).gt.ptlepton(i)) then
3901+ temp=ptlepton(i)
3902+ ptlepton(i)=ptlepton(j)
3903+ ptlepton(j)=temp
3904+ endif
3905+ enddo
3906+ enddo
3907+ if(debug) write (*,*) 'ordered ',nleptons,' ',ptlepton
3908+
3909+ if(nleptons.gt.0) then
3910+
3911+ notgood = .false.
3912+ do i=1,min(nleptons,4)
3913+ if(debug) write (*,*) i,ptlepton(i), ' ',ptlmin4(i),':',ptlmax4(i)
3914+c--- if one of the leptons does not pass, the event is rejected
3915+ notgood=notgood.or.
3916+ $ (ptlmax4(i).ge.0d0.and.ptlepton(i).gt.ptlmax4(i)).or.
3917+ $ (ptlepton(i).lt.ptlmin4(i))
3918+ if(debug) write (*,*) i,' notgood total:', notgood
3919+ enddo
3920+
3921+
3922+ if (notgood) then
3923+ if(debug) write (*,*) i, ' multiple pt -> fails'
3924+ passcuts=.false.
3925+ return
3926+ endif
3927+ endif
3928+ endif
3929+C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
3930+C SPECIAL CUTS
3931+C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
3932+
3933+C REQUIRE AT LEAST ONE JET WITH PT>XPTJ
3934+
3935+ IF(xptj.gt.0d0) THEN
3936+ xvar=0
3937+ do i=nincoming+1,nexternal
3938+ if(is_a_j(i)) xvar=max(xvar,pt(p(0,i)))
3939+ enddo
3940+ if (xvar .lt. xptj) then
3941+ passcuts=.false.
3942+ return
3943+ endif
3944+ ENDIF
3945+
3946+C REQUIRE AT LEAST ONE PHOTON WITH PT>XPTA
3947+
3948+ IF(xpta.gt.0d0) THEN
3949+ xvar=0
3950+ do i=nincoming+1,nexternal
3951+ if(is_a_a(i)) xvar=max(xvar,pt(p(0,i)))
3952+ enddo
3953+ if (xvar .lt. xpta) then
3954+ passcuts=.false.
3955+ return
3956+ endif
3957+ ENDIF
3958+
3959+C REQUIRE AT LEAST ONE B WITH PT>XPTB
3960+
3961+ IF(xptb.gt.0d0) THEN
3962+ xvar=0
3963+ do i=nincoming+1,nexternal
3964+ if(is_a_b(i)) xvar=max(xvar,pt(p(0,i)))
3965+ enddo
3966+ if (xvar .lt. xptb) then
3967+ passcuts=.false.
3968+ return
3969+ endif
3970+ ENDIF
3971+
3972+C REQUIRE AT LEAST ONE LEPTON WITH PT>XPTL
3973+
3974+ IF(xptl.gt.0d0) THEN
3975+ xvar=0
3976+ do i=nincoming+1,nexternal
3977+ if(is_a_l(i)) xvar=max(xvar,pt(p(0,i)))
3978+ enddo
3979+ if (xvar .lt. xptl) then
3980+ passcuts=.false.
3981+ if(debug) write (*,*) ' xptl -> fails'
3982+ return
3983+ endif
3984+ ENDIF
3985+C
3986+C WBF CUTS: TWO TYPES
3987+C
3988+C FIRST TYPE: implemented by FM
3989+C
3990+C 1. FIND THE 2 HARDEST JETS
3991+C 2. REQUIRE |RAP(J)|>XETAMIN
3992+C 3. REQUIRE RAP(J1)*ETA(J2)<0
3993+C
3994+C SECOND TYPE : added by Simon de Visscher 1-08-2007
3995+C
3996+C 1. FIND THE 2 HARDEST JETS
3997+C 2. REQUIRE |RAP(J1)-RAP(J2)|>DELTAETA
3998+C 3. REQUIRE RAP(J1)*RAP(J2)<0
3999+C
4000+C
4001+ hardj1=0
4002+ hardj2=0
4003+ ptmax1=0d0
4004+ ptmax2=0d0
4005+
4006+C-- START IF AT LEAST ONE OF THE CUTS IS ACTIVATED
4007+
4008+ IF(XETAMIN.GT.0D0.OR.DELTAETA.GT.0D0) THEN
4009+
4010+C-- FIND THE HARDEST JETS
4011+
4012+ do i=nincoming+1,nexternal
4013+ if(is_a_j(i)) then
4014+c write (*,*) i,pt(p(0,i))
4015+ if(pt(p(0,i)).gt.ptmax1) then
4016+ hardj2=hardj1
4017+ ptmax2=ptmax1
4018+ hardj1=i
4019+ ptmax1=pt(p(0,i))
4020+ elseif(pt(p(0,i)).gt.ptmax2) then
4021+ hardj2=i
4022+ ptmax2=pt(p(0,i))
4023+ endif
4024+c write (*,*) hardj1,hardj2,ptmax1,ptmax2
4025+ endif
4026+ enddo
4027+
4028+ if (hardj2.eq.0) goto 21 ! bypass vbf cut since not enough jets
4029+
4030+C-- NOW APPLY THE CUT I
4031+
4032+ if (abs(rap(p(0,hardj1))) .lt. xetamin
4033+ & .or.abs(rap(p(0,hardj2))) .lt. xetamin
4034+ & .or.rap(p(0,hardj1))*rap(p(0,hardj2)) .gt.0d0) then
4035+ passcuts=.false.
4036+ return
4037+ endif
4038+
4039+
4040+C-- NOW APPLY THE CUT II
4041+
4042+ if (abs(rap(p(0,hardj1))-rap(p(0,hardj2))) .lt. deltaeta) then
4043+ passcuts=.false.
4044+ return
4045+ endif
4046+
4047+c write (*,*) hardj1,hardj2,rap(p(0,hardj1)),rap(p(0,hardj2))
4048+
4049+ ENDIF
4050+
4051+c Begin photon isolation
4052+c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
4053+c Use is made of parton cm frame momenta. If this must be
4054+c changed, pQCD used below must be redefined
4055+c NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE NOTE
4056+c If we do not require a mimimum jet energy, there's no need to apply
4057+c jet clustering and all that.
4058+ if (ptgmin.ne.0d0) then
4059+
4060+c Put all (light) QCD partons in momentum array for jet clustering.
4061+c From the run_card.dat, maxjetflavor defines if b quark should be
4062+c considered here (via the logical variable 'is_a_jet'). nQCD becomes
4063+c the number of (light) QCD partons at the real-emission level (i.e. one
4064+c more than the Born).
4065+
4066+ nQCD=0
4067+ do j=nincoming+1,nexternal
4068+ if (is_a_j(j)) then
4069+ nQCD=nQCD+1
4070+ do i=0,3
4071+ pQCD(i,nQCD)=p(i,j) ! Use C.o.M. frame momenta
4072+ enddo
4073+ endif
4074+ enddo
4075+
4076+ nph=0
4077+ do j=nincoming+1,nexternal
4078+ if (is_a_a(j)) then
4079+ nph=nph+1
4080+ do i=0,3
4081+ pgamma(i,nph)=p(i,j) ! Use C.o.M. frame momenta
4082+ enddo
4083+ endif
4084+ enddo
4085+ if(nph.eq.0) goto 444
4086+
4087+ if(isoEM)then
4088+ nem=nph
4089+ do k=1,nem
4090+ do i=0,3
4091+ pem(i,k)=pgamma(i,k)
4092+ enddo
4093+ enddo
4094+ do j=nincoming+1,nexternal
4095+ if (is_a_l(j)) then
4096+ nem=nem+1
4097+ do i=0,3
4098+ pem(i,nem)=p(i,j) ! Use C.o.M. frame momenta
4099+ enddo
4100+ endif
4101+ enddo
4102+ endif
4103+
4104+ alliso=.true.
4105+
4106+ j=0
4107+ dowhile(j.lt.nph.and.alliso)
4108+c Loop over all photons
4109+ j=j+1
4110+
4111+ ptg=pt(pgamma(0,j))
4112+ if(ptg.lt.ptgmin)then
4113+ passcuts=.false.
4114+ return
4115+ endif
4116+
4117+c Isolate from hadronic energy
4118+ do i=1,nQCD
4119+ drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
4120+ enddo
4121+ call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
4122+ Etsum(0)=0.d0
4123+ nin=0
4124+ do i=1,nQCD
4125+ if(dble(drlist(isorted(i))).le.R0gamma)then
4126+ nin=nin+1
4127+ Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
4128+ endif
4129+ enddo
4130+ do i=1,nin
4131+ alliso=alliso .and.
4132+ # Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
4133+ # R0gamma,xn,epsgamma,ptg)
4134+ enddo
4135+
4136+c Isolate from EM energy
4137+ if(isoEM.and.nem.gt.1)then
4138+ do i=1,nem
4139+ drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
4140+ enddo
4141+ call sortzv(drlist,isorted,nem,ismode,isway,izero)
4142+c First of list must be the photon: check this, and drop it
4143+ if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
4144+ write(*,*)'Error #1 in photon isolation'
4145+ write(*,*)j,isorted(1),drlist(isorted(1))
4146+ stop
4147+ endif
4148+ Etsum(0)=0.d0
4149+ nin=0
4150+ do i=2,nem
4151+ if(dble(drlist(isorted(i))).le.R0gamma)then
4152+ nin=nin+1
4153+ Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
4154+ endif
4155+ enddo
4156+ do i=1,nin
4157+ alliso=alliso .and.
4158+ # Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
4159+ # R0gamma,xn,epsgamma,ptg)
4160+ enddo
4161+
4162+ endif
4163+
4164+c End of loop over photons
4165+ enddo
4166+
4167+ if(.not.alliso)then
4168+ passcuts=.false.
4169+ return
4170+ endif
4171+ endif
4172+
4173+ 444 continue
4174+c End photon isolation
4175+
4176+c
4177+c call the dummy_cuts function to check plugin/user defined cuts
4178+c
4179+
4180+ if(.not.dummy_cuts(P))then
4181+ passcuts=.false.
4182+ return
4183+ endif
4184+
4185+
4186+C...Set couplings if event passed cuts
4187+
4188+ 21 if(.not.fixed_ren_scale) then
4189+ call set_ren_scale(P,scale)
4190+ if(scale.gt.0) G = SQRT(4d0*PI*ALPHAS(scale))
4191+ endif
4192+
4193+ if(.not.fixed_fac_scale1.or..not.fixed_fac_scale2) then
4194+ call set_fac_scale(P,q2fact)
4195+ endif
4196+
4197+c
4198+c Here we cluster event and reset factorization and renormalization
4199+c scales on an event-by-event basis, as well as check xqcut for jets
4200+c
4201+c Note the following condition is the first line of setclscales
4202+
4203+ IF (FIRSTTIME2) THEN
4204+ FIRSTTIME2=.FALSE.
4205+ write(6,*) 'alpha_s for scale ',scale,' is ', G**2/(16d0*atan(1d0))
4206+ ENDIF
4207+
4208+ if(debug) write (*,*) '============================='
4209+ if(debug) write (*,*) ' EVENT PASSED THE CUTS '
4210+ if(debug) write (*,*) '============================='
4211+
4212+ CUTSPASSED=.TRUE.
4213+
4214+ RETURN
4215+ END
4216+
4217+
4218+C
4219+C FUNCTION FOR ISOLATION
4220+C
4221+
4222+ function iso_getdrv40(p1,p2)
4223+ implicit none
4224+ real*8 iso_getdrv40,p1(0:3),p2(0:3)
4225+ real*8 iso_getdr
4226+c
4227+ iso_getdrv40=iso_getdr(p1(0),p1(1),p1(2),p1(3),
4228+ # p2(0),p2(1),p2(2),p2(3))
4229+ return
4230+ end
4231+
4232+
4233+ function iso_getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
4234+ implicit none
4235+ real*8 iso_getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
4236+ # iso_getpseudorap,iso_getdelphi
4237+c
4238+ deta=iso_getpseudorap(en1,ptx1,pty1,pl1)-
4239+ # iso_getpseudorap(en2,ptx2,pty2,pl2)
4240+ dphi=iso_getdelphi(ptx1,pty1,ptx2,pty2)
4241+ iso_getdr=sqrt(dphi**2+deta**2)
4242+ return
4243+ end
4244+
4245+
4246+ function iso_getpseudorap(en,ptx,pty,pl)
4247+ implicit none
4248+ real*8 iso_getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
4249+ parameter (tiny=1.d-5)
4250+c
4251+ pt=sqrt(ptx**2+pty**2)
4252+ if(pt.lt.tiny.and.abs(pl).lt.tiny)then
4253+ eta=sign(1.d0,pl)*1.d8
4254+ else
4255+ th=atan2(pt,pl)
4256+ eta=-log(tan(th/2.d0))
4257+ endif
4258+ iso_getpseudorap=eta
4259+ return
4260+ end
4261+
4262+
4263+ function iso_getdelphi(ptx1,pty1,ptx2,pty2)
4264+ implicit none
4265+ real*8 iso_getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
4266+ parameter (tiny=1.d-5)
4267+c
4268+ pt1=sqrt(ptx1**2+pty1**2)
4269+ pt2=sqrt(ptx2**2+pty2**2)
4270+ if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
4271+ tmp=ptx1*ptx2+pty1*pty2
4272+ tmp=tmp/(pt1*pt2)
4273+ if(abs(tmp).gt.1.d0+tiny)then
4274+ write(*,*)'Cosine larger than 1'
4275+ stop
4276+ elseif(abs(tmp).ge.1.d0)then
4277+ tmp=sign(1.d0,tmp)
4278+ endif
4279+ tmp=acos(tmp)
4280+ else
4281+ tmp=1.d8
4282+ endif
4283+ iso_getdelphi=tmp
4284+ return
4285+ end
4286+
4287+ function chi_gamma_iso(dr,R0,xn,epsgamma,pTgamma)
4288+c Eq.(3.4) of Phys.Lett. B429 (1998) 369-374 [hep-ph/9801442]
4289+ implicit none
4290+ real*8 chi_gamma_iso,dr,R0,xn,epsgamma,pTgamma
4291+ real*8 tmp,axn
4292+c
4293+ axn=abs(xn)
4294+ tmp=epsgamma*pTgamma
4295+ if(axn.ne.0.d0)then
4296+ tmp=tmp*( (1-cos(dr))/(1-cos(R0)) )**axn
4297+ endif
4298+ chi_gamma_iso=tmp
4299+ return
4300+ end
4301+
4302+
4303+*
4304+* $Id: sortzv.F,v 1.1.1.1 1996/02/15 17:49:50 mclareni Exp $
4305+*
4306+* $Log: sortzv.F,v $
4307+* Revision 1.1.1.1 1996/02/15 17:49:50 mclareni
4308+* Kernlib
4309+*
4310+*
4311+c$$$#include "kerngen/pilot.h"
4312+ SUBROUTINE SORTZV (A,INDEX,N1,MODE,NWAY,NSORT)
4313+C
4314+C CERN PROGLIB# M101 SORTZV .VERSION KERNFOR 3.15 820113
4315+C ORIG. 02/10/75
4316+C
4317+ DIMENSION A(N1),INDEX(N1)
4318+C
4319+C
4320+ N = N1
4321+ IF (N.LE.0) RETURN
4322+ IF (NSORT.NE.0) GO TO 2
4323+ DO 1 I=1,N
4324+ 1 INDEX(I)=I
4325+C
4326+ 2 IF (N.EQ.1) RETURN
4327+ IF (MODE) 10,20,30
4328+ 10 STOP 5 ! CALL SORTTI (A,INDEX,N)
4329+ GO TO 40
4330+C
4331+ 20 STOP 5 ! CALL SSORTTC(A,INDEX,N)
4332+ GO TO 40
4333+C
4334+ 30 CALL SORTTF (A,INDEX,N)
4335+C
4336+ 40 IF (NWAY.EQ.0) GO TO 50
4337+ N2 = N/2
4338+ DO 41 I=1,N2
4339+ ISWAP = INDEX(I)
4340+ K = N+1-I
4341+ INDEX(I) = INDEX(K)
4342+ 41 INDEX(K) = ISWAP
4343+ 50 RETURN
4344+ END
4345+* ========================================
4346+ SUBROUTINE SORTTF (A,INDEX,N1)
4347+C
4348+ DIMENSION A(N1),INDEX(N1)
4349+C
4350+ N = N1
4351+ DO 3 I1=2,N
4352+ I3 = I1
4353+ I33 = INDEX(I3)
4354+ AI = A(I33)
4355+ 1 I2 = I3/2
4356+ IF (I2) 3,3,2
4357+ 2 I22 = INDEX(I2)
4358+ IF (AI.LE.A (I22)) GO TO 3
4359+ INDEX (I3) = I22
4360+ I3 = I2
4361+ GO TO 1
4362+ 3 INDEX (I3) = I33
4363+ 4 I3 = INDEX (N)
4364+ INDEX (N) = INDEX (1)
4365+ AI = A(I3)
4366+ N = N-1
4367+ IF (N-1) 12,12,5
4368+ 5 I1 = 1
4369+ 6 I2 = I1 + I1
4370+ IF (I2.LE.N) I22= INDEX(I2)
4371+ IF (I2-N) 7,9,11
4372+ 7 I222 = INDEX (I2+1)
4373+ IF (A(I22)-A(I222)) 8,9,9
4374+ 8 I2 = I2+1
4375+ I22 = I222
4376+ 9 IF (AI-A(I22)) 10,11,11
4377+ 10 INDEX(I1) = I22
4378+ I1 = I2
4379+ GO TO 6
4380+ 11 INDEX (I1) = I3
4381+ GO TO 4
4382+ 12 INDEX (1) = I3
4383+ RETURN
4384+ END
4385+* ========================================
4386+ SUBROUTINE SORTTI (A,INDEX,N1)
4387+C
4388+ INTEGER A,AI
4389+ DIMENSION A(N1),INDEX(N1)
4390+C
4391+ N = N1
4392+ DO 3 I1=2,N
4393+ I3 = I1
4394+ I33 = INDEX(I3)
4395+ AI = A(I33)
4396+ 1 I2 = I3/2
4397+ IF (I2) 3,3,2
4398+ 2 I22 = INDEX(I2)
4399+ IF (AI.LE.A (I22)) GO TO 3
4400+ INDEX (I3) = I22
4401+ I3 = I2
4402+ GO TO 1
4403+ 3 INDEX (I3) = I33
4404+ 4 I3 = INDEX (N)
4405+ INDEX (N) = INDEX (1)
4406+ AI = A(I3)
4407+ N = N-1
4408+ IF (N-1) 12,12,5
4409+ 5 I1 = 1
4410+ 6 I2 = I1 + I1
4411+ IF (I2.LE.N) I22= INDEX(I2)
4412+ IF (I2-N) 7,9,11
4413+ 7 I222 = INDEX (I2+1)
4414+ IF (A(I22)-A(I222)) 8,9,9
4415+ 8 I2 = I2+1
4416+ I22 = I222
4417+ 9 IF (AI-A(I22)) 10,11,11
4418+ 10 INDEX(I1) = I22
4419+ I1 = I2
4420+ GO TO 6
4421+ 11 INDEX (I1) = I3
4422+ GO TO 4
4423+ 12 INDEX (1) = I3
4424+ RETURN
4425+ END
4426+* ========================================
4427+ SUBROUTINE SORTTC (A,INDEX,N1)
4428+C
4429+ INTEGER A,AI
4430+ DIMENSION A(N1),INDEX(N1)
4431+C
4432+ N = N1
4433+ DO 3 I1=2,N
4434+ I3 = I1
4435+ I33 = INDEX(I3)
4436+ AI = A(I33)
4437+ 1 I2 = I3/2
4438+ IF (I2) 3,3,2
4439+ 2 I22 = INDEX(I2)
4440+ IF(ICMPCH(AI,A(I22)))3,3,21
4441+ 21 INDEX (I3) = I22
4442+ I3 = I2
4443+ GO TO 1
4444+ 3 INDEX (I3) = I33
4445+ 4 I3 = INDEX (N)
4446+ INDEX (N) = INDEX (1)
4447+ AI = A(I3)
4448+ N = N-1
4449+ IF (N-1) 12,12,5
4450+ 5 I1 = 1
4451+ 6 I2 = I1 + I1
4452+ IF (I2.LE.N) I22= INDEX(I2)
4453+ IF (I2-N) 7,9,11
4454+ 7 I222 = INDEX (I2+1)
4455+ IF (ICMPCH(A(I22),A(I222))) 8,9,9
4456+ 8 I2 = I2+1
4457+ I22 = I222
4458+ 9 IF (ICMPCH(AI,A(I22))) 10,11,11
4459+ 10 INDEX(I1) = I22
4460+ I1 = I2
4461+ GO TO 6
4462+ 11 INDEX (I1) = I3
4463+ GO TO 4
4464+ 12 INDEX (1) = I3
4465+ RETURN
4466+ END
4467+* ========================================
4468+ FUNCTION ICMPCH(IC1,IC2)
4469+C FUNCTION TO COMPARE TWO 4 CHARACTER EBCDIC STRINGS - IC1,IC2
4470+C ICMPCH=-1 IF HEX VALUE OF IC1 IS LESS THAN IC2
4471+C ICMPCH=0 IF HEX VALUES OF IC1 AND IC2 ARE THE SAME
4472+C ICMPCH=+1 IF HEX VALUES OF IC1 IS GREATER THAN IC2
4473+ I1=IC1
4474+ I2=IC2
4475+ IF(I1.GE.0.AND.I2.GE.0)GOTO 40
4476+ IF(I1.GE.0)GOTO 60
4477+ IF(I2.GE.0)GOTO 80
4478+ I1=-I1
4479+ I2=-I2
4480+ IF(I1-I2)80,70,60
4481+ 40 IF(I1-I2)60,70,80
4482+ 60 ICMPCH=-1
4483+ RETURN
4484+ 70 ICMPCH=0
4485+ RETURN
4486+ 80 ICMPCH=1
4487+ RETURN
4488+ END
4489+
4490+c************************************************************************
4491+c Returns pTLund (i.e. the Pythia8 evolution variable) between two
4492+c particles with momenta prad and pemt with momentum prec as spectator
4493+c************************************************************************
4494+
4495+ DOUBLE PRECISION FUNCTION RHOPYTHIA(PRAD, PEMT, PREC, FLAVRAD,
4496+ & FLAVEMT, RADTYPE, RECTYPE)
4497+
4498+ IMPLICIT NONE
4499+c
4500+c Arguments
4501+c
4502+ DOUBLE PRECISION PRAD(0:3),PEMT(0:3), PREC(0:3)
4503+ INTEGER FLAVRAD, FLAVEMT, RADTYPE,RECTYPE
4504+c
4505+c Local
4506+c
4507+ DOUBLE PRECISION Q(0:3),SUM(0:3), qBR(0:3), qAR(0:3)
4508+ DOUBLE PRECISION Q2, m2Rad, m2Emt, m2Dip, qBR2, qAR2, x1, x2, z
4509+ DOUBLE PRECISION m2RadAft, m2EmtAft, m2Rec, m2RadBef, m2ar, rescale
4510+ DOUBLE PRECISION TEMP, lambda13, k1, k3, m2Final
4511+ DOUBLE PRECISION m0u, m0d, m0c, m0s, m0t, m0b, m0w, m0z, m0x
4512+ DOUBLE PRECISION PRECAFT(0:3)
4513+ INTEGER emtsign
4514+ INTEGER idRadBef
4515+ LOGICAL allowed
4516+
4517+c-----
4518+c Begin Code
4519+c-----
4520+
4521+C Set masses. Currently no way of getting those?
4522+ m0u = 0.0
4523+ m0d = 0.0
4524+ m0c = 1.5
4525+ m0s = 0.0
4526+ m0b = 4.7
4527+ m0t = 172.5
4528+ m0w = 80.4
4529+ m0z = 91.188
4530+ m0x = 400.0
4531+
4532+C Store recoiler momentum (since FI splittings require recoiler
4533+C rescaling)
4534+ PRECAFT(0) = PREC(0)
4535+ PRECAFT(1) = PREC(1)
4536+ PRECAFT(2) = PREC(2)
4537+ PRECAFT(3) = PREC(3)
4538+C Get sign of emitted momentum
4539+ emtsign = 1
4540+ if(radtype .eq. -1) emtsign = -1
4541+
4542+C Get virtuality
4543+ Q(0) = pRad(0) + emtsign*pEmt(0)
4544+ Q(1) = pRad(1) + emtsign*pEmt(1)
4545+ Q(2) = pRad(2) + emtsign*pEmt(2)
4546+ Q(3) = pRad(3) + emtsign*pEmt(3)
4547+ Q2 = emtsign * ( Q(0)**2 - Q(1)**2 - Q(2)**2 - Q(3)**2 );
4548+
4549+C Reset allowed
4550+ allowed = .true.
4551+
4552+C Splitting not possible for negative virtuality.
4553+ if ( Q2 .lt. 0.0 ) allowed = .false.
4554+
4555+C Try to reconstruct flavour of radiator before emission.
4556+ idRadBef = 0
4557+C gluon radiation: idBef = idAft
4558+ if (abs(flavEmt) .eq. 21 .or. abs(flavEmt) .eq. 22 ) idRadBef=flavRad
4559+C final state gluon splitting: idBef = 21
4560+ if (radtype .eq. 1 .and. flavEmt .eq. -flavRad) idRadBef=21
4561+C final state quark -> gluon conversion
4562+ if (radtype .eq. 1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=flavEmt
4563+C initial state gluon splitting: idBef = -idEmt
4564+ if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. 21) idRadBef=-flavEmt
4565+C initial state gluon -> quark conversion
4566+ if (radtype .eq. -1 .and. abs(flavEmt) .lt. 10 .and. flavRad .eq. flavEmt) idRadBef=21
4567+C W-boson radiation
4568+ if (flavEmt .eq. 24) idRadBef = flavRad+1
4569+ if (flavEmt .eq. -24) idRadBef = flavRad-1
4570+
4571+C Set particle masses.
4572+ m2RadAft = 0.0
4573+ m2EmtAft = 0.0
4574+ m2Rec = 0.0
4575+ m2RadBef = 0.0
4576+
4577+ m2RadAft = pRad(0)*pRad(0)-pRad(1)*pRad(1)-pRad(2)*pRad(2)-pRad(3)*pRad(3)
4578+ m2EmtAft = pEmt(0)*pEmt(0)-pEmt(1)*pEmt(1)-pEmt(2)*pEmt(2)-pEmt(3)*pEmt(3)
4579+ m2Rec = pRec(0)*pRec(0)-pRec(1)*pRec(1)-pRec(2)*pRec(2)-pRec(3)*pRec(3)
4580+
4581+ if (m2RadAft .lt. 1d-4) m2RadAft = 0.0
4582+ if (m2EmtAft .lt. 1d-4) m2EmtAft = 0.0
4583+ if (m2Rec .lt. 1d-4) m2Rec = 0.0
4584+
4585+ if (abs(flavRad) .ne. 21 .and. abs(flavRad) .ne. 22 .and.
4586+ & abs(flavEmt) .ne. 24 .and.
4587+ & abs(flavRad) .ne. abs(flavEmt)) then
4588+ m2RadBef = m2RadAft
4589+ else if (abs(flavEmt) .eq. 24) then
4590+ if (idRadBef .ne. 0) then
4591+ if( abs(idRadBef) .eq. 4 ) m2RadBef = m0c**2
4592+ if( abs(idRadBef) .eq. 5 ) m2RadBef = m0b**2
4593+ if( abs(idRadBef) .eq. 6 ) m2RadBef = m0t**2
4594+ if( abs(idRadBef) .eq. 9000001 ) m2RadBef = m0x**2
4595+ endif
4596+ else if (radtype .eq. -1) then
4597+ if (abs(flavRad) .eq. 21 .and. abs(flavEmt) .eq. 21) m2RadBef = m2EmtAft
4598+ endif
4599+
4600+C Calculate dipole mass for final-state radiation.
4601+ m2Final = 0.0
4602+ m2Final = m2Final + (pRad(0) + pRec(0) + pEmt(0))**2
4603+ m2Final = m2Final - (pRad(1) + pRec(1) + pEmt(1))**2
4604+ m2Final = m2Final - (pRad(2) + pRec(2) + pEmt(2))**2
4605+ m2Final = m2Final - (pRad(3) + pRec(3) + pEmt(3))**2
4606+
4607+C Final state splitting not possible for negative dipole mass.
4608+ if ( radtype .eq. 1 .and. m2Final .lt. 0.0 ) allowed = .false.
4609+
4610+C Rescale recoiler for final-intial splittings.
4611+ rescale = 1.0
4612+ if (radtype .eq. 1 .and. rectype .eq. -1) then
4613+ m2ar = m2Final - 2.0*Q2 + 2.0*m2RadBef
4614+ rescale = (1.0 - (Q2 - m2RadBef) / (m2ar-m2RadBef))
4615+ & /(1.0 + (Q2 - m2RadBef) / (m2ar-m2RadBef))
4616+ pRecAft(0) = pRecAft(0)*rescale
4617+ pRecAft(1) = pRecAft(1)*rescale
4618+ pRecAft(2) = pRecAft(2)*rescale
4619+ pRecAft(3) = pRecAft(3)*rescale
4620+ endif
4621+
4622+C Final-initial splitting not possible for negative rescaling.
4623+ if ( rescale .lt. 0.0 ) allowed = .false.
4624+
4625+C Construct dipole momentum for FSR.
4626+ sum(0) = pRad(0) + pRecAft(0) + pEmt(0)
4627+ sum(1) = pRad(1) + pRecAft(1) + pEmt(1)
4628+ sum(2) = pRad(2) + pRecAft(2) + pEmt(2)
4629+ sum(3) = pRad(3) + pRecAft(3) + pEmt(3)
4630+ m2Dip = sum(0)**2 - sum(1)**2 - sum(2)**2 - sum(3)**2
4631+
4632+C Construct 2->3 variables for FSR
4633+ x1 = 2. * ( sum(0)*pRad(0) - sum(1)*pRad(1)
4634+ & - sum(2)*pRad(2) - sum(3)*pRad(3) ) / m2Dip
4635+ x2 = 2. * ( sum(0)*pRecAft(0) - sum(1)*pRecAft(1)
4636+ & - sum(2)*pRecAft(2) - sum(3)*pRecAft(3) ) / m2Dip
4637+
4638+C Final state splitting not possible for ill-defined
4639+C 3-body-variables.
4640+ if ( radtype .eq. 1 .and. x1 .lt. 0.0 ) allowed = .false.
4641+ if ( radtype .eq. 1 .and. x1 .gt. 1.0 ) allowed = .false.
4642+ if ( radtype .eq. 1 .and. x2 .lt. 0.0 ) allowed = .false.
4643+ if ( radtype .eq. 1 .and. x2 .gt. 1.0 ) allowed = .false.
4644+
4645+C Auxiliary variables for massive FSR
4646+ lambda13 = DSQRT( (Q2 - m2RadAft - m2EmtAft )**2 - 4.0 * m2RadAft*m2EmtAft)
4647+ k1 = ( Q2 - lambda13 + (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
4648+ k3 = ( Q2 - lambda13 - (m2EmtAft - m2RadAft ) ) / ( 2.0 * Q2)
4649+
4650+C Construct momenta of dipole before/after splitting for ISR
4651+ qBR(0) = pRad(0) + pRec(0) - pEmt(0)
4652+ qBR(1) = pRad(1) + pRec(1) - pEmt(1)
4653+ qBR(2) = pRad(2) + pRec(2) - pEmt(2)
4654+ qBR(3) = pRad(3) + pRec(3) - pEmt(3)
4655+ qBR2 = qBR(0)**2 - qBR(1)**2 - qBR(2)**2 - qBR(3)**2
4656+
4657+ qAR(0) = pRad(0) + pRec(0)
4658+ qAR(1) = pRad(1) + pRec(1)
4659+ qAR(2) = pRad(2) + pRec(2)
4660+ qAR(3) = pRad(3) + pRec(3)
4661+ qAR2 = qAR(0)**2 - qAR(1)**2 - qAR(2)**2 - qAR(3)**2
4662+
4663+C Calculate z of splitting, different for FSR and ISR
4664+ z = 1.0 / (1.0 - k1 -k3) * ( x1 / (2.0-x2) - k3)
4665+ if(radtype .eq. -1 ) z = qBR2 / qAR2;
4666+
4667+C Splitting not possible for ill-defined energy sharing.
4668+ if ( z .lt. 0.0 .or. z .gt. 1.0 ) allowed = .false.
4669+
4670+C pT^2 = separation * virtuality (corrected with mass for FSR)
4671+ if (radtype .eq. 1) temp = z*(1-z)*(Q2 - m2RadBef)
4672+ if (radtype .eq. -1) temp = (1-z)*Q2
4673+
4674+C Check threshold in ISR
4675+ if (radtype .ne. 1) then
4676+ if ((abs(flavRad) .eq. 4 .or. abs(flavEmt) .eq. 4)
4677+ & .and. dsqrt(temp) .le. 2.0*m0c**2 ) temp = (1.-z)*(Q2+m0c**2)
4678+ if ((abs(flavRad) .eq. 5 .or. abs(flavEmt) .eq. 5)
4679+ & .and. dsqrt(temp) .le. 2.0*m0b**2 ) temp = (1.-z)*(Q2+m0b**2)
4680+ endif
4681+
4682+C Kinematically impossible splittings should not be included in the
4683+C pT definition!
4684+ if( .not. allowed) temp = 1d15
4685+
4686+ if(temp .lt. 0.0) temp = 0.0
4687+
4688+C Return pT
4689+ rhoPythia = dsqrt(temp);
4690+
4691+ RETURN
4692+ END
4693
4694 logical function cut_bw(p)
4695 include 'madweight_param.inc'
4696
4697=== modified file 'tests/input_files/full_sm/particles.py'
4698--- tests/input_files/full_sm/particles.py 2015-10-01 16:00:08 +0000
4699+++ tests/input_files/full_sm/particles.py 2021-08-21 21:19:44 +0000
4700@@ -269,4 +269,4 @@
4701 'propagating': True,
4702 'is_part': True,
4703 'self_antipart': True
4704-}]
4705\ No newline at end of file
4706+}]
4707
4708=== modified file 'tests/input_files/mssm/particles.py'
4709--- tests/input_files/mssm/particles.py 2015-10-01 16:00:08 +0000
4710+++ tests/input_files/mssm/particles.py 2021-08-21 21:19:44 +0000
4711@@ -719,4 +719,4 @@
4712 'propagating': True,
4713 'is_part': True,
4714 'self_antipart': True
4715-}]
4716\ No newline at end of file
4717+}]
4718
4719=== modified file 'tests/input_files/sm/particles.py'
4720--- tests/input_files/sm/particles.py 2015-10-01 16:00:08 +0000
4721+++ tests/input_files/sm/particles.py 2021-08-21 21:19:44 +0000
4722@@ -373,4 +373,4 @@
4723 'apx_decaywidth_err': 0.00,
4724 'decay_amplitudes': {},
4725 '2body_massdiff': 0.00
4726-}]
4727\ No newline at end of file
4728+}]
4729
4730=== modified file 'tests/parallel_tests/decay_comparator.py'
4731--- tests/parallel_tests/decay_comparator.py 2021-02-14 22:27:00 +0000
4732+++ tests/parallel_tests/decay_comparator.py 2021-08-21 21:19:44 +0000
4733@@ -531,4 +531,4 @@
4734 start = time.time()
4735 print('comparing decay for %s' % name)
4736 self.assertEqual('True', decay_framework.has_same_decay(name))
4737- print('done in %s s' % (time.time() - start))
4738\ No newline at end of file
4739+ print('done in %s s' % (time.time() - start))
4740
4741=== modified file 'tests/parallel_tests/me_comparator.py'
4742--- tests/parallel_tests/me_comparator.py 2020-02-14 10:21:44 +0000
4743+++ tests/parallel_tests/me_comparator.py 2021-08-21 21:19:44 +0000
4744@@ -161,7 +161,7 @@
4745 raise IOError("Path %s for test already exist" % \
4746 str(os.path.join(mg4_path, temp_dir)))
4747
4748- shutil.copytree(os.path.join(mg4_path, 'Template'),
4749+ misc.copytree(os.path.join(mg4_path, 'Template'),
4750 os.path.join(mg4_path, temp_dir))
4751
4752 self.temp_dir_name = temp_dir
4753
4754=== modified file 'tests/parallel_tests/test_ML5MSSMQCD.py'
4755--- tests/parallel_tests/test_ML5MSSMQCD.py 2019-04-17 14:39:47 +0000
4756+++ tests/parallel_tests/test_ML5MSSMQCD.py 2021-08-21 21:19:44 +0000
4757@@ -71,7 +71,7 @@
4758 def setUp(self):
4759 """ Copy the model form the test input_files to the default directory."""
4760 if not os.path.exists(os.path.join(_mg5_path,'models','loop_MSSM')):
4761- shutil.copytree(os.path.join(_mg5_path,'tests','input_files','loop_MSSM'),
4762+ misc.copytree(os.path.join(_mg5_path,'tests','input_files','loop_MSSM'),
4763 os.path.join(_mg5_path,'models','loop_MSSM'))
4764
4765 # The tests below probe one quite long process at a time individually, so
4766
4767=== modified file 'tests/test_manager.py'
4768--- tests/test_manager.py 2021-05-25 08:29:08 +0000
4769+++ tests/test_manager.py 2021-08-21 21:19:44 +0000
4770@@ -361,7 +361,7 @@
4771 if update and path.isdir(_hc_comparison_files):
4772 if path.isdir(hc_comparison_files_BackUp):
4773 shutil.rmtree(hc_comparison_files_BackUp)
4774- shutil.copytree(_hc_comparison_files,hc_comparison_files_BackUp)
4775+ misc.copytree(_hc_comparison_files,hc_comparison_files_BackUp)
4776
4777 IOTestManager.testFolders_filter = arg.split('/')[0].split('&')
4778 IOTestManager.testNames_filter = arg.split('/')[1].split('&')
4779@@ -498,7 +498,7 @@
4780 else:
4781 if path.isdir(hc_comparison_files_BackUp):
4782 shutil.rmtree(_hc_comparison_files)
4783- shutil.copytree(hc_comparison_files_BackUp,_hc_comparison_files)
4784+ misc.copytree(hc_comparison_files_BackUp,_hc_comparison_files)
4785 print(colored%(32,"INFO:: No modifications applied."))
4786 else:
4787 print(colored%(31,
4788
4789=== modified file 'tests/unit_tests/__init__.py'
4790--- tests/unit_tests/__init__.py 2019-04-17 14:39:47 +0000
4791+++ tests/unit_tests/__init__.py 2021-08-21 21:19:44 +0000
4792@@ -78,4 +78,4 @@
4793 # tests.NBTEST += 1
4794 # unittest.TestCase.assertFalse(self, *arg, **opt)
4795
4796-
4797\ No newline at end of file
4798+
4799
4800=== modified file 'tests/unit_tests/loop/test_loop_exporters.py'
4801--- tests/unit_tests/loop/test_loop_exporters.py 2019-04-17 14:39:47 +0000
4802+++ tests/unit_tests/loop/test_loop_exporters.py 2021-08-21 21:19:44 +0000
4803@@ -319,4 +319,4 @@
4804 glob.glob(pjoin(self.IOpath,'helas_calls_*.f'))+\
4805 glob.glob(pjoin(self.IOpath,'loop_CT_calls*.f')):
4806 base_name = '.'.join(file.split('.')[:-1])+'_%s.'+file.split('.')[-1]
4807- shutil.move(file,base_name%name)
4808\ No newline at end of file
4809+ shutil.move(file,base_name%name)
4810
4811=== modified file 'tests/unit_tests/madweight/test_permutation.py'
4812--- tests/unit_tests/madweight/test_permutation.py 2019-04-17 14:39:47 +0000
4813+++ tests/unit_tests/madweight/test_permutation.py 2021-08-21 21:19:44 +0000
4814@@ -121,4 +121,4 @@
4815 self.assertEqual(len(output), 24)
4816 solution = [[1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 3, 4, 5, 6, 8, 7], [1, 2, 3, 4, 5, 7, 6, 8], [1, 2, 3, 4, 5, 7, 8, 6], [1, 2, 3, 4, 5, 8, 6, 7], [1, 2, 3, 4, 5, 8, 7, 6], [1, 6, 3, 4, 5, 2, 7, 8], [1, 6, 3, 4, 5, 2, 8, 7], [1, 6, 3, 4, 5, 7, 2, 8], [1, 6, 3, 4, 5, 7, 8, 2], [1, 6, 3, 4, 5, 8, 2, 7], [1, 6, 3, 4, 5, 8, 7, 2], [1, 7, 3, 4, 5, 2, 6, 8], [1, 7, 3, 4, 5, 2, 8, 6], [1, 7, 3, 4, 5, 6, 2, 8], [1, 7, 3, 4, 5, 6, 8, 2], [1, 7, 3, 4, 5, 8, 2, 6], [1, 7, 3, 4, 5, 8, 6, 2], [1, 8, 3, 4, 5, 2, 6, 7], [1, 8, 3, 4, 5, 2, 7, 6], [1, 8, 3, 4, 5, 6, 2, 7], [1, 8, 3, 4, 5, 6, 7, 2], [1, 8, 3, 4, 5, 7, 2, 6], [1, 8, 3, 4, 5, 7, 6, 2]]
4817 for sol in solution:
4818- self.assertTrue(sol in output)
4819\ No newline at end of file
4820+ self.assertTrue(sol in output)
4821
4822=== modified file 'tests/unit_tests/various/test_banner.py'
4823--- tests/unit_tests/various/test_banner.py 2021-02-12 15:28:58 +0000
4824+++ tests/unit_tests/various/test_banner.py 2021-08-21 21:19:44 +0000
4825@@ -23,6 +23,7 @@
4826 import models
4827 import six
4828 StringIO = six
4829+import sys
4830 from madgraph import MG5DIR
4831
4832
4833@@ -610,6 +611,37 @@
4834 run_card3.write(fsock2)
4835 fsock2.close()
4836 self.assertEqual(open(fsock.name).read(), open(fsock2.name).read())
4837+
4838+ def test_check_valid_LO(self):
4839+ """ensure that some handling are done correctly"""
4840+
4841+ run_card = bannermod.RunCardLO()
4842+ run_card['dsqrt_q2fact1'] = 10
4843+ run_card['dsqrt_q2fact2'] = 20
4844+
4845+ # check that if fixed_fac_scale is on False and lpp1=1 fixed_fac_scale1 is False but if lpp=2 then fixed_fac_scale1 is True
4846+ run_card.set('fixed_fac_scale', False, user=True)
4847+ run_card.set('lpp1', 2, user=True)
4848+ run_card.check_validity()
4849+ self.assertEqual(run_card['fixed_fac_scale1'], True)
4850+ self.assertEqual(run_card['fixed_fac_scale2'], False)
4851+ run_card.set('lpp1', 1, user=True)
4852+ run_card.check_validity()
4853+ self.assertEqual(run_card['fixed_fac_scale1'], False)
4854+ self.assertEqual(run_card['fixed_fac_scale2'], False)
4855+
4856+ # check that for elastisc a a collision we force to use fixed_fac_scale1/2
4857+ run_card.set('lpp1', 2, user=True)
4858+ run_card.set('lpp2', 2, user=True)
4859+ with self.assertRaises(bannermod.InvalidRunCard):
4860+ run_card.check_validity()
4861+ run_card.set('fixed_fac_scale1', False, user=True)
4862+ run_card.set('fixed_fac_scale2', False, user=True)
4863+ run_card.check_validity() # no crashing anymore
4864+ self.assertEqual(run_card['fixed_fac_scale1'], False)
4865+ self.assertEqual(run_card['fixed_fac_scale2'], False)
4866+
4867+
4868
4869
4870 MadLoopParam = bannermod.MadLoopParam
4871
4872=== modified file 'tests/unit_tests/various/test_cmd.py'
4873--- tests/unit_tests/various/test_cmd.py 2021-02-12 15:28:58 +0000
4874+++ tests/unit_tests/various/test_cmd.py 2021-08-21 21:19:44 +0000
4875@@ -58,4 +58,4 @@
4876 are different. This probably fine but please check it before any release."""
4877 if text1 != text2:
4878 print(warning)
4879-
4880\ No newline at end of file
4881+
4882
4883=== modified file 'tests/unit_tests/various/test_misc.py'
4884--- tests/unit_tests/various/test_misc.py 2019-04-17 14:39:47 +0000
4885+++ tests/unit_tests/various/test_misc.py 2021-08-21 21:19:44 +0000
4886@@ -119,4 +119,4 @@
4887 self.assertFalse(eq(1e-8, 0, zero_limit=False))
4888 self.assertFalse(eq(1e-7, 0, zero_limit=False))
4889 self.assertFalse(eq(1e-6, 0, zero_limit=False))
4890- self.assertFalse(eq(1e-1, 0, zero_limit=False))
4891\ No newline at end of file
4892+ self.assertFalse(eq(1e-1, 0, zero_limit=False))
4893
4894=== modified file 'tests/unit_tests/various/test_rambo.py'
4895--- tests/unit_tests/various/test_rambo.py 2019-04-17 14:39:47 +0000
4896+++ tests/unit_tests/various/test_rambo.py 2021-08-21 21:19:44 +0000
4897@@ -547,4 +547,4 @@
4898 values = wavefunctions.vxxxxx(P, M, NHEL, NRST)
4899 self.assertEqual(len(values), len(solution))
4900 for i in range(len(values)):
4901- self.assertAlmostEqual(values[i], solution[i])
4902\ No newline at end of file
4903+ self.assertAlmostEqual(values[i], solution[i])

Subscribers

People subscribed via source and target branches

to all changes: