~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/PY8Analyzer/mcatnlo_pyan_pp_ttx_hepmc.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:
 
1
c
 
2
c Example analysis for "p p > t t~ [QCD]" process.
 
3
c
1
4
C----------------------------------------------------------------------
2
5
      SUBROUTINE RCLOS()
3
6
C     DUMMY IF HBOOK IS USED
10
13
C     USER'S ROUTINE FOR INITIALIZATION
11
14
C----------------------------------------------------------------------
12
15
      INCLUDE 'HEPMC.INC'
 
16
      include 'reweight0.inc'
13
17
      REAL*8 pi
14
 
      integer j,k
 
18
      integer j,kk,l,i
15
19
      PARAMETER (PI=3.14159265358979312D0)
16
20
      character*5 cc(2)
17
21
      data cc/'     ',' cuts'/
 
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
18
28
c
19
29
      call inihist
20
 
      do j=1,2
21
 
        k=(j-1)*50
22
 
        call mbook(k+ 1,'tt pt'//cc(j),2.d0,0.d0,100.d0)
23
 
        call mbook(k+ 2,'tt log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
24
 
        call mbook(k+ 3,'tt inv m'//cc(j),10.d0,300.d0,1000.d0)
25
 
        call mbook(k+ 4,'tt azimt'//cc(j),pi/20.d0,0.d0,pi)
26
 
        call mbook(k+ 5,'tt del R'//cc(j),pi/20.d0,0.d0,3*pi)
27
 
        call mbook(k+ 6,'tb pt'//cc(j),5.d0,0.d0,500.d0)
28
 
        call mbook(k+ 7,'tb log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
29
 
        call mbook(k+ 8,'t pt'//cc(j),5.d0,0.d0,500.d0)
30
 
        call mbook(k+ 9,'t log[pt]'//cc(j),0.05d0,0.1d0,5.d0)
31
 
        call mbook(k+10,'tt delta eta'//cc(j),0.2d0,-4.d0,4.d0)
32
 
        call mbook(k+11,'y_tt'//cc(j),0.1d0,-4.d0,4.d0)
33
 
        call mbook(k+12,'delta y'//cc(j),0.2d0,-4.d0,4.d0)
34
 
        call mbook(k+13,'tt azimt'//cc(j),pi/60.d0,2*pi/3,pi)
35
 
        call mbook(k+14,'tt del R'//cc(j),pi/60.d0,2*pi/3,4*pi/3)
36
 
        call mbook(k+15,'y_tb'//cc(j),0.1d0,-4.d0,4.d0)
37
 
        call mbook(k+16,'y_t'//cc(j),0.1d0,-4.d0,4.d0)
38
 
        call mbook(k+17,'tt log[pi-azimt]'//cc(j),0.05d0,-4.d0,0.1d0)
 
30
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
31
c To be changed !!
 
32
      nwgt=1
 
33
      weights_info(nwgt)="central value  "
 
34
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
35
      nwgt_analysis=nwgt
 
36
      do kk=1,nwgt_analysis
 
37
      do i=1,2
 
38
        l=(kk-1)*40+(i-1)*20
 
39
        call mbook(l+ 1,'tt pt            '
 
40
     &       //weights_info(kk)//cc(i),2.d0,0.d0,100.d0)
 
41
        call mbook(l+ 2,'tt log[pt]       '
 
42
     &       //weights_info(kk)//cc(i),0.05d0,0.1d0,5.d0)
 
43
        call mbook(l+ 3,'tt inv m         '
 
44
     &       //weights_info(kk)//cc(i),10.d0,300.d0,1000.d0)
 
45
        call mbook(l+ 4,'tt azimt         '
 
46
     &       //weights_info(kk)//cc(i),pi/20.d0,0.d0,pi)
 
47
        call mbook(l+ 5,'tt del R         '
 
48
     &       //weights_info(kk)//cc(i),pi/20.d0,0.d0,3*pi)
 
49
        call mbook(l+ 6,'tb pt            '
 
50
     &       //weights_info(kk)//cc(i),5.d0,0.d0,500.d0)
 
51
        call mbook(l+ 7,'tb log[pt]       '
 
52
     &       //weights_info(kk)//cc(i),0.05d0,0.1d0,5.d0)
 
53
        call mbook(l+ 8,'t pt             '
 
54
     &       //weights_info(kk)//cc(i),5.d0,0.d0,500.d0)
 
55
        call mbook(l+ 9,'t log[pt]        '
 
56
     &       //weights_info(kk)//cc(i),0.05d0,0.1d0,5.d0)
 
57
        call mbook(l+10,'tt delta eta     '
 
58
     &       //weights_info(kk)//cc(i),0.2d0,-4.d0,4.d0)
 
59
        call mbook(l+11,'y_tt             '
 
60
     &       //weights_info(kk)//cc(i),0.1d0,-4.d0,4.d0)
 
61
        call mbook(l+12,'delta y          '
 
62
     &       //weights_info(kk)//cc(i),0.2d0,-4.d0,4.d0)
 
63
        call mbook(l+13,'tt azimt         '
 
64
     &       //weights_info(kk)//cc(i),pi/60.d0,2*pi/3,pi)
 
65
        call mbook(l+14,'tt del R         '
 
66
     &       //weights_info(kk)//cc(i),pi/60.d0,2*pi/3,4*pi/3)
 
67
        call mbook(l+15,'y_tb             '
 
68
     &       //weights_info(kk)//cc(i),0.1d0,-4.d0,4.d0)
 
69
        call mbook(l+16,'y_t              '
 
70
     &       //weights_info(kk)//cc(i),0.1d0,-4.d0,4.d0)
 
71
        call mbook(l+17,'tt log[pi-azimt] '
 
72
     &       //weights_info(kk)//cc(i),0.05d0,-4.d0,0.1d0)
 
73
        call mbook(l+18,'tt pt            '
 
74
     &       //weights_info(kk)//cc(i),20.d0,80.d0,2000.d0)
 
75
        call mbook(l+19,'tb pt            '
 
76
     &       //weights_info(kk)//cc(i),20.d0,400.d0,2400.d0)
 
77
        call mbook(l+20,'t pt             '
 
78
     &       //weights_info(kk)//cc(i),20.d0,400.d0,2400.d0)
39
79
      enddo
40
 
      do j=1,2
41
 
        k=(j-1)*50
42
 
        call mbook(k+18,'tt pt'//cc(j),20.d0,80.d0,2000.d0)
43
 
        call mbook(k+19,'tb pt'//cc(j),20.d0,400.d0,2400.d0)
44
 
        call mbook(k+20,'t pt'//cc(j),20.d0,400.d0,2400.d0)
45
80
      enddo
46
81
      END
47
82
 
52
87
C----------------------------------------------------------------------
53
88
      INCLUDE 'HEPMC.INC'
54
89
      REAL*8 XNORM
55
 
      INTEGER I,J,K
 
90
      INTEGER I,J,KK,l,nwgt_analysis
 
91
      integer NPL
 
92
      parameter(NPL=15000)
 
93
      common/c_analysis/nwgt_analysis
56
94
      OPEN(UNIT=99,FILE='PYTQQ.TOP',STATUS='UNKNOWN')
57
95
C XNORM IS SUCH THAT THE CROSS SECTION PER BIN IS IN PB, SINCE THE HERWIG 
58
96
C WEIGHT IS IN NB, AND CORRESPONDS TO THE AVERAGE CROSS SECTION
59
97
      XNORM=IEVTTOT/DFLOAT(NEVHEP)
60
 
      DO I=1,100              
 
98
      DO I=1,NPL
61
99
        CALL MFINAL3(I)             
62
 
        CALL MCOPY(I,I+100)
63
 
        CALL MOPERA(I+100,'F',I+100,I+100,(XNORM),0.D0)
64
 
        CALL MFINAL3(I+100)             
 
100
        CALL MCOPY(I,I+NPL)
 
101
        CALL MOPERA(I+NPL,'F',I+NPL,I+NPL,(XNORM),0.D0)
 
102
        CALL MFINAL3(I+NPL)             
65
103
      ENDDO                          
66
104
C
67
 
      do j=1,2
68
 
        k=(j-1)*50
69
 
        call multitop(100+k+ 1,99,2,3,'tt pt',' ','LOG')
70
 
        call multitop(100+k+ 2,99,2,3,'tt log[pt]',' ','LOG')
71
 
        call multitop(100+k+ 3,99,2,3,'tt inv m',' ','LOG')
72
 
        call multitop(100+k+ 4,99,2,3,'tt azimt',' ','LOG')
73
 
        call multitop(100+k+ 5,99,2,3,'tt del R',' ','LOG')
74
 
        call multitop(100+k+ 6,99,2,3,'tb pt',' ','LOG')
75
 
        call multitop(100+k+ 7,99,2,3,'tb log[pt]',' ','LOG')
76
 
        call multitop(100+k+ 8,99,2,3,'t pt',' ','LOG')
77
 
        call multitop(100+k+ 9,99,2,3,'t log[pt]',' ','LOG')
78
 
        call multitop(100+k+10,99,2,3,'tt Delta eta',' ','LOG')
79
 
        call multitop(100+k+11,99,2,3,'y_tt',' ','LOG')
80
 
        call multitop(100+k+12,99,2,3,'tt Delta y',' ','LOG')
81
 
        call multitop(100+k+13,99,2,3,'tt azimt',' ','LOG')
82
 
        call multitop(100+k+14,99,2,3,'tt del R',' ','LOG')
83
 
        call multitop(100+k+15,99,2,3,'tb y',' ','LOG')
84
 
        call multitop(100+k+16,99,2,3,'t y',' ','LOG')
85
 
        call multitop(100+k+17,99,2,3,'tt log[pi-azimt]',' ','LOG')
86
 
      enddo
87
 
      do j=1,2
88
 
        k=(j-1)*50
89
 
        call multitop(100+k+18,99,2,3,'tt pt',' ','LOG')
90
 
        call multitop(100+k+19,99,2,3,'tb pt',' ','LOG')
91
 
        call multitop(100+k+20,99,2,3,'t pt',' ','LOG')
92
 
      enddo
93
 
c
 
105
      do kk=1,nwgt_analysis
 
106
      do i=1,2
 
107
        l=(kk-1)*40+(i-1)*20
 
108
        call multitop(NPL+l+ 1,NPL-1,2,3,'tt pt',' ','LOG')
 
109
        call multitop(NPL+l+ 2,NPL-1,2,3,'tt log[pt]',' ','LOG')
 
110
        call multitop(NPL+l+ 3,NPL-1,2,3,'tt inv m',' ','LOG')
 
111
        call multitop(NPL+l+ 4,NPL-1,2,3,'tt azimt',' ','LOG')
 
112
        call multitop(NPL+l+ 5,NPL-1,2,3,'tt del R',' ','LOG')
 
113
        call multitop(NPL+l+ 6,NPL-1,2,3,'tb pt',' ','LOG')
 
114
        call multitop(NPL+l+ 7,NPL-1,2,3,'tb log[pt]',' ','LOG')
 
115
        call multitop(NPL+l+ 8,NPL-1,2,3,'t pt',' ','LOG')
 
116
        call multitop(NPL+l+ 9,NPL-1,2,3,'t log[pt]',' ','LOG')
 
117
        call multitop(NPL+l+10,NPL-1,2,3,'tt Delta eta',' ','LOG')
 
118
        call multitop(NPL+l+11,NPL-1,2,3,'y_tt',' ','LOG')
 
119
        call multitop(NPL+l+12,NPL-1,2,3,'tt Delta y',' ','LOG')
 
120
        call multitop(NPL+l+13,NPL-1,2,3,'tt azimt',' ','LOG')
 
121
        call multitop(NPL+l+14,NPL-1,2,3,'tt del R',' ','LOG')
 
122
        call multitop(NPL+l+15,NPL-1,2,3,'tb y',' ','LOG')
 
123
        call multitop(NPL+l+16,NPL-1,2,3,'t y',' ','LOG')
 
124
        call multitop(NPL+l+17,NPL-1,2,3,'tt log[pi-azimt]',' ','LOG')
 
125
        call multitop(NPL+l+18,NPL-1,2,3,'tt pt',' ','LOG')
 
126
        call multitop(NPL+l+19,NPL-1,2,3,'tb pt',' ','LOG')
 
127
        call multitop(NPL+l+20,NPL-1,2,3,'t pt',' ','LOG')
 
128
      enddo
 
129
      enddo
94
130
      CLOSE(99)
95
131
      END
96
132
 
100
136
C     USER'S ROUTINE TO ANALYSE DATA FROM EVENT
101
137
C----------------------------------------------------------------------
102
138
      INCLUDE 'HEPMC.INC'
 
139
      include 'reweight0.inc'
103
140
      DOUBLE PRECISION HWVDOT,PSUM(4)
104
141
      INTEGER ICHSUM,ICHINI,IHEP
105
142
      LOGICAL DIDSOF,flcuts,siq1flag,siq2flag,ddflag
117
154
      REAL*8 PI
118
155
      PARAMETER (PI=3.14159265358979312D0)
119
156
      REAL*8 WWW0
120
 
      INTEGER KK,IVLEP1,IVLEP2
 
157
      INTEGER KK,IVLEP1,IVLEP2,i,l
121
158
      COMMON/VVLIN/IVLEP1,IVLEP2
 
159
      integer nwgt_analysis,max_weight
 
160
      common/c_analysis/nwgt_analysis
 
161
      parameter (max_weight=maxscales*maxscales+maxpdfs+1)
 
162
      double precision ww(max_weight),www(max_weight)
 
163
      common/cww/ww
 
164
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
165
c To be changed !!
 
166
      ww(1)=1d0
 
167
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
122
168
c
123
169
      IF(MOD(NEVHEP,10000).EQ.0)RETURN
124
 
c
 
170
      IF (WW(1).EQ.0D0) THEN
 
171
         WRITE(*,*)'WW(1) = 0. Stopping'
 
172
         STOP
 
173
      ENDIF
125
174
C INCOMING PARTONS MAY TRAVEL IN THE SAME DIRECTION: IT'S A POWER-SUPPRESSED
126
175
C EFFECT, SO THROW THE EVENT AWAY
127
176
 
131
180
        WRITE(*,*)PHEP(3,1),PHEP(3,2)
132
181
        GOTO 999
133
182
      ENDIF
134
 
      WWW0=EVWGT
 
183
      DO I=1,nwgt_analysis
 
184
         WWW(I)=EVWGT*ww(i)/ww(1)
 
185
      ENDDO
135
186
      ICHSUM=0
136
187
      DIDSOF=.FALSE.
137
188
      IQ1=0
184
235
      siq2flag=ptq2.gt.ptcut.and.abs(yq2).lt.ycut
185
236
      ddflag=siq1flag.and.siq2flag
186
237
c-------------------------------------------------------------
187
 
      call mfill(1,(ptg),(WWW0))
188
 
      call mfill(18,(ptg),(WWW0))
189
 
      if(ptg.gt.0) call mfill(2,(log10(ptg)),(WWW0))
190
 
      call mfill(3,(qqm),(WWW0))
191
 
      call mfill(4,(azi),(WWW0))
192
 
      call mfill(13,(azi),(WWW0))
193
 
      if(azinorm.gt.0)call mfill(17,(log10(azinorm)),(WWW0))
194
 
      call mfill(5,(dr),(WWW0))
195
 
      call mfill(14,(dr),(WWW0))
196
 
      call mfill(10,(etaq1-etaq2),(WWW0))
197
 
      call mfill(11,(yqq),(WWW0))
198
 
      call mfill(12,(yq1-yq2),(WWW0))
199
 
      call mfill(6,(ptq2),(WWW0))
200
 
      call mfill(19,(ptq2),(WWW0))
201
 
      if(ptq2.gt.0) call mfill(7,(log10(ptq2)),(WWW0))
202
 
      call mfill(15,(yq2),(WWW0))
203
 
      call mfill(8,(ptq1),(WWW0))
204
 
      call mfill(20,(ptq1),(WWW0))
205
 
      if(ptq1.gt.0) call mfill(9,(log10(ptq1)),(WWW0))
206
 
      call mfill(16,(yq1),(WWW0))
 
238
      do kk=1,nwgt_analysis
 
239
         l=(kk-1)*40
 
240
         call mfill(l+1,ptg,WWW(kk))
 
241
         call mfill(l+18,ptg,WWW(kk))
 
242
         if(ptg.gt.0) call mfill(l+2,log10(ptg),WWW(kk))
 
243
         call mfill(l+3,qqm,WWW(kk))
 
244
         call mfill(l+4,azi,WWW(kk))
 
245
         call mfill(l+13,azi,WWW(kk))
 
246
         if(azinorm.gt.0)call mfill(l+17,log10(azinorm),WWW(kk))
 
247
         call mfill(l+5,dr,WWW(kk))
 
248
         call mfill(l+14,dr,WWW(kk))
 
249
         call mfill(l+10,etaq1-etaq2,WWW(kk))
 
250
         call mfill(l+11,yqq,WWW(kk))
 
251
         call mfill(l+12,yq1-yq2,WWW(kk))
 
252
         call mfill(l+6,ptq2,WWW(kk))
 
253
         call mfill(l+19,ptq2,WWW(kk))
 
254
         if(ptq2.gt.0) call mfill(l+7,log10(ptq2),WWW(kk))
 
255
         call mfill(l+15,yq2,WWW(kk))
 
256
         call mfill(l+8,ptq1,WWW(kk))
 
257
         call mfill(l+20,ptq1,WWW(kk))
 
258
         if(ptq1.gt.0) call mfill(l+9,log10(ptq1),WWW(kk))
 
259
         call mfill(l+16,yq1,WWW(kk))
207
260
c
208
261
c***************************************************** with cuts
209
262
c
210
 
      kk=50
211
 
      if(ddflag)then
212
 
        call mfill(kk+1,(ptg),(WWW0))
213
 
        call mfill(kk+18,(ptg),(WWW0))
214
 
        if(ptg.gt.0) call mfill(kk+2,(log10(ptg)),(WWW0))
215
 
        call mfill(kk+3,(qqm),(WWW0))
216
 
        call mfill(kk+4,(azi),(WWW0))
217
 
        call mfill(kk+13,(azi),(WWW0))
218
 
        if(azinorm.gt.0) 
219
 
     #    call mfill(kk+17,(log10(azinorm)),(WWW0))
220
 
        call mfill(kk+5,(dr),(WWW0))
221
 
        call mfill(kk+14,(dr),(WWW0))
222
 
        call mfill(kk+10,(etaq1-etaq2),(WWW0))
223
 
        call mfill(kk+11,(yqq),(WWW0))
224
 
        call mfill(kk+12,(yq1-yq2),(WWW0))
225
 
      endif
226
 
      if(abs(yq2).lt.ycut)then
227
 
        call mfill(kk+6,(ptq2),(WWW0))
228
 
        call mfill(kk+19,(ptq2),(WWW0))
229
 
        if(ptq2.gt.0) call mfill(kk+7,(log10(ptq2)),(WWW0))
230
 
      endif
231
 
      if(ptq2.gt.ptcut)call mfill(kk+15,(yq2),(WWW0))
232
 
      if(abs(yq1).lt.ycut)then
233
 
        call mfill(kk+8,(ptq1),(WWW0))
234
 
        call mfill(kk+20,(ptq1),(WWW0))
235
 
        if(ptq1.gt.0) call mfill(kk+9,(log10(ptq1)),(WWW0))
236
 
      endif
237
 
      if(ptq1.gt.ptcut)call mfill(kk+16,(yq1),(WWW0))
238
 
 
 
263
         l=l+20
 
264
c
 
265
         if(ddflag)then
 
266
            call mfill(l+1,ptg,WWW(kk))
 
267
            call mfill(l+18,ptg,WWW(kk))
 
268
            if(ptg.gt.0) call mfill(l+2,log10(ptg),WWW(kk))
 
269
            call mfill(l+3,qqm,WWW(kk))
 
270
            call mfill(l+4,azi,WWW(kk))
 
271
            call mfill(l+13,azi,WWW(kk))
 
272
            if(azinorm.gt.0) call mfill(l+17,log10(azinorm),WWW(kk))
 
273
            call mfill(l+5,dr,WWW(kk))
 
274
            call mfill(l+14,dr,WWW(kk))
 
275
            call mfill(l+10,etaq1-etaq2,WWW(kk))
 
276
            call mfill(l+11,yqq,WWW(kk))
 
277
            call mfill(l+12,yq1-yq2,WWW(kk))
 
278
         endif
 
279
         if(abs(yq2).lt.ycut)then
 
280
            call mfill(l+6,ptq2,WWW(kk))
 
281
            call mfill(l+19,ptq2,WWW(kk))
 
282
            if(ptq2.gt.0) call mfill(l+7,log10(ptq2),WWW(kk))
 
283
         endif
 
284
         if(ptq2.gt.ptcut)call mfill(l+15,yq2,WWW(kk))
 
285
         if(abs(yq1).lt.ycut)then
 
286
            call mfill(l+8,ptq1,WWW(kk))
 
287
            call mfill(l+20,ptq1,WWW(kk))
 
288
            if(ptq1.gt.0) call mfill(l+9,log10(ptq1),WWW(kk))
 
289
         endif
 
290
         if(ptq1.gt.ptcut)call mfill(l+16,yq1,WWW(kk))
 
291
      enddo
239
292
 999  return
240
293
      end
241
294