-
Notifications
You must be signed in to change notification settings - Fork 0
/
djangoh_l.f
6711 lines (6305 loc) · 205 KB
/
djangoh_l.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
c
c***********************************************************************
c *
c djangoh version 1.8 *
c last mod by hs, 20.02.2016 *
c - to correct nucleon scattering *
c modifications by hs, 24.11.97 *
c update for LEPTO version 6.5.1 *
c *
c***********************************************************************
c
chs
chs DJANGO6: LEPTO65 part + new routines of DJANGO (CQP)
chs -----------------------------------------------
c
C*********************************************************************
C List of new routines or modified ones from LEPTO651
C
C List of subprograms in order of appearance, with main purpose
C (S = subroutine, F = function, B = block data)
C *
C B LEPTOD give default values to switches and parameters *
C S DJGINIT initialization (orig.LINIT) *
C S DJGCHC read input values *
C S HSWCUT check if enough W *
C S DJGVAR copy some kinem. var. from HS to Lepto Common blocks *
C S DJGEVT administer event generation (orig.LEPTO) *
C S LEPTOX select process and choose kinematical variables *
C F LKINEM mainly for LWEITS; to calculate some kinem.var. *
C S LWEITS to calculate QCD weights; *
C S DJGBEG add boson,scattered electron and photon from HS *
C S DJGLEV restore event when hadronization has failed *
C S DFRAME in double precision (orig.LFRAME) *
C S DUROBO rotation+boost in DP (and DUDBRB entry) (orig.LUROBO)*
C S LQEV generate quark event *
C S LSHOWR parton shower *
C S LMEPS parton cascades and matrix element *
C S LYREMN parton cascades (remnant) *
C*********************************************************************
ck
ck unmodified routines of LEPTO651
ck
ck (S) LTIMEX
ck (S) LQCDPR
ck (S) LXP
ck (S) LZP
ck (F) LQMCUT
ck (F) DSIGMA
ck (F) DQCD
ck (F) DQCDI
ck (S) LREMH
ck (S) LPRIKT
ck (S) LPRWTS
ck (S) LSCALE
ck (S) LYSPLI
ch (S) LNSTRF
ch (S) LFLAV
ch (S) LSCI
ch (S) LEASWI
ch (S) LECSWI
ch (S) LSMALL
ch
ch the following routines use double precision DUROBO instead of LUROBO
ch (otherwise unchanged)
ck (S) LQGEV
ck (S) LQQBEV
ck (S) LAZIMU
ck (S) LMEPS
ck (S) LYSSPA
ch
ch the following routines originally contained in LEPTO are part of
ch HERACLES 4.6
ch (S) FLTABL modified
ch (S) FLIPOL modified
ch (S) FLINTG unchanged
ch (F) FLQINT unchanged
ch (F) FLGINT unchanged
ch (F) FLTINT unchanged
ch (S) GADAP unchanged
ch (S) LYSTFU modified
ch
ch the following routines originally contained in LEPTO are removed
ch (S) LWBB
ch (S) LSIGMX
ch (S) LXSECT
ch (S) RIWIBD
ch (S) DVNOPT
ch (F) DFUN
ch (F) RIWFUN
ch (F) DCROSS
ch (F) DLOWER
ch (F) DUPPER
ch also removed are the routines from the MINUIT package
ch
ch
ck*******************************************************************
ckc..from L61
chs..updated for L65
BLOCK DATA LEPTOD
C...Give sensible default values to switches and parameters.
COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U
COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17)
COMMON /LFLMIX/ CABIBO(4,4)
COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),OPTW2(4),COMFAC
COMMON /LGRID/ NXX,NWW,XX(31),WW(21),PQG(31,21,3),PQQB(31,21,2),
&QGMAX(31,21,3),QQBMAX(31,21,2),YCUT(31,21),XTOT(31,21),NP
COMMON /FLGRID/ NFX,NFQ,XR(2),QR(2),FLQT(41,16),FLGT(41,16),
&FLMT(41,16)
COMMON /LYPARA/ IPY(80),PYPAR(80),PYVAR(80)
COMMON /LMINUI/ XKIN(4),UKIN(4),WKIN(4),AIN(4),BIN(4),
&MAXFIN,RELUP,RELERR,RELER2,FCNMAX
COMMON /LMINUC/ NAMKIN(4),NAM(30)
CHARACTER*10 NAMKIN,NAM
C...LEPTOU: Cuts, basic switches and parameters.
DATA CUT/1.E-04,1.,0.,1.,4.,1.E+08,5.,1.E+08,1.,1.E+08,1.,1.E+08,
&0.,3.1416/
Ckc..differences with L65 (do not change !!)
Ckc..LST(1)=2 - (x,y) as independent variables
Ckc..LST(2)=3 - (x,y) supplied by user via LEPTOU, no cuts applied
Ckc..LST(6)=0 - no phi rotation,
Ckc..LST(17)=1 - variable energies of initial particles
Ckc..LST(18)=0 - no running alpha, is already done in HERACLES
C... 0 1 2 3 4 5 6 7 8 9
DATA LST/ 2, 3, 2, 1, 3, 0, 1, 12, 5,
1 1, 0, 6, 5, 4, 9, 1, 1, 0, -10,
2 5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3 0, 0, 0, 0, 2, 1, 5*0/
C... 0 1 2 3 4 5 6 7 8 9
DATA PARL/ 1., 1., 0.44,0.75,.2319,0., 0.5, .04, 4.,
1 0.44, 0.01,0.01,0.1, 0.35,0.01,7.29735E-03,
& 1.16639E-05,0.044,0.03,
2 0.1,10*0./
C...Internally used variables.
DATA PARI/40*0./
DATA QC/-.33333,.66667,-.33333,.66667,-.33333,.66667,
& -.33333,.66667/
DATA CABIBO/.95,.05,2*0.,.05,.948,.002,2*0.,.002,.998,4*0.,1./
DATA OPTX/1.,3*0./,OPTY/1.,3*0./,OPTQ2/1.,3*0./,OPTW2/1.,3*0./
DATA NXX,NWW/31,21/
DATA PQG,PQQB,QGMAX,QQBMAX/6510*0./,YCUT/651*0./,XTOT/651*0./
DATA NFX,NFQ/41,16/,FLQT,FLGT,FLMT/1968*0./
DATA XKIN/1.,2.,3.,4./,UKIN,WKIN,AIN,BIN/16*0./,MAXFIN/2000/
DATA RELUP,RELERR,RELER2/0.1,0.05,0.05/
DATA NAMKIN/' x',' ',' ',' '/
DATA IPY/
1 0, 0, 2, 2, 6, 1, 1, 6, 3, 1,
2 3, 1, 1, 2, 1, 1, 4, 1, 1, 1,
3 0, 1, 1, 1, 1, 1, 1, 0, 0, 0,
4 1, 2, 1, 1, 30, 33, 1, 1, 7, 0,
5 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
6 0, 0, 0, 1, 100, 0, 0, 0, 0, 0,
7 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
8 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/
DATA (PYPAR(I),I=1,40)/
1 7.299E-03, 2.290E-01, 2.000E-01, 2.500E-01, 4.000E+00,
1 1.000E+00, 4.400E-01, 4.400E-01, 7.500E-02, 0.000E+00,
2 2.000E+00, 2.000E+00, 1.000E+00, 0.000E+00, 3.000E+00,
2 1.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.000E+00,
3 2.500E-01, 1.000E+00, 2.000E+00, 1.000E-03, 1.000E+00,
3 1.000E+00, 1.000E+00, -2.000E-02, -1.000E-02, 0.000E+00,
4 0.000E+00, 1.600E+00, 0.500E+00, 0.200E+00, 3.894E-01,
4 1.000E+00, 3.300E-01, 6.600E-01, 0.000E+00, 1.000E+00/
DATA (PYPAR(I),I=41,80)/
5 2.260E+00, 1.000E+04, 1.000E-04, 0.000E+00, 0.000E+00,
5 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,
6 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,
6 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,
7 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,
7 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,
8 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,
8 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
DATA PYVAR/80*0./
END
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
SUBROUTINE DJGINIT(LEPIN,PLZ,PPZ,INTER)
chs..from LINIT, L61
chs..updated for L62, L63, L65
C...Initialize for an incoming lepton (type LEPIN, momentum pz=PLZ)
C...and target nucleon (momentum pz=PPZ) to interact via INTER.
C...Find maximum of differential cross section, calculate QCD event
C...probabilities or read them from logical file LFILE.
C...Numerical integration to obtain total cross-section.
ckc..LEPIN: e-(11), e+(-11)
ckc..LEPTO commons & declarations
COMMON /LINTRL/ PSAVE(3,4,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX,
&Q2MIN,Q2MAX,W2MIN,W2MAX,ILEP,INU,IG,IZ
COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U
COMMON /LINTER/ PARI(40),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17)
COMMON /LGRID/ NXX,NWW,XX(31),WW(21),PQG(31,21,3),PQQB(31,21,2),
&QGMAX(31,21,3),QQBMAX(31,21,2),YCUT(31,21),XTOT(31,21),NP
COMMON /LOPTIM/ OPTX(4),OPTY(4),OPTQ2(4),OPTW2(4),COMFAC
COMMON /LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
COMMON /LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
COMMON /LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
COMMON /LBOOST/ DBETA(2,3),DTHETA(2),DPHI(2),PB(5),PHIR
COMMON /LPFLAG/ LST3
COMMON /LYPARA/ IPY(80),PYPAR(80),PYVAR(80)
c...hs upgrade to Pythia6
double precision parp,paripy
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARIPY(200)
DOUBLE PRECISION DTHETA,DPHI,DBETA
DATA NCALL/0/
ckc..DJANGO
COMMON /DJPASS/ NTOT,NPASS,NFAILL
ckc..HERACLES
DOUBLE PRECISION
& HXMIN,HXMAX,HQ2MIN,HQ2MAX,HYMIN,HYMAX,HWMIN,HGMIN
& ,POLARI
& ,SW,CW,SW2,CW2
& ,MW,MZ,MH,ME,MMY,MTAU,MU,MD,MS,MC,MB,MT
& ,MW2,MZ2,MH2,ME2,MMY2,MTAU2,MU2,MD2,MS2,MC2,MB2,MT2
& ,PI,ALPHA,ALP1PI,ALP2PI,ALP4PI,E,GF,SXNORM,SX1NRM
& ,DELTAR,AGF0,DRHOT,DALPMZ,XGMT,ALPQCD,BTOP4,DRPIW2
& ,SIGTOT,SIGTRR,SIGG,SIGGRR
& ,PLZ,PPZ
COMMON /HSCUTS/ HXMIN,HXMAX,HQ2MIN,HQ2MAX,HYMIN,HYMAX,HWMIN,HGMIN
COMMON /HSPARL/ LPAR(20),LPARIN(12)
COMMON /HSPARM/ POLARI,HPOLAR,LLEPT,LQUA
COMMON /HSGSW/ SW,CW,SW2,CW2
* ,MW,MZ,MH,ME,MMY,MTAU,MU,MD,MS,MC,MB,MT
* ,MW2,MZ2,MH2,ME2,MMY2,MTAU2,MU2,MD2,MS2,MC2,MB2,MT2
COMMON /HSKNST/ PI,ALPHA,ALP1PI,ALP2PI,ALP4PI,E,GF,SXNORM,SX1NRM
COMMON /HSDELR/ DELTAR,AGF0,DRHOT,DALPMZ,XGMT,ALPQCD,BTOP4,DRPIW2
COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20)
COMMON /HSPDFO/ IPDFOP,IFLOPT,LQCD,LTM,LHT
chs-20.02.2016: Do not initialize for nucleon
c COMMON /HSNUCL/ HNA,HNZ,INUMOD
c DOUBLE PRECISION HNA,HNZ
C...for call to ariadne
LOGICAL LARIAD
COMMON /HSARIA/ LARIAD
C...initialize Block data for LEPTO and PYTHIA
EXTERNAL LUDATA
EXTERNAL PYDATA
ckc..temporarily
LFILE=0
chs..variable energy
LST(17)=1
chs..x and y are independent variables, transferred from HERACLES
LST(1)=2
LST(2)=3
C...Read users choice for LEPTO
ch CALL DJGCHC <-- moved to heracles (14/12/93)
C...Initialize ariadne
IF (LST(8).EQ.9) THEN
CALL ARINIT('LEPTO')
CALL ARPRDA
LARIAD=.TRUE.
ENDIF
C...F_L - from HS44 input
LST(11)=IFLOPT
ch...initialization for parton distributions done in HERACLES
C...top mass
PMAS(6,1)=MT
C...Z mass
PMAS(23,1)=MZ
C...W mass
PMAS(24,1)=MW
C...sin^2(\theta_W)
PARL(5)=SW2
C...lepton polarization
PARL(6)=POLARI
IF (ABS(POLARI).GE.0.99) PARL(6)=0.99d0*SIGN(1d0,POLARI)
C...alpha ew
PARL(16)=ALPHA
C...G_F
PARL(17)=GF
C...\Delta r
PARL(18)=DELTAR
C...Cross section within user defined cuts (in picobarn)
PARL(23)=SIGTOT*1E3
PARL(24)=SIGTOT*1E3
C...Relative accuracy of cross section
ACCUR=SIGTRR/SIGTOT
ckc..LEPTO
NCALL=NCALL+1
LST3=LST(3)
C...Couplings between Z0 and left/right-handed leptons and quarks.
ZL(1,1)=-.5+PARL(5)
ZL(1,2)=PARL(5)
ZL(2,1)=ZL(1,2)
ZL(2,2)=ZL(1,1)
ZL(1,3)=0.5
ZL(2,3)=0.
ZL(1,4)=0.
ZL(2,4)=0.5
DO 10 IFL=1,8
ZQ(1,IFL)=SIGN(0.5,QC(IFL))-QC(IFL)*PARL(5)
10 ZQ(2,IFL)=-QC(IFL)*PARL(5)
C...Set initial state.
chs-20.02.2016: Always initialize for proton, even for nuclei
chs scattering. Nucleon is chosen later
KSAVE(2)=2212
ckc..target nucleon
c PARL(1)=HNA
c PARL(2)=HNZ
LST(22)=1
LST(23)=INTER
KSAVE(1)=LEPIN
cC incoming proton:
c IF (INT(HNA).EQ.1.AND.INT(HNZ).EQ.1) THEN
c KSAVE(2)=2212
cC incoming neutron:
c ELSEIF (INT(HNA).EQ.1.AND.INT(HNZ).EQ.0) THEN
c KSAVE(2)=2112
cC other nuclei
c ELSE
c KSAVE(2)=1000000000+HNZ*10000+HNA*10
c ENDIF
K(1,1)=21
K(1,2)=KSAVE(1)
K(1,3)=0
K(1,4)=0
K(1,5)=0
K(2,1)=21
K(2,2)=KSAVE(2)
K(2,3)=0
K(2,4)=0
K(2,5)=0
ckc
P(1,1)=0.
P(1,2)=0.
P(1,3)=PLZ
P(1,5)=ULMASS(KSAVE(1))
P(1,4)=SQRT(P(1,3)**2+P(1,5)**2)
P(2,1)=0.
P(2,2)=0.
P(2,3)=PPZ
P(2,5)=ULMASS(KSAVE(2))
P(2,4)=SQRT(P(2,3)**2+P(2,5)**2)
N=2
LST(28)=3
C...Save momentum vectors of incoming particles
DO 20 I=1,2
DO 20 J=1,5
20 PSAVE(3,I,J)=P(I,J)
C...Dot-product of initial particles, cms energy
PARL(21)=2.*(DBLE(P(1,4))*DBLE(P(2,4))-DBLE(P(1,3))*DBLE(P(2,3)))
ROOTS=SQRT((DBLE(P(1,4))+DBLE(P(2,4)))**2
& -(DBLE(P(1,3))+DBLE(P(2,3)))**2)
IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1)) WRITE(6,1000)
&LEPIN,(P(1,J),J=1,3),PARL(1),PARL(2),(P(2,J),J=1,3),INTER,ROOTS
IF(PLZ*PPZ.GT.0.1) THEN
WRITE(6,1010)
STOP
ENDIF
C...Reduced header for Jetset/Pythia
MSTU(12)=0
MSTP(122)=0
IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1))
&WRITE(6,1020) MSTU(181),MSTU(182),MSTP(181),MSTP(182)
C...If JETSET version before 7.402, problem with azimuthal dependence
C...in LUSHOW solved by chosing flat azimuthal dependence.
IF(MSTU(181).LE.7.AND.MSTU(182).LT.402) THEN
MSTJ(46)=0
WRITE(6,1030) MSTJ(46)
ENDIF
C...Initialize PYTHIA for parton densities.
CALL HSPVER(0.1D0,100D0)
C...Set switches and parameters for parton densities in PYSTFU.
chs IF(LST(15).GT.0) THEN
chs initialization of PYINIT already done in HSPVER
chs MSTP(51)=LST(15)
chs MSTP(52)=LST(16)
chs MSTP(58)=LST(12)
chs ENDIF
chs CALL PYINIT('NONE','e-','p',ROOTS)
chs PARL(26)=PARP(1)
CAE-- use Lambda from parton densities in initial cascade
PYPAR(21)=PARP(1)
C...Reset PYTHIA 4.8 parameters from LEPTO parameters.
IF(MOD(LST(8),10).EQ.3.OR.MOD(LST(8),10).EQ.5) IPY(13)=0
IF(LST(35).EQ.0.AND.
&(MOD(LST(8),10).EQ.4.OR.MOD(LST(8),10).EQ.5)) IPY(14)=0
IPY(8)=LST(12)
IF(PSAVE(3,1,3).LT.0.) THEN
C...Flip event to have initial lepton along +z axis
P(1,3)=-P(1,3)
P(2,3)=-P(2,3)
ENDIF
C...Boost parameters to cms of incoming particles
DBETA(1,1)=0.D0
DBETA(1,2)=0.D0
DBETA(1,3)=(DBLE(P(1,3))+DBLE(P(2,3)))/(DBLE(P(1,4))+DBLE(P(2,4)))
DPHI(1)=0.D0
DTHETA(1)=0.D0
IF(LST(17).NE.0) THEN
C...For varying beam energies, transform to cms, lepton along +z axis.
CALL DUDBRB(0,0,0.D0,0.D0,0.D0,0.D0,-DBETA(1,3))
DPHI(1)=ULANGL(P(1,1),P(1,2))
CALL DUDBRB(0,0,0.D0,-DPHI(1),0.D0,0.D0,0.D0)
DTHETA(1)=ULANGL(P(1,3),P(1,1))
CALL DUDBRB(0,0,-DTHETA(1),0.D0,0.D0,0.D0,0.D0)
LST(28)=2
ENDIF
C...Effective limits on kinematic variables x, y, Q**2, W**2
PM2=P(2,5)**2
S=PARL(21)
Ckc..cuts on leptonic x,y,q2, taken from HERACLES
CUT(1)=HXMIN
CUT(2)=HXMAX
CUT(3)=HYMIN
CUT(4)=HYMAX
CUT(5)=HQ2MIN
CUT(6)=S
CUT(7)=HWMIN**2
CUT(8)=S+PM2
CUT(9)=0.0
CUT(10)=S/(2.*P(2,5))
CUT(11)=0.0
CUT(12)=PSAVE(3,2,4)
CUT(13)=0.0
CUT(14)=PI
chs...for use in LEPTO routines
XMIN=MAX(CUT(1),0.)
XMAX=MIN(CUT(2),1.)
YMIN=MAX(CUT(3),0.)
YMAX=MIN(CUT(4),1.)
Q2MIN=MAX(CUT(5),0.)
Q2MAX=MIN(CUT(6),S)
W2MIN=MAX(CUT(7),0.)
W2MAX=MIN(CUT(8),S+PM2)
UMIN=MAX(CUT(9),0.)
UMAX=MIN(CUT(10),S/(2.*P(2,5)))
chs..for documentation only
Q2MIN=MAX(CUT(5),S*CUT(1)*CUT(3))
Q2MAX=MIN(CUT(6),S*CUT(2)*CUT(4),S*CUT(4)-CUT(7)+PM2)
W2MIN=MAX(W2MIN,S*(1.-XMAX)*YMIN+PM2,Q2MIN*(1.-XMAX)/XMAX+PM2,
&S*YMIN-Q2MAX+PM2,2.*P(2,5)*UMIN*(1.-XMAX)+PM2)
W2MAX=MIN(W2MAX,S*(1.-XMIN)*YMAX+PM2,
&Q2MAX*(1.-XMIN)/MAX(XMIN,1.E-22)+PM2,
&S*YMAX-Q2MIN+PM2,2.*P(2,5)*UMAX*(1.-XMIN)+PM2)
ckc..final cuts written out
IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1)) WRITE(6,1050)
&CUT,XMIN,XMAX,YMIN,YMAX,Q2MIN,Q2MAX,W2MIN,W2MAX,UMIN,UMAX
IF(XMAX.LT.XMIN.OR.YMAX.LT.YMIN.OR.Q2MAX.LT.Q2MIN.OR.
&W2MAX.LT.W2MIN) THEN
IF(LST(3).GE.1) WRITE(6,1100)
IF(LST(3).GE.2) THEN
WRITE(6,1900)
STOP
ENDIF
ENDIF
IF(XMIN.LT.1.E-10.OR.Q2MIN.LT.1.E-01) THEN
IF(LST(3).GE.1) WRITE(6,1110)
IF(LST(3).GE.2) THEN
WRITE(6,1900)
STOP
ENDIF
ENDIF
IF (SQRT(W2MIN).LT.5D0) THEN
WRITE(6,1051) HWMIN
ENDIF
PARI(11)=(PARL(1)-PARL(2))/PARL(1)
KSAVE(4)=LEPIN
ILEP=1
IF(LEPIN.LT.0) ILEP=2
INU=0
IF(IABS(LEPIN).EQ.12.OR.IABS(LEPIN).EQ.14
&.OR.IABS(LEPIN).EQ.16) INU=1
IF(INU.EQ.1) THEN
C...Set full polarisation for incoming neutrino.
PARL(6)=-1.
IF(LEPIN.LT.0) PARL(6)=1.
ENDIF
IF(LST(23).EQ.1.AND.INU.EQ.0) THEN
C...Electromagnetic interaction.
KSAVE(3)=22
IG=1
IZ=0
ELSEIF(LST(23).EQ.2) THEN
C...Weak charged current, only one helicity state contributes.
IF(KSAVE(1).LT.0.AND.PARL(6).LT.-0.99
& .OR.KSAVE(1).GT.0.AND.PARL(6).GT.0.99) THEN
IF(LST(3).GE.1) WRITE(6,1150) LEPIN,PARL(6)
IF(LST(3).GE.2) THEN
WRITE(6,1900)
STOP
ENDIF
ENDIF
IF(MOD(IABS(LEPIN),2).EQ.0) THEN
KSAVE(3)=ISIGN(24,LEPIN)
KSAVE(4)=ISIGN(IABS(LEPIN)-1,LEPIN)
ELSE
KSAVE(3)=ISIGN(24,-LEPIN)
KSAVE(4)=ISIGN(IABS(LEPIN)+1,LEPIN)
ENDIF
ELSEIF(LST(23).EQ.3.OR.(LST(23).EQ.4.AND.INU.EQ.1)) THEN
C...Weak neutral current.
KSAVE(3)=23
IG=0
IZ=1
ELSEIF(LST(23).EQ.4.AND.INU.EQ.0) THEN
C...Neutral current, electromagnetic and weak with interference.
KSAVE(3)=23
IG=1
IZ=1
ELSE
IF(LST(3).GE.1) WRITE(6,1200) INTER,LEPIN
IF(LST(3).GE.2) THEN
WRITE(6,1900)
STOP
ENDIF
ENDIF
C...Choice of independent variables.
IF(LST(1).EQ.0) THEN
LST(31)=1
IF(INTER.EQ.2.OR.INTER.EQ.3) LST(31)=2
ELSE
LST(31)=IABS(LST(1))
ENDIF
IF(LST(23).EQ.2) THEN
C...Constant factor GF**2/pi for CC, transformation to picobarn.
PARI(31)=PARL(17)**2/PI*0.39E+09
ELSE
C...Constant factor 2*pi*alpha**2 for NC, transformation to picobarn.
PARI(31)=2.*PI*PARL(16)**2*0.39E+09
ENDIF
ckc..settings of lepto flags printed out
IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1))
&WRITE(6,1250) (I,LST(I),LST(I+10),PARL(I),PARL(I+10),I=1,10)
ckc..from DJANGO
WRITE(6,889) PARL(23),ACCUR
ckc..this will be needed for ME option
IF((LST(19).GE.0.OR.LST(19).EQ.-10).AND.
&(LST(8).EQ.1.OR.LST(8)/10.EQ.1.OR.MOD(LST(8),10).EQ.9)) THEN
C...Calculate weights if 1st order QCD requested.
CALL LTIMEX(TI1)
CALL LWEITS(LFILE)
CALL LTIMEX(TI2)
IF(LST(3).GE.4.OR.(LST(3).EQ.3.AND.NCALL.EQ.1))
& WRITE(6,1510) TI2-TI1
ENDIF
C...Reset counters
NFAILL=0
NPASS=0
NTOT=0
C...Reset counters to zero for Monte Carlo estimate of cross section.
LST(32)=0
RETURN
1000 FORMAT(' ',//,5X,
&'A MONTE CARLO GENERATOR FOR DEEP INELASTIC LEPTON-'
&,'NUCLEON SCATTERING',/,5X,68('='),//,
&25X,'LEPTO version 6.5.1, October 31, 1996',//,
&' Lepton: type =',I3,5X,'momentum (px,py,pz) =',3F8.1,
&' GeV',//,' Target: A, Z =',2F3.0,2X,
&'momentum (px,py,pz) =',3F8.1,' GeV',//,
&' Interaction :',I3,14X,' CMS energy =',1PG12.4,' GeV',/)
1010 FORMAT(' Warning: lepton and nucleon momenta in same direction',
&' not allowed.',/,10X,'Execution stopped.')
1020 FORMAT(/,' JETSET version ',I3,'.',I3,' is used.',/,
&' Parton densities in PYTHIA version ',I3,'.',I3,' are used.',/)
1030 FORMAT(' Warning (LINIT): JETSET version before 7.402, MSTJ(46)'
& ,' set to',I4,/,18X,'to avoid mismatch LEPTO<-->LUSHOW.',/)
1050 FORMAT(/,' User applied cuts (+ phase space) : ',1P,
& G12.4,' < x < ',G12.4,
&/,37X,G12.4,' < y < ',G12.4,
&/,37X,G12.4,' < Q**2 < ',G12.4,
&/,37X,G12.4,' < W**2 < ',G12.4,
&/,37X,G12.4,' < nu < ',G12.4,
&/,37X,G12.4,' < E'' < ',G12.4,
&/,37X,G12.4,' < theta < ',G12.4,/,
&/, ' Effective ranges (from above cuts): ',
& G12.4,' < x < ',G12.4,
&/,37X,G12.4,' < y < ',G12.4,
&/,37X,G12.4,' < Q**2 < ',G12.4,
&/,37X,G12.4,' < W**2 < ',G12.4,
&/,37X,G12.4,' < nu < ',G12.4)
1051 FORMAT(' ********************** Warning ***********************'
& /,' Minimum of hadronic final state mass (',F5.1,' GeV)',/
& ,' is less than 5 GeV. ',/
& ,' Phase space might be reduced and the total cross section',/
& ,' possibly be renormalized.',/
& ,' The correct value for the total cross section is stored in'
& ,' ---> PARL(24) <---',/
& ,' ********************************************************',/)
1100 FORMAT(' Warning: effective upper limit of kinematical ',
&'variable(s) smaller than corresponding lower limit.')
1110 FORMAT(' Warning: lower limit in x and/or Q2 too small for ',
&'DIS formalism.')
1150 FORMAT(' Warning: weak charged current cross section zero for ',
&'specified lepton helicity; LEPIN, PARL(6) =',I3,F5.2)
1200 FORMAT(' Warning: unrecognized interaction in LINIT call: ',
&'INTER = ',I5,' for lepton LEPIN =',I5)
c 1210 FORMAT(' Warning: unallowed value of LST(1) =',I3,
c &' and/or LST(31) =',I3)
c 1220 FORMAT(/,' User-defined optimization parameters:',
c &/,5X,'OPTX(1...4) =',4G11.3,/,5X,'OPTY(1...4) =',4G11.3,
c &/,5X,'OPYQ2(1...4) =',4G11.3,/,5X,'OPTW2(1...4) =',4G11.3,/)
1250 FORMAT(/,' Parameter values:',//,9X,'I',4X,'LST(I)',1X,
&'LST(I+10)',8X,'PARL(I)',5X,'PARL(I+10)',1P,
&/,5X,55('-'),10(/,3I10,2G15.4),/)
889 FORMAT(' Cross-section from HERACLES =',
&G12.5,' pb',
&/,' Estimated relative error = ',G10.2,/)
1510 FORMAT(/,' Time for calculating QCD weights =',F5.1,' seconds',/)
1900 FORMAT(' Execution stopped ',/)
END
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
SUBROUTINE DJGCHC(IHSONL)
CHARACTER*80 TITLE
CHARACTER*10 CODE,CODEWD
DIMENSION CODE(20)
COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),X,Y,W2,Q2,U
COMMON /HSUNTS/ LUNTES,LUNDAT,LUNIN,LUNOUT,LUNRND
ck..Ariadne
COMMON /ARDAT1/ PARA(40),MSTA(40)
cs..Sophia
double precision WSOPHIA
COMMON /SOPHCT/ WSOPHIA
C
C-----------------------------------------------------------------------
DATA CODE/
1'TITLE ','OUT-LEP ','FRAME ','FRAG ','CASCADES ',
2'BARYON ','MAX-VIRT ','KT-PARTON ','DIQUARK ','KT-REMNANT',
3'AR-REMNANT',' ',' ',' ',' ',
4'SOPHIA ',' ',' ',' ','CONTINUE '/
DATA TITLE/' '/
C
IHSONL=0
C---PRINT 2nd TITLE
WRITE(LUNOUT,9)
9 FORMAT(1X,///
1'**************************************************',
2'*****************************',
3/,10X,' '
4/,10X,' DJANGOH '
5/,10X,' version 1.9 '
6/,10X,' '
7/,10X,' HERACLES 4.6.13 + LEPTO 6.5.1 '
8/,10X,' '
9/,10X,' author: H.Spiesberger '
9/,10X,' with major contributions from '
9/,10X,' G.A.Schuler and K.Charchula '
1/,10X,' '
2/,10X,' last change: 27 May 2016 (by HS) '
3/,10X,' '
4/,' **************************************************',
5'****************************',//)
C
C***********************************************************************
C READ INPUT DATA
C
C STRUCTURE OF INPUT:
C 1) CODEWD (A10)
C 2) CORRESPONDING DATA (FORMAT FREE)
C***********************************************************************
C
1 CONTINUE
READ(LUNIN,90,END=4) CODEWD
WRITE(LUNOUT,91) CODEWD
DO 2 ISW=1,20
IF(CODEWD.EQ.CODE(ISW))GO TO 3
2 CONTINUE
WRITE(LUNOUT,92)
GO TO 1
3 GO TO(
C------------------------------------------------------------------
C TITLE , OUT-LEP , FRAME , FRAG , CASCADES ,
1 100 , 200 , 300 , 400 , 500 ,
C
C------------------------------------------------------------------
C BARYON , MAX-VIRT , KT-PARTON , DIQUARK , KT-REMNANT,
2 600 , 700 , 800 , 900 , 1000 ,
C
C------------------------------------------------------------------
C AR-REMNANT , , , , ,
3 1100 , 1200 , 1300 , 1400 , 1500 ,
C
C------------------------------------------------------------------
C SOPHIA , , , , CONTINUE ,
4 1600 , 1700 , 1800 , 1900 , 2000 )
C
C------------------------------------------------------------------
9,ISW
GO TO 1
4 CONTINUE
IHSONL=1
WRITE(6,93)
GOTO 2000
C
90 FORMAT(A10)
91 FORMAT(//' *****NEXT CONTROL CARD ***** ',A10/)
92 FORMAT(/,' UNKNOWN CODEWORD - CONTROL CARD IGNORED')
93 FORMAT(/,' NO INPUT FOR DJANGO6 - RUN HERACLES ONLY.')
C
C***********************************************************************
C CONTROL CARD: CODEWD = TITLE
C DEFINES THE TITLE OF THE JOB
C***********************************************************************
100 CONTINUE
READ(LUNIN,190) TITLE
WRITE(6,191) TITLE
GO TO 1
190 FORMAT(A80)
191 FORMAT(//,5X,A80,//)
C
C***********************************************************************
C CONTROL CARD: CODEWD = OUT-LEP
C Regulates information in event record (D=1)
C LST(4) = I_lepton + 10*I_shower
C I_lepton = 0/1 inactive/active scattered electron
C I_shower = 0/1 excluded/included interm. partons
C***********************************************************************
200 CONTINUE
READ(LUNIN,*) LST(4)
WRITE(LUNOUT,'(5X,A,I3)')
* ' LST(4)=',LST(4)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = FRAME
C Choice of frame for event (D=3)
C LST(5) =1 : hadronic CM frame, z-axis along boson
C =2 : lepton-nucleon CM frame, z-axis along lepton
C =3 : lab system as defined by HS (e in + z direction)
C =4 : as 3 but z-axis along exchanged boson
C
C***********************************************************************
300 CONTINUE
READ(LUNIN,*) LST(5)
WRITE(LUNOUT,'(5X,A,I3)')
* ' LST(5)=',LST(5)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = FRAG
C Fragmentation (D=1)
C LST(7) = 0 : only parton level
C = 1 : as 0 + hadronization and decays
C***********************************************************************
400 CONTINUE
READ(LUNIN,*) LST(7)
WRITE(LUNOUT,'(5X,A,I3)')
* ' LST(7)=',LST(7)
IF (LST(7).LT.0) THEN
IHSONL=1
WRITE(LUNOUT,'(5X,A,I3)')
* ' RUN ONLY HERACLES, ADDITIONAL DJANGO6-INPUT IGNORED'
GOTO 2000
ENDIF
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = CASCADES
C Simulation of QCD events (D=12)
C LST(8)=0 : QCD events switched off
C =1 : 1st order QCD processes (gluon radiation and boson-
C gluon fusion) according to exact matrix elements.
C =2 : QCD parton cascade evolution of initial and final quark
C =3 : QCD parton cascade evolution of initial quark only.
C =4 : QCD parton cascade evolution of final quark only.
C =5 : QCD switched off, but target treatment exactly as in
C parton cascade case.
C =9 : ARIADNE: simulating parton emmision in the
C colour dipole model
C =12-15: ME+PS (as 2-5)
C
C***********************************************************************
500 CONTINUE
READ(LUNIN,*) LST(8)
WRITE(LUNOUT,'(5X,A,I3)')
* ' LST(8)=',LST(8)
IF(LST(8).NE.0.AND.LST(8).NE.1.AND.
& LST(8).NE.2.AND.LST(8).NE.3.AND.LST(8).NE.4.AND.LST(8).NE.5
&.AND.LST(8).NE.9.AND.
& LST(8).NE.12.AND.LST(8).NE.13.AND.LST(8).NE.14.AND.LST(8).NE.15)
& THEN
LST(8)=12
WRITE(LUNOUT,'(5X,A,I3)')
* 'Warning: LST(8) out of range, set to: ',LST(8)
ENDIF
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = BARYON
C Treatment of targer remant (D=1)
C LST(14) = 0 : as an anti-parton
C = 1 : into baryon
C = 2 & 3 : as 1 but with different probability distributions
C***********************************************************************
600 CONTINUE
READ(LUNIN,*) LST(14)
WRITE(LUNOUT,'(5X,A,I3)')
* ' LST(14)=',LST(14)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = MAX-VIRT
C max virtuality in parton cascades (D=5)
C LST(9) =1: Q^2
C =2: W^2
C =3: W*Q
C =4: Q^2*(1-x)
C =5: Q^2*(1-x)*max(1,ln1/x)
C =9: W^4/3, i.e. similar as in dipole model
C***********************************************************************
700 CONTINUE
READ(LUNIN,*) LST(9)
WRITE(LUNOUT,'(5X,A,I3)')
* ' LST(9)=',LST(9)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = KT-PARTON
C width of gaussian primordial transverse momentum distr.
C
C PARL(3) = 0.44 GeV (Default)
C***********************************************************************
800 CONTINUE
READ(LUNIN,*) PARL(3)
WRITE(LUNOUT,'(5X,A,1PE13.4)')
* ' PARL(3)=',PARL(3)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = DIQUARK
C Probability that a ud-diquark in target remnant has
C spin and isospin equal zero, i.e. I=S=0
C PARL(4) = 0.75 (Default)
C***********************************************************************
900 CONTINUE
READ(LUNIN,*) PARL(4)
WRITE(LUNOUT,'(5X,A,1PE13.4)')
* ' PARL(4)=',PARL(4)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = KT-REMNANT
C width of gaussian pt when non-trivial target remnant is
C split into a particle and a jet
C PARL(14) = 0.44 GeV (Default)
C***********************************************************************
1000 CONTINUE
READ(LUNIN,*) PARL(14)
WRITE(LUNOUT,'(5X,A,1PE13.4)')
* ' PARL(14)=',PARL(14)
GO TO 1
C
C***********************************************************************
C CONTROL CARD: CODEWD = AR-REMNANT
C Ariadne: regulates remnant vs struck quark dipole (D=1)
C MSTA(30) = 0 : struck quark pointlike, mu=PARA(11)
C = 1 : as 1 mu=PARA(11)/(1-x)
C = 2 : as 1 , mu=Q
C***********************************************************************
1100 CONTINUE
READ(LUNIN,*) MSTA(30)
WRITE(LUNOUT,'(5X,A,I3)')
* ' MSTA(30)=',MSTA(30)
GO TO 1
1200 CONTINUE
1300 CONTINUE
1400 CONTINUE
1500 CONTINUE
C
C***********************************************************************
C CONTROL CARD: CODEWD = SOPHIA
C for Sophia: maximal value of Whad for which
C the hadronic final state is generated by Sophia
C WSOPHIA = 1.5 GeV (Default)
C***********************************************************************
1600 CONTINUE
READ(LUNIN,*) WSOPHIA
WRITE(LUNOUT,'(5X,A,1PE13.4)')
* ' WSOPHIA=',WSOPHIA
GO TO 1
1700 CONTINUE
1800 CONTINUE
1900 CONTINUE
GO TO 1
C
C***********************************************************************
2000 CONTINUE
RETURN
C
END
c
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
SUBROUTINE HSWCUT(W,IQFC,IQFRC,IFAIL)
ckc..from DJANGO; modified
ckc..IQF in old quark convention; should be modified
IMPLICIT double precision (A-H,M,O-Z)
ckc..from HS
COMMON /HSNUME/ SIGTOT,SIGTRR,SIGG(20),SIGGRR(20),NEVENT,NEVE(20)
LOGICAL LFIRST
DATA LFIRST/.TRUE./
ckc..from L61
COMMON /LEPTOU/CUT(14),LST(40),PARL(30),XHAD,YHAD,W2HAD,Q2HAD,UHAD
REAL CUT,PARL,ULMASS,XHAD,YHAD,W2HAD,Q2HAD,UHAD,WR
C...L61 (L52 also) uses old PYTHIA common/PYPARA/
COMMON /LYPARA/ IPY(80),PYPAR(80),PYVAR(80)
REAL PYPAR ,PYVAR
C
ckc..
WR=W
IQF=IQFC
IQFR=IQFRC
ckc..
IF (LFIRST) THEN
LFIRST=.FALSE.
NEVTRN=0
NREJCT=0
ENDIF
NEVTRN=NEVTRN+1
IFAIL=0
ckc..change of u,d code
IF(ABS(IQF).EQ.1) THEN
IQF=IQF+ISIGN(1,IQF)
ELSEIF(ABS(IQF).EQ.2) THEN
IQF=IQF-ISIGN(1,IQF)
ENDIF
cd IF(LST(14).EQ.0.OR.IQFR.GT.10.OR.LST(8).GE.2) THEN
IF(LST(14).EQ.0.OR.IQFR.GT.10
&.OR.(LST(8).GE.2.AND.MOD(LST(8),10).NE.9)) THEN
C...Check if energy in jet system is enough for fragmentation.
IF(WR.LT.ULMASS(2212)+ULMASS(IQF)+PYPAR(12)
& +SQRT(PARL(3)**2+ULMASS(IQF)**2)) THEN
IFAIL=1
ckc..here flavour of the lightest quark should be given
IF (IQF.EQ.2) THEN
NREJCT=NREJCT+1
chs PARL(24)=SIGTOT*(1D0-DFLOAT(NREJCT)/NEVTRN)*1E3
chs..moved correction total cross section do HERACLES
ENDIF
ENDIF
ENDIF
RETURN
END
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
SUBROUTINE DJGVAR(ICHNN,X,Y,Q2)
C---
C Transfer event characteristics from HERACLES to LEPTO
C---
IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
COMMON /HSELAB/ SP,EELE,PELE,EPRO,PPRO
COMMON /HSIKP/ S,T,U,SS,TS,US,DKP,DKPS,DKQ,DKQS
COMMON /HSGSW1/ MEI,MEF,MQI,MQF,MEI2,MEF2,MQI2,MQF2,MPRO,MPRO2
PARAMETER (NMXHEP=2000)
COMMON /HEPEVT/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
& JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),
& PHEP(5,NMXHEP),VHKK(4,NMXHEP)
COMMON /DJRSCL/ HSSCH,HQ2SCH,HW2SCH,HXSCH,HYSCH
& ,SLIX,SLIY,SLIZ,SLIE,SLIM
& ,SLFX,SLFY,SLFZ,SLFE,SLFM
COMMON /LEPTOU/CUT(14),LST(40),PARL(30),XSCH,YSCH,W2SCH,Q2SCH,USCH
C---declarations for LEPTO
REAL CUT,PARL,XSCH,YSCH,W2SCH,Q2SCH,USCH
C---event characteristics in common HEPEVT
C electron: IDHEP(1), PHEP(i,1)
C quark: IDHEP(2), PHEP(i,2)
C photon: IDHEP(3), PHEP(i,3)
C---notation for rescaled 4-momenta
C SLIX,SLIY,SLIZ: rescaled initial lepton momentum
C SLIE: energy of initial lepton
C SLIM: off-shell mass of the initial lepton
C---
ck..25/01 neutrino mass set to zero
DATA MNU/0./
IF (ICHNN.LE.2) THEN
C---non-radiative channels