← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lplm.f

  • Committer: Marco Zaro
  • Date: 2014-01-27 16:54:10 UTC
  • mfrom: (78.124.55 MG5_aMC_2.1)
  • Revision ID: marco.zaro@gmail.com-20140127165410-5lma8c2hzbzm426j
merged with lp:~maddevelopers/madgraph5/MG5_aMC_2.1 r 267

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/PYAnalyzer/mcatnlo_pyan_pp_lplm.f
 
1
c
 
2
c Example analysis for "p p > l+ l- [QCD]" process.
 
3
c
1
4
C----------------------------------------------------------------------
2
5
      SUBROUTINE RCLOS()
3
6
C     DUMMY IF HBOOK IS USED
9
12
      SUBROUTINE PYABEG
10
13
C     USER'S ROUTINE FOR INITIALIZATION
11
14
C----------------------------------------------------------------------
12
 
      real * 8 xm0,gam,xmlow,xmupp,bin
13
 
      real * 8 xmi,xms,pi
 
15
      implicit none
 
16
      include 'reweight0.inc'
 
17
      real * 8 bin,xmi,xms,pi
14
18
      PARAMETER (PI=3.14159265358979312D0)
15
 
      integer j,k,jpr
 
19
      integer j,kk,l
16
20
      character*5 cc(2)
17
21
      data cc/'     ',' cuts'/
18
 
c
 
22
      integer nwgt,max_weight,nwgt_analysis
 
23
      common/cnwgt/nwgt
 
24
      common/c_analysis/nwgt_analysis
 
25
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
26
      character*15 weights_info(max_weight)
 
27
      common/cwgtsinfo/weights_info
19
28
      call inihist
20
 
      do j=1,2
21
 
      k=(j-1)*50
 
29
      nwgt_analysis=nwgt
22
30
c
23
31
      xmi=50.d0
24
32
      xms=130.d0
25
33
      bin=0.8d0
