2
C**************************************************************************
4
C**************************************************************************
10
include 'nexternal.inc'
18
parameter (ZERO = 0d0)
20
parameter( Pi = 3.14159265358979323846d0 )
27
integer icollider,detail_level
28
logical do_cuts(nexternal)
35
double precision pmass(nexternal)
37
double precision spole(maxinvar),swidth(maxinvar),bwjac
38
common/to_brietwigner/spole ,swidth ,bwjac
40
double precision etmin(nincoming+1:nexternal)
41
double precision etamax(nincoming+1:nexternal)
42
double precision emin(nincoming+1:nexternal)
43
double precision r2min(nincoming+1:nexternal,nincoming+1:nexternal)
44
double precision s_min(nexternal,nexternal)
45
double precision etmax(nincoming+1:nexternal)
46
double precision etamin(nincoming+1:nexternal)
47
double precision emax(nincoming+1:nexternal)
48
double precision r2max(nincoming+1:nexternal,nincoming+1:nexternal)
49
double precision s_max(nexternal,nexternal)
50
common/to_cuts/ etmin, emin, etamax, r2min, s_min,
51
$ etmax, emax, etamin, r2max, s_max
53
double precision ptjmin4(4),ptjmax4(4),htjmin4(2:4),htjmax4(2:4)
55
common/to_jet_cuts/ ptjmin4,ptjmax4,htjmin4,htjmax4,jetor
58
c les houches accord stuff to identify neutrinos
61
parameter (maxflow=999)
62
integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
63
& icolup(2,nexternal,maxflow)
64
c include 'leshouche.inc'
65
common /c_leshouche_inc/idup,mothup,icolup
67
LOGICAL IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
68
LOGICAL IS_A_PH(NEXTERNAL)
69
COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
74
parameter (maxpara=100)
75
character*20 param(maxpara),value(maxpara)
78
c setup masses for the final-state particles (fills the /to_mass/ common block)
86
c No pdfs for decay processes - set here since here we know the nincoming
87
c Also set stot here, and use mass of incoming particle for ren scale
89
if(nincoming.eq.1)then
95
fixed_ren_scale=.true.
97
c-check consistency of maxjetflavor in the run_card and with Nf
98
c specified in coupl.inc
99
if (maxjetflavor.gt.int(Nf)) then
100
write(*,*) "WARNING: the value of maxjetflavor specified in"//
101
1 " the run_card is inconsistent with the number of light"//
102
1 " flavours in the model."
103
write(*,*)"Hence it will be set to ", int(Nf)
104
maxjetflavor = int(Nf)
109
do i=nincoming+1,nexternal
116
if (abs(idup(i,1)).le.min(maxjetflavor,5)) then
119
if (abs(idup(i,1)).eq.21) is_a_j(i)=.true. ! gluon is a jet
122
if (idup(i,1).eq.11) is_a_lm(i)=.true. ! e-
123
if (idup(i,1).eq.13) is_a_lm(i)=.true. ! mu-
124
if (idup(i,1).eq.15) is_a_lm(i)=.true. ! ta-
125
if (idup(i,1).eq.-11) is_a_lp(i)=.true. ! e-
126
if (idup(i,1).eq.-13) is_a_lp(i)=.true. ! mu-
127
if (idup(i,1).eq.-15) is_a_lp(i)=.true. ! ta-
130
if (idup(i,1).eq.22) is_a_ph(i)=.true. ! photon
137
subroutine set_tau_min()
138
c Sets the lower bound for tau=x1*x2, using information on particle
139
c masses and on the jet minimum pt, as entered in run_card.dat,
142
double precision zero,vtiny
143
parameter (zero=0.d0,vtiny=1d-8)
147
include 'nexternal.inc'
149
include 'nFKSconfigs.inc'
150
include "fks_info.inc"
151
LOGICAL IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
152
LOGICAL IS_A_PH(NEXTERNAL)
153
COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
155
double precision pmass(-nexternal:0,lmaxconfigs)
156
double precision pwidth(-nexternal:0,lmaxconfigs)
157
integer pow(-nexternal:0,lmaxconfigs)
158
integer itree(2,-max_branch:-1),iconfig
159
common /to_itree/itree,iconfig
161
COMMON/C_NFKSPROCESS/NFKSPROCESS
162
double precision taumin(fks_configs),taumin_s(fks_configs)
163
& ,taumin_j(fks_configs),stot
164
save taumin,taumin_s,taumin_j,stot
165
integer i,j,k,d1,d2,iFKS
166
double precision xm(-nexternal:nexternal),xm1,xm2,xmi
167
double precision xw(-nexternal:nexternal),xw1,xw2,xwi
169
double precision tau_Born_lower_bound,tau_lower_bound_resonance
171
common/ctau_lower_bound/tau_Born_lower_bound
172
& ,tau_lower_bound_resonance,tau_lower_bound
174
double precision mass_min(-nexternal:nexternal)
175
integer cBW_FKS_level_max(fks_configs),
176
& cBW_FKS(fks_configs,-nexternal:-1),
177
& cBW_FKS_level(fks_configs,-nexternal:-1)
178
double precision cBW_FKS_mass(fks_configs,-nexternal:-1,-1:1),
179
& cBW_FKS_width(fks_configs,-nexternal:-1,-1:1)
180
save cBW_FKS_level_max,cBW_FKS,cBW_FKS_level,cBW_FKS_mass
182
integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
183
double precision cBW_mass(-nexternal:-1,-1:1),
184
& cBW_width(-nexternal:-1,-1:1)
185
common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
187
double precision s_mass(-nexternal:-1)
188
$ ,s_mass_FKS(fks_configs,-nexternal:nexternal)
190
common/to_phase_space_s_channel/s_mass
192
real*8 emass(nexternal)
195
data firsttime /.true./
196
include "born_props.inc"
198
if(.not.IS_A_J(NEXTERNAL))then
199
write(*,*)'Fatal error in set_tau_min'
202
c The following assumes that light QCD particles are at the end of the
203
c list. Exclude one of them to set tau bound at the Born level This
204
c sets a hard cut in the minimal shat of the Born phase-space
207
c The contribution from ptj should be treated only as a 'soft lower
208
c bound' if j_fks is initial state: the real-emission i_fks parton is
209
c not necessarily the softest. Therefore, it could be that even though
210
c the Born does not have enough energy to pass the cuts set by ptj, the
214
do iFKS=1,fks_configs
219
do i=nincoming+1,nexternal
220
c Add the minimal jet pTs to tau
221
if(IS_A_J(i) .and. i.ne.nexternal)then
222
if (abs(emass(i)).gt.vtiny) then
223
write (*,*) 'Error in set_tau_min in setcuts.f:'
224
write (*,*) 'mass of a jet should be zero',i
228
if (j_fks.gt.nincoming .and. j_fks.lt.nexternal) then
229
taumin(iFKS)=taumin(iFKS)+ptj
230
taumin_s(iFKS)=taumin_s(iFKS)+ptj
231
taumin_j(iFKS)=taumin_j(iFKS)+ptj
232
elseif (j_fks.ge.1 .and. j_fks.le.nincoming) then
233
taumin_s(iFKS)=taumin_s(iFKS)+ptj
234
taumin_j(iFKS)=taumin_j(iFKS)+ptj
235
elseif (j_fks.eq.nexternal) then
237
& 'ERROR, j_fks cannot be the final parton'
241
write (*,*) 'ERROR, j_fks not correctly defined'
246
c Add the minimal photon pTs to tau
247
elseif(IS_A_PH(i))then
248
if (abs(emass(i)).gt.vtiny) then
249
write (*,*) 'Error in set_tau_min in setcuts.f:'
250
write (*,*) 'mass of a photon should be zero',i
254
if (j_fks.gt.nincoming)
255
& taumin(iFKS)=taumin(iFKS)+ptgmin
256
taumin_s(iFKS)=taumin_s(iFKS)+ptgmin
257
taumin_j(iFKS)=taumin_j(iFKS)+ptgmin
258
xm(i)=emass(i)+ptgmin
259
elseif (is_a_lp(i)) then
260
c Add the postively charged lepton pTs to tau
261
taumin(iFKS)=taumin(iFKS)+emass(i)
262
if (j_fks.gt.nincoming)
263
& taumin(iFKS)=taumin(iFKS)+ptl
264
taumin_s(iFKS)=taumin_s(iFKS)+emass(i)+ptl
265
taumin_j(iFKS)=taumin_j(iFKS)+emass(i)+ptl
267
c Add the lepton invariant mass to tau if there is at least another
268
c lepton of opposite charge. (Only add half of it, i.e. 'the part
269
c contributing from this lepton'). Remove possible overcounting with the
271
do j=nincoming+1,nexternal
273
if (j_fks.gt.nincoming)
274
& taumin(iFKS)= taumin(iFKS)-ptl-emass(i) +
275
& max(mll/2d0,ptl+emass(i))
276
taumin_s(iFKS) = taumin_s(iFKS)-ptl-emass(i)
277
$ + max(mll/2d0,ptl+emass(i))
278
taumin_j(iFKS) = taumin_j(iFKS)-ptl-emass(i)
279
$ + max(mll/2d0,ptl+emass(i))
280
xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
285
elseif (is_a_lm(i)) then
286
c Add the negatively charged lepton pTs to tau
287
taumin(iFKS)=taumin(iFKS)+emass(i)
288
if (j_fks.gt.nincoming)
289
& taumin(iFKS)=taumin(iFKS)+ptl
290
taumin_s(iFKS)=taumin_s(iFKS)+emass(i)+ptl
291
taumin_j(iFKS)=taumin_j(iFKS)+emass(i)+ptl
293
c Add the lepton invariant mass to tau if there is at least another
294
c lepton of opposite charge. (Only add half of it, i.e. 'the part
295
c contributing from this lepton'). Remove possible overcounting with the
297
do j=nincoming+1,nexternal
299
if (j_fks.gt.nincoming)
300
& taumin(iFKS) = taumin(iFKS)-ptl-emass(i) +
301
& max(mll/2d0,ptl+emass(i))
302
taumin_s(iFKS) = taumin_s(iFKS)-ptl-emass(i)
303
$ + max(mll/2d0,ptl+emass(i))
304
taumin_j(iFKS) = taumin_j(iFKS)-ptl-emass(i)
305
$ + max(mll/2d0,ptl+emass(i))
306
xm(i)=xm(i)-ptl-emass(i)+max(mll/2d0,ptl
312
taumin(iFKS)=taumin(iFKS)+emass(i)
313
taumin_s(iFKS)=taumin_s(iFKS)+emass(i)
314
taumin_j(iFKS)=taumin_j(iFKS)+emass(i)
319
stot = 4d0*ebeam(1)*ebeam(2)
320
tau_Born_lower_bound=taumin(iFKS)**2/stot
321
tau_lower_bound=taumin_j(iFKS)**2/stot
323
c Also find the minimum lower bound if all internal s-channel particles
326
do i=-1,-(nexternal-3),-1 ! All propagators
327
if ( itree(1,i) .eq. 1 .or. itree(1,i) .eq. 2 ) tsign=1
328
if (tsign.eq.-1) then ! Only s-channels
331
c If daughter is a jet, we should treat the ptj as a mass. Except if
332
c d1=nexternal, because we check the Born, so final parton should be
333
c skipped. [This is already done above; also for the leptons]
338
c On-shell mass of the intermediate resonance
340
c Width of the intermediate resonance
341
xwi=pwidth(i,iconfig)
342
c Set the intermediate mass equal to the max of its actual mass and
343
c the sum of the masses of the two daugters.
344
if (xmi.gt.xm1+xm2) then
349
xw(i)=xw1+xw2 ! just sum the widths
351
c Add the new mass to the bound. To avoid double counting, we should
352
c subtract the daughters, because they are already included above or in
353
c the previous iteration of the loop
354
taumin_s(iFKS)=taumin_s(iFKS)+xm(i)-xm1-xm2
360
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
361
c Determine the "minimal" s-channel invariant masses
362
do i=nincoming+1,nexternal-1
363
s_mass_FKS(iFKS,i)=xm(i)**2
365
do i=-1,-(nexternal-3),-1 ! All propagators
366
if ( itree(1,i) .eq. 1 .or. itree(1,i) .eq. 2 ) exit ! only s-channels
367
s_mass_FKS(iFKS,i)=(sqrt(s_mass_FKS(iFKS,itree(1,i)))
368
$ +sqrt(s_mass_FKS(iFKS,itree(2,i))))**2
371
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
372
c Determine the conflicting Breit-Wigner's. Note that xm(i) contains the
374
do i=nincoming+1,nexternal-1
375
mass_min(i)=xm(i) ! minimal allowed resonance mass (including masses set by cuts)
377
cBW_FKS_level_max(iFKS)=0
378
do i=-1,-(nexternal-3),-1 ! All propagators
379
cBW_FKS_mass(iFKS,i,1)=0d0
380
cBW_FKS_width(iFKS,i,1)=0d0
381
cBW_FKS_mass(iFKS,i,-1)=0d0
382
cBW_FKS_width(iFKS,i,-1)=0d0
383
if ( itree(1,i) .eq. 1 .or. itree(1,i) .eq. 2 ) exit ! only s-channels
384
mass_min(i)=mass_min(itree(1,i))+mass_min(itree(2,i))
385
if (xm(i).lt.mass_min(i)) then
387
$ 'ERROR in the determination of conflicting BW',i
391
if (pmass(i,iconfig).lt.xm(i)) then
392
c Possible conflict in BW
393
if (pmass(i,iconfig).lt.mass_min(i)) then
394
c Resonance can never go on-shell due to the kinematics of the event
396
cBW_FKS_level(iFKS,i)=0
397
elseif(pmass(i,iconfig).lt.xm(i)) then
398
c Conflicting Breit-Wigner
400
cBW_FKS_level(iFKS,i)=1
401
cBW_FKS_level_max(iFKS)=max(cBW_FKS_level_max(iFKS)
402
$ ,cBW_FKS_level(iFKS,i))
403
c Set here the mass (and width) of the alternative mass; it's the
404
c sum of daughter masses. (3rd argument is '1', because this
405
c alternative mass is LARGER than the resonance mass).
406
cBW_FKS_mass(iFKS,i,1)=xm(i)
407
cBW_FKS_width(iFKS,i,1)=xw(i)
409
c set the daughters also as conflicting (recursively)
411
if (cBW_FKS(iFKS,j).ne.0) then
412
do k=1,2 ! loop over the 2 daughters
413
if (itree(k,j).lt.0) then
414
if(cBW_FKS(iFKS,itree(k,j)).ne.2) then
415
cBW_FKS(iFKS,itree(k,j))=1
416
cBW_FKS_level(iFKS,itree(k,j))=
417
& cBW_FKS_level(iFKS,j)+1
418
cBW_FKS_level_max(iFKS)=
419
$ max(cBW_FKS_level_max(iFKS)
420
$ ,cBW_FKS_level(iFKS,itree(k,j)))
421
c Set here the mass (and width) of the alternative mass; it's the
422
c difference between the mother and the sister masses. (3rd argument
423
c is '-1', because this alternative mass is SMALLER than the
425
cBW_FKS_mass(iFKS,itree(k,j),-1)=
426
& pmass(j,iconfig)-xm(itree(3-k,j)) ! mass difference
427
cBW_FKS_width(iFKS,itree(k,j),-1)=
428
& pwidth(j,iconfig)+xw(itree(3-k,j)) ! sum of widths
435
c Normal Breit-Wigner
439
c Conflicting BW's determined. They are saved in cBW_FKS
440
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
442
c If the lower bound found here is smaller than the hard bound,
443
c simply set the soft bound equal to the hard bound.
445
& max(taumin_j(iFKS),taumin_s(iFKS))
447
c For the bound, we have to square and divide by stot.
448
tau_lower_bound_resonance=taumin_s(iFKS)**2/stot
450
write (*,'(a,i3,a,3(e12.5,x))') 'nFKSprocess:',iFKS
451
& ,'. Absolute lower bound for tau at the Born is'
452
& ,tau_Born_lower_bound,taumin(iFKS),dsqrt(stot)
453
if (j_fks.le.nincoming) then
454
write (*,'(a,i3,a,3(e12.5,x))') 'nFKSprocess:',iFKS
455
& ,'. Lower bound for tau is',tau_lower_bound
456
& ,taumin_j(iFKS),dsqrt(stot)
458
write (*,'(a,i3,a,3(e12.5,x))') 'nFKSprocess:',iFKS
459
& ,'. Lower bound for tau is (taking resonances'/
460
& /' into account)' ,tau_lower_bound_resonance
461
& ,taumin_s(iFKS) ,dsqrt(stot)
464
tau_Born_lower_bound=taumin(nFKSprocess)**2/stot
465
tau_lower_bound=taumin_j(nFKSprocess)**2/stot
466
tau_lower_bound_resonance=taumin_s(nFKSprocess)**2/stot
468
cBW(i)=cBW_FKS(nFKSprocess,i)
469
cBW_level(i)=cBW_FKS_level(nFKSprocess,i)
471
cBW_mass(i,j)=cBW_FKS_mass(nFKSprocess,i,j)
472
cBW_width(i,j)=cBW_FKS_width(nFKSprocess,i,j)
474
s_mass(i)=s_mass_FKS(nFKSprocess,i)
476
cBW_level_max=cBW_FKS_level_max(nFKSprocess)
added
removed
Template/NLO/SubProcesses/setcuts.f