26
 
      call mbook(k+ 1,'V pt'//cc(j),2.d0,0.d0,200.d0)
27
 
      call mbook(k+ 2,'V pt'//cc(j),10.d0,0.d0,1000.d0)
28
 
      call mbook(k+ 3,'V log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
29
 
      call mbook(k+ 4,'V y'//cc(j),0.2d0,-9.d0,9.d0)
30
 
      call mbook(k+ 5,'V eta'//cc(j),0.2d0,-9.d0,9.d0)
31
 
      call mbook(k+ 6,'mV'//cc(j),(bin),xmi,xms)
32
 
c
33
 
      call mbook(k+ 7,'l pt'//cc(j),2.d0,0.d0,200.d0)
34
 
      call mbook(k+ 8,'l pt'//cc(j),10.d0,0.d0,1000.d0)
35
 
      call mbook(k+ 9,'l log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
36
 
      call mbook(k+10,'l eta'//cc(j),0.2d0,-9.d0,9.d0)
37
 
      call mbook(k+11,'lb pt'//cc(j),2.d0,0.d0,200.d0)
38
 
      call mbook(k+12,'lb pt'//cc(j),10.d0,0.d0,1000.d0)
39
 
      call mbook(k+13,'lb log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
40
 
      call mbook(k+14,'lb eta'//cc(j),0.2d0,-9.d0,9.d0)
41
 
c
42
 
      call mbook(k+15,'llb delta eta'//cc(j),0.2d0,-9.d0,9.d0)
43
 
      call mbook(k+16,'llb azimt'//cc(j),pi/20.d0,0.d0,pi)
44
 
      call mbook(k+17,'llb log[pi-azimt]'//cc(j),0.05d0,-4.d0,0.1d0)
45
 
      call mbook(k+18,'llb inv m'//cc(j),(bin),xmi,xms)
46
 
      call mbook(k+19,'llb pt'//cc(j),2.d0,0.d0,200.d0)
47
 
      call mbook(k+20,'llb log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
48
 
c
49
 
      call mbook(k+21,'total'//cc(j),1.d0,-1.d0,1.d0)
 
34
      do kk=1,nwgt_analysis
 
35
      do j=1,2
 
36
      l=(kk-1)*42+(j-1)*21
 
37
      call mbook(l+ 1,'V pt      '//weights_info(kk)//cc(j)
 
38
     &     ,2.d0,0.d0,200.d0)
 
39
      call mbook(l+ 2,'V pt      '//weights_info(kk)//cc(j)
 
40
     &     ,10.d0,0.d0,1000.d0)
 
41
      call mbook(l+ 3,'V log[pt] '//weights_info(kk)//cc(j)
 
42
     &     ,0.05d0,0.1d0,5.d0)
 
43
      call mbook(l+ 4,'V y       '//weights_info(kk)//cc(j)
 
44
     &     ,0.2d0,-9.d0,9.d0)
 
45
      call mbook(l+ 5,'V eta     '//weights_info(kk)//cc(j)
 
46
     &     ,0.2d0,-9.d0,9.d0)
 
47
      call mbook(l+ 6,'mV        '//weights_info(kk)//cc(j)
 
48
     &     ,bin,xmi,xms)
 
49
c
 
50
      call mbook(l+ 7,'lm pt      '//weights_info(kk)//cc(j)
 
51
     &     ,2.d0,0.d0,200.d0)
 
52
      call mbook(l+ 8,'lm pt      '//weights_info(kk)//cc(j)
 
53
     &     ,10.d0,0.d0,1000.d0)
 
54
      call mbook(l+ 9,'lm log[pt] '//weights_info(kk)//cc(j)
 
55
     &     ,0.05d0,0.1d0,5.d0)
 
56
      call mbook(l+10,'lm eta     '//weights_info(kk)//cc(j)
 
57
     &     ,0.2d0,-9.d0,9.d0)
 
58
      call mbook(l+11,'lp pt      '//weights_info(kk)//cc(j)
 
59
     &     ,2.d0,0.d0,200.d0)
 
60
      call mbook(l+12,'lp pt      '//weights_info(kk)//cc(j)
 
61
     &     ,10.d0,0.d0,1000.d0)
 
62
      call mbook(l+13,'lp log[pt] '//weights_info(kk)//cc(j)
 
63
     &     ,0.05d0,0.1d0,5.d0)
 
64
      call mbook(l+14,'lp eta     '//weights_info(kk)//cc(j)
 
65
     &     ,0.2d0,-9.d0,9.d0)
 
66
c
 
67
      call mbook(l+15,'lmlp delta eta     '//weights_info(kk)//cc(j)
 
68
     $     ,0.2d0,-9.d0,9.d0)
 
69
      call mbook(l+16,'lmlp azimt         '//weights_info(kk)//cc(j)
 
70
     $     ,pi/20.d0,0.d0,pi)
 
71
      call mbook(l+17,'lmlp log[pi-azimt] '//weights_info(kk)//cc(j)
 
72
     $     ,0.05d0,-4.d0,0.1d0)
 
73
      call mbook(l+18,'lmlp inv m         '//weights_info(kk)//cc(j)
 
74
     $     ,bin,xmi,xms)
 
75
      call mbook(l+19,'lmlp pt            '//weights_info(kk)//cc(j)
 
76
     $     ,2.d0,0.d0,200.d0)
 
77
      call mbook(l+20,'lmlp log[pt]       '//weights_info(kk)//cc(j)
 
78
     $     ,0.05d0,0.1d0,5.d0)
 
79
c
 
80
      call mbook(l+21,'total'//weights_info(kk)//cc(j),1.d0,-1.d0,1.d0)
 
81
      enddo
50
82
      enddo
51
83
 999  END
52
84
 
56
88
C     USER'S ROUTINE FOR TERMINAL CALCULATIONS, HISTOGRAM OUTPUT, ETC
57
89
C----------------------------------------------------------------------
58
90
      REAL*8 XNORM
59
 
      INTEGER I,J,K
 
91
      INTEGER I,J,KK,l,nwgt_analysis
 
92
      integer NPL
 
93
      parameter(NPL=15000)
 
94
      common/c_analysis/nwgt_analysis
60
95
      OPEN(UNIT=99,FILE='PYTLL.TOP',STATUS='UNKNOWN')
61
96
      XNORM=1.D0/IEVT
62
 
      DO I=1,100              
 
97
      DO I=1,NPL
63
98
        CALL MFINAL3(I)             
64
 
        CALL MCOPY(I,I+100)
65
 
        CALL MOPERA(I+100,'F',I+100,I+100,(XNORM),0.D0)
66
 
        CALL MFINAL3(I+100)             
 
99
        CALL MCOPY(I,I+NPL)
 
100
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
101
        CALL MFINAL3(I+NPL)             
67
102
      ENDDO                          
68
103
C
69
 
      do j=1,2
70
 
      k=(j-1)*50
71
 
      call multitop(100+k+ 1,99,3,2,'V pt',' ','LOG')
72
 
      call multitop(100+k+ 2,99,3,2,'V pt',' ','LOG')
73
 
      call multitop(100+k+ 3,99,3,2,'V log[pt]',' ','LOG')
74
 
      call multitop(100+k+ 4,99,3,2,'V y',' ','LOG')
75
 
      call multitop(100+k+ 5,99,3,2,'V eta',' ','LOG')
76
 
      call multitop(100+k+ 6,99,3,2,'mV',' ','LOG')
 
104
      do kk=1,nwgt_analysis
 
105
      do i=1,2
 
106
      l=(kk-1)*42+(i-1)*21
 
107
C
 
108
      call multitop(NPL+l+ 1,NPL-1,3,2,'V pt',' ','LOG')
 
109
      call multitop(NPL+l+ 2,NPL-1,3,2,'V pt',' ','LOG')
 
110
      call multitop(NPL+l+ 3,NPL-1,3,2,'V log[pt]',' ','LOG')
 
111
      call multitop(NPL+l+ 4,NPL-1,3,2,'V y',' ','LOG')
 
112
      call multitop(NPL+l+ 5,NPL-1,3,2,'V eta',' ','LOG')
 
113
      call multitop(NPL+l+ 6,NPL-1,3,2,'mV',' ','LOG')
 
114
c
 
115
      call multitop(NPL+l+ 7,NPL-1,3,2,'lm pt',' ','LOG')
 
116
      call multitop(NPL+l+ 8,NPL-1,3,2,'lm pt',' ','LOG')
 
117
      call multitop(NPL+l+ 9,NPL-1,3,2,'lm log[pt]',' ','LOG')
 
118
      call multitop(NPL+l+10,NPL-1,3,2,'lm eta',' ','LOG')
 
119
      call multitop(NPL+l+11,NPL-1,3,2,'lm pt',' ','LOG')
 
120
      call multitop(NPL+l+12,NPL-1,3,2,'lm pt',' ','LOG')
 
121
      call multitop(NPL+l+13,NPL-1,3,2,'lm log[pt]',' ','LOG')
 
122
      call multitop(NPL+l+14,NPL-1,3,2,'lm eta',' ','LOG')
 
123
c
 
124
      call multitop(NPL+l+15,NPL-1,3,2,'lmlp deta',' ','LOG')
 
125
      call multitop(NPL+l+16,NPL-1,3,2,'lmlp azi',' ','LOG')
 
126
      call multitop(NPL+l+17,NPL-1,3,2,'lmlp azi',' ','LOG')
 
127
      call multitop(NPL+l+18,NPL-1,3,2,'lmlp inv m',' ','LOG')
 
128
      call multitop(NPL+l+19,NPL-1,3,2,'lmlp pt',' ','LOG')
 
129
      call multitop(NPL+l+20,NPL-1,3,2,'lmlp pt',' ','LOG')
 
130
c
 
131
      call multitop(NPL+l+21,NPL-1,3,2,'total',' ','LOG')
77
132
      enddo
78
 
c
79
 
      do j=1,2
80
 
      k=(j-1)*50
81
 
      call multitop(100+k+ 7,99,3,2,'l pt',' ','LOG')
82
 
      call multitop(100+k+ 8,99,3,2,'l pt',' ','LOG')
83
 
      call multitop(100+k+ 9,99,3,2,'l log[pt]',' ','LOG')
84
 
      call multitop(100+k+10,99,3,2,'l eta',' ','LOG')
85
 
      call multitop(100+k+11,99,3,2,'l pt',' ','LOG')
86
 
      call multitop(100+k+12,99,3,2,'l pt',' ','LOG')
87
 
      call multitop(100+k+13,99,3,2,'l log[pt]',' ','LOG')
88
 
      call multitop(100+k+14,99,3,2,'l eta',' ','LOG')
89
 
c
90
 
      call multitop(100+k+15,99,3,2,'llb deta',' ','LOG')
91
 
      call multitop(100+k+16,99,3,2,'llb azi',' ','LOG')
92
 
      call multitop(100+k+17,99,3,2,'llb azi',' ','LOG')
93
 
      call multitop(100+k+18,99,3,2,'llb inv m',' ','LOG')
94
 
      call multitop(100+k+19,99,3,2,'llb pt',' ','LOG')
95
 
      call multitop(100+k+20,99,3,2,'llb pt',' ','LOG')
96
 
c
97
 
      call multitop(100+k+21,99,3,2,'total',' ','LOG')
98
133
      enddo
99
134
c
100
135
      CLOSE(99)
106
141
C----------------------------------------------------------------------
107
142
      implicit double precision(a-h, o-z)
108
143
      implicit integer(i-n)
 
144
      include 'reweight0.inc'
109
145
      DOUBLE PRECISION HWVDOT,PSUM(4),PPV(5),YCUT,XMV,PTV,YV,THV,ETAV,
110
146
     #  PPL(5),PPLB(5),PTL,YL,THL,ETAL,PLL,ENL,PTLB,YLB,THLB,ETALB,
111
147
     #  PLLB,ENLB,PTPAIR,DLL,CLL,AZI,AZINORM,XMLL,DETALLB,PPV0(5),
112
 
     #  PPL0(5),PPLB0(5)
 
148
     #  PPL0(5),PPLB0(5),P1(4),P2(4),PIHEP(4)
113
149
      INTEGER ICHSUM,ICHINI,IHEP,IV,IFV,IST,ID,IJ,ID1,JPR,IDENT,
114
150
     #  ILL,ILLB,IHRD
115
151
      integer pychge
120
156
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
121
157
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
122
158
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
123
 
      LOGICAL DIDSOF,TEST1,TEST2,TEST3,TEST4,TEST5,TEST6,TEST7,flag
124
 
      REAL*8 PI,wmass,wgamma,bwcutoff,getinvm,getdelphi,getpseudorap
 
159
      LOGICAL DIDSOF,flag,ISLP,ISLM,FOUNDP,FOUNDM
 
160
      REAL*8 PI,wmass,wgamma,bwcutoff,getinvm,getdelphi,getrapidity,
 
161
     &getpseudorap
125
162
      PARAMETER (PI=3.14159265358979312D0)
126
 
      REAL*8 WWW0,TINY,p1(4),p2(4),pihep(4)
127
 
      INTEGER KK
 
163
      REAL*8 WWW0,TINY
 
164
      INTEGER KK,i,l
128
165
      DATA TINY/.1D-5/
129
 
      DOUBLE PRECISION EVWEIGHT,GETRAPIDITY
 
166
      DOUBLE PRECISION EVWEIGHT
130
167
      COMMON/CEVWEIGHT/EVWEIGHT
131
168
      INTEGER IFAIL
132
169
      COMMON/CIFAIL/IFAIL
133
 
 
134
170
      SAVE INOBOSON,INOLEPTON,INOLEPTONB
135
 
 
 
171
      integer nwgt_analysis,max_weight
 
172
      common/c_analysis/nwgt_analysis
 
173
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
174
      double precision ww(max_weight),www(max_weight)
 
175
      common/cww/ww
136
176
      IF(IFAIL.EQ.1)RETURN
137
 
C
138
 
C--DECIDE IDENTITY OF THE VECTOR BOSON, ACCORDING TO THE PDG CODE
139
 
      IDENT=24
 
177
      IF (WW(1).EQ.0D0) THEN
 
178
         WRITE(*,*)'WW(1) = 0. Stopping'
 
179
         STOP
 
180
      ENDIF
 
181
C CHOOSE IDENT = 11 FOR ELECTRON PAIRS
 
182
C        IDENT = 13 FOR MUON PAIRS
 
183
C        IDENT = 15 FOR TAU PAIRS
 
184
      IDENT=13
140
185
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
141
186
C EFFECT, SO THROW THE EVENT AWAY
142
187
      IF(SIGN(1.D0,P(3,3)).EQ.SIGN(1.D0,P(4,3)))THEN
143
188
         WRITE(*,*)'WARNING 111 IN PYANAL'
144
189
         GOTO 999
145
190
      ENDIF
146
 
      WWW0=EVWEIGHT
 
191
      DO I=1,nwgt_analysis
 
192
         WWW(I)=EVWEIGHT*ww(i)/ww(1)
 
193
      ENDDO
147
194
      DO I=1,4
148
195
         P1(I)=P(1,I)
149
196
         P2(I)=P(2,I)
154
201
      KF1=K(1,2)
155
202
      KF2=K(2,2)
156
203
      ICHINI=PYCHGE(KF1)+PYCHGE(KF2)
157
 
      IFV =0
158
 
      IFL =0
159
 
      IFLB=0
160
 
      IV0 =10000
161
 
      IL0 =10000
162
 
      ILB0=10000
163
 
C
 
204
      FOUNDP=.FALSE.
 
205
      FOUNDM=.FALSE.
164
206
      DO 100 IHEP=1,N
165
207
        DO J=1,4
166
208
          PIHEP(J)=P(IHEP,J)
172
214
          CALL VVSUM(4,PIHEP,PSUM,PSUM)
173
215
          ICHSUM=ICHSUM+PYCHGE(ID1)
174
216
        ENDIF
175
 
        TEST1=IORI.EQ.0
176
 
        TEST2=ID1.EQ.IDENT
177
 
        TEST3=ID1.GT.0.AND.ABS(ID1).GE.11.AND.ABS(ID1).LE.16
178
 
        TEST4=ID1.LT.0.AND.ABS(ID1).GE.11.AND.ABS(ID1).LE.16
179
 
        IF(TEST1.AND.TEST2)IV0=IHEP
180
 
        TEST5=IORI.EQ.IV0
181
 
        IF(TEST5)THEN
182
 
           IF(TEST2)THEN
183
 
              IV1=IHEP
184
 
              IFV=IFV+1
185
 
           ELSEIF(TEST3)THEN
186
 
              IL0 =IHEP
187
 
           ELSEIF(TEST4)THEN
188
 
              ILB0=IHEP
189
 
           ENDIF
190
 
        ENDIF
191
 
        TEST6=IORI.EQ.IL0
192
 
        TEST7=IORI.EQ.ILB0
193
 
        IF(TEST6.AND.TEST3)THEN
194
 
           IL =IHEP
195
 
           IFL=IFL+1
196
 
        ENDIF
197
 
        IF(TEST7.AND.TEST4)THEN
198
 
           ILB=IHEP
199
 
           IFLB=IFLB+1
 
217
        ISLP=ID1.EQ. IDENT
 
218
        ISLM=ID1.EQ.-IDENT
 
219
        IF(IORI.LE.10.AND.ISLM)THEN
 
220
           ILL=IHEP
 
221
           FOUNDM=.TRUE.
 
222
        ENDIF
 
223
        IF(IORI.LE.10.AND.ISLP)THEN
 
224
           ILLB=IHEP
 
225
           FOUNDP=.TRUE.
200
226
        ENDIF
201
227
 100  CONTINUE
202
 
C
203
 
      DO IJ=1,5
204
 
        PPV(IJ) =P(IV1,IJ)
205
 
        PPL(IJ) =P(IL, IJ)
206
 
        PPLB(IJ)=P(ILB,IJ)
207
 
      ENDDO
208
 
C CHECK MULTIPLICITIES
209
 
      IF(IFV.EQ.0) THEN
210
 
         INOBOSON=INOBOSON+1
211
 
         WRITE(*,*)'WARNING IN PYANAL: NO WEAK BOSON ',INOBOSON
212
 
         GOTO 999
213
 
      ELSEIF(IFV.GT.1) THEN
214
 
         WRITE(*,*)'WARNING IN PYANAL: TOO MANY WEAK BOSONS ',IFV
215
 
         GOTO 999
216
 
      ENDIF
217
 
      IF(IFL.EQ.0) THEN
218
 
         INOLEPTON=INOLEPTON+1
219
 
         WRITE(*,*)'WARNING IN PYANAL: NO LEPTON ',INOLEPTON
220
 
         GOTO 999
221
 
      ELSEIF(IFL.GT.1) THEN
222
 
         WRITE(*,*)'WARNING IN PYANAL: TOO MANY LEPTONS ',IFL
223
 
         GOTO 999
224
 
      ENDIF
225
 
      IF(IFLB.EQ.0) THEN
226
 
         INOLEPTONB=INOLEPTONB+1
227
 
         WRITE(*,*)'WARNING IN PYANAL: NO ANTILEPTON ',INOLEPTONB
228
 
         GOTO 999
229
 
      ELSEIF(IFLB.GT.1) THEN
230
 
         WRITE(*,*)'WARNING IN PYANAL: TOO MANY ANTILEPTONS ',IFLB
231
 
         GOTO 999
232
 
      ENDIF
233
 
C CHECK THAT THE LEPTONS ARE FINAL-STATE LEPTONS
234
 
      IF( (K(IL,1).NE.1).OR.K(ILB,1).NE.1 ) THEN
235
 
         WRITE(*,*)'WARNING 505 IN PYANAL'
236
 
         GOTO 999
 
228
      IF(.NOT.FOUNDP.OR..NOT.FOUNDM)THEN
 
229
         WRITE(*,*)'NO LEPTONS FOUND.'
 
230
         WRITE(*,*)'CURRENTLY THIS ANALYSIS LOOKS FOR'
 
231
         IF(IDENT.EQ.11)WRITE(*,*)'ELECTRON PAIRS'
 
232
         IF(IDENT.EQ.13)WRITE(*,*)'MUON PAIRS'
 
233
         IF(IDENT.EQ.15)WRITE(*,*)'TAU PAIRS'
 
234
         WRITE(*,*)'IF THIS IS NOT MEANT,'
 
235
         WRITE(*,*)'PLEASE CHANGE THE VALUE OF IDENT IN THIS FILE.'
 
236
         STOP
237
237
      ENDIF
238
238
C CHECK MOMENTUM AND CHARGE CONSERVATION
239
239
      IF (VDOT(3,PSUM,PSUM).GT.1.E-4*P(1,4)**2) THEN
244
244
         WRITE(*,*)'WARNING 113 IN PYANAL'
245
245
         GOTO 999
246
246
      ENDIF
 
247
C FIND THE LEPTONS
 
248
      DO IJ=1,5
 
249
        PPL(IJ) =P(ILL, IJ)
 
250
        PPLB(IJ)=P(ILLB,IJ)
 
251
        PPV(IJ) =PPL(IJ)+PPLB(IJ)
 
252
      ENDDO
247
253
C FILL THE HISTOS
248
254
      YCUT=2.5D0
249
255
C Variables of the vector boson
266
272
      xmll=xmv
267
273
      detallb=etal-etalb
268
274
c
269
 
      kk=0
270
 
      wmass=80.419d0
271
 
      wgamma=2.046d0
272
 
      bwcutoff=15.d0
273
 
      flag=(xmv.ge.wmass-wgamma*bwcutoff.and.
274
 
     &      xmv.le.wmass+wgamma*bwcutoff)
275
 
      if(flag)then
276
 
      call mfill(kk+1,(ptv),(WWW0))
277
 
      call mfill(kk+2,(ptv),(WWW0))
278
 
      if(ptv.gt.0.d0)call mfill(kk+3,(log10(ptv)),(WWW0))
279
 
      call mfill(kk+4,(yv),(WWW0))
280
 
      call mfill(kk+5,(etav),(WWW0))
281
 
      call mfill(kk+6,(xmv),(WWW0))
282
 
c
283
 
      call mfill(kk+7,(ptl),(WWW0))
284
 
      call mfill(kk+8,(ptl),(WWW0))
285
 
      if(ptl.gt.0.d0)call mfill(kk+9,(log10(ptl)),(WWW0))
286
 
      call mfill(kk+10,(etal),(WWW0))
287
 
      call mfill(kk+11,(ptlb),(WWW0))
288
 
      call mfill(kk+12,(ptlb),(WWW0))
289
 
      if(ptlb.gt.0.d0)call mfill(kk+13,(log10(ptlb)),(WWW0))
290
 
      call mfill(kk+14,(etalb),(WWW0))
291
 
c
292
 
      call mfill(kk+15,(detallb),(WWW0))
293
 
      call mfill(kk+16,(azi),(WWW0))
 
275
      do kk=1,nwgt_analysis
 
276
      l=(kk-1)*42
 
277
      call mfill(l+1,(ptv),(WWW(kk)))
 
278
      call mfill(l+2,(ptv),(WWW(kk)))
 
279
      if(ptv.gt.0.d0)call mfill(l+3,(log10(ptv)),(WWW(kk)))
 
280
      call mfill(l+4,(yv),(WWW(kk)))
 
281
      call mfill(l+5,(etav),(WWW(kk)))
 
282
      call mfill(l+6,(xmv),(WWW(kk)))
 
283
c
 
284
      call mfill(l+7,(ptl),(WWW(kk)))
 
285
      call mfill(l+8,(ptl),(WWW(kk)))
 
286
      if(ptl.gt.0.d0)call mfill(l+9,(log10(ptl)),(WWW(kk)))
 
287
      call mfill(l+10,(etal),(WWW(kk)))
 
288
      call mfill(l+11,(ptlb),(WWW(kk)))
 
289
      call mfill(l+12,(ptlb),(WWW(kk)))
 
290
      if(ptlb.gt.0.d0)call mfill(l+13,(log10(ptlb)),(WWW(kk)))
 
291
      call mfill(l+14,(etalb),(WWW(kk)))
 
292
c
 
293
      call mfill(l+15,(detallb),(WWW(kk)))
 
294
      call mfill(l+16,(azi),(WWW(kk)))
294
295
      if(azinorm.gt.0.d0)
295
 
     #  call mfill(kk+17,(log10(azinorm)),(WWW0))
296
 
      call mfill(kk+18,(xmll),(WWW0))
297
 
      call mfill(kk+19,(ptpair),(WWW0))
298
 
      if(ptpair.gt.0)call mfill(kk+20,(log10(ptpair)),(WWW0))
299
 
      call mfill(kk+21,(0d0),(WWW0))
 
296
     #  call mfill(l+17,(log10(azinorm)),(WWW(kk)))
 
297
      call mfill(l+18,(xmll),(WWW(kk)))
 
298
      call mfill(l+19,(ptpair),(WWW(kk)))
 
299
      if(ptpair.gt.0)call mfill(l+20,(log10(ptpair)),(WWW(kk)))
 
300
      call mfill(l+21,(0d0),(WWW(kk)))
300
301
c
301
 
      kk=50
 
302
      l=l+21
 
303
 
302
304
      if(abs(etav).lt.ycut)then
303
 
        call mfill(kk+1,(ptv),(WWW0))
304
 
        call mfill(kk+2,(ptv),(WWW0))
305
 
        if(ptv.gt.0.d0)call mfill(kk+3,(log10(ptv)),(WWW0))
 
305
        call mfill(l+1,(ptv),(WWW(kk)))
 
306
        call mfill(l+2,(ptv),(WWW(kk)))
 
307
        if(ptv.gt.0.d0)call mfill(l+3,(log10(ptv)),(WWW(kk)))
306
308
      endif
307
309
      if(ptv.gt.20.d0)then
308
 
        call mfill(kk+4,(yv),(WWW0))
309
 
        call mfill(kk+5,(etav),(WWW0))
 
310
        call mfill(l+4,(yv),(WWW(kk)))
 
311
        call mfill(l+5,(etav),(WWW(kk)))
310
312
      endif
311
313
      if(abs(etav).lt.ycut.and.ptv.gt.20.d0)then
312
 
         call mfill(kk+6,(xmv),(WWW0))
313
 
         call mfill(kk+21,(0d0),(WWW0))
 
314
         call mfill(l+6,(xmv),(WWW(kk)))
 
315
         call mfill(l+21,(0d0),(WWW(kk)))
314
316
      endif
315
317
c
316
318
      if(abs(etal).lt.ycut)then
317
 
        call mfill(kk+7,(ptl),(WWW0))
318
 
        call mfill(kk+8,(ptl),(WWW0))
319
 
        if(ptl.gt.0.d0)call mfill(kk+9,(log10(ptl)),(WWW0))
 
319
        call mfill(l+7,(ptl),(WWW(kk)))
 
320
        call mfill(l+8,(ptl),(WWW(kk)))
 
321
        if(ptl.gt.0.d0)call mfill(l+9,(log10(ptl)),(WWW(kk)))
320
322
      endif
321
 
      if(ptl.gt.20.d0)call mfill(kk+10,(etal),(WWW0))
 
323
      if(ptl.gt.20.d0)call mfill(l+10,(etal),(WWW(kk)))
322
324
      if(abs(etalb).lt.ycut)then
323
 
        call mfill(kk+11,(ptlb),(WWW0))
324
 
        call mfill(kk+12,(ptlb),(WWW0))
325
 
        if(ptlb.gt.0.d0)call mfill(kk+13,(log10(ptlb)),(WWW0))
 
325
        call mfill(l+11,(ptlb),(WWW(kk)))
 
326
        call mfill(l+12,(ptlb),(WWW(kk)))
 
327
        if(ptlb.gt.0.d0)call mfill(l+13,(log10(ptlb)),(WWW(kk)))
326
328
      endif
327
 
      if(ptlb.gt.20.d0)call mfill(kk+14,(etalb),(WWW0))
 
329
      if(ptlb.gt.20.d0)call mfill(l+14,(etalb),(WWW(kk)))
328
330
c
329
331
      if( abs(etal).lt.ycut.and.abs(etalb).lt.ycut .and.
330
332
     #    ptl.gt.20.d0.and.ptlb.gt.20.d0)then
331
 
        call mfill(kk+15,(detallb),(WWW0))
332
 
        call mfill(kk+16,(azi),(WWW0))
 
333
        call mfill(l+15,(detallb),(WWW(kk)))
 
334
        call mfill(l+16,(azi),(WWW(kk)))
333
335
        if(azinorm.gt.0.d0)
334
 
     #    call mfill(kk+17,(log10(azinorm)),(WWW0))
335
 
        call mfill(kk+18,(xmll),(WWW0))
336
 
        call mfill(kk+19,(ptpair),(WWW0))
 
336
     #    call mfill(l+17,(log10(azinorm)),(WWW(kk)))
 
337
        call mfill(l+18,(xmll),(WWW(kk)))
 
338
        call mfill(l+19,(ptpair),(WWW(kk)))
337
339
        if(ptpair.gt.0) 
338
 
     #    call mfill(kk+20,(log10(ptpair)),(WWW0))
339
 
      endif
340
 
      endif
 
340
     #    call mfill(l+20,(log10(ptpair)),(WWW(kk)))
 
341
      endif
 
342
 
 
343
      enddo
 
344
 
341
345
 999  END
342
346
 
343
 
 
344
347
C-----------------------------------------------------------------------
345
348
      SUBROUTINE VVSUM(N,P,Q,R)
346
349
C-----------------------------------------------------------------------
367
370
   10 Q(I)=C*P(I)
368
371
      END
369
372
 
370
 
 
371
 
 
372
373
C-----------------------------------------------------------------------
373
374
      FUNCTION VDOT(N,P,Q)
374
375
C-----------------------------------------------------------------------