← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to Template/NLO/FixedOrderAnalysis/analysis_td_pp_lplm.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/FixedOrderAnalysis/analysis_td_pp_lplm.f
 
1
c
 
2
c Example analysis for "p p > l+ l- [QCD]" process.
 
3
c
 
4
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
5
      subroutine analysis_begin(nwgt,weights_info)
 
6
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
7
      implicit none
 
8
      integer nwgt
 
9
      character*(*) weights_info(*)
 
10
      integer j,kk,l,nwgt_analysis
 
11
      common/c_analysis/nwgt_analysis
 
12
      character*5 cc(2)
 
13
      data cc/'     ',' cuts'/
 
14
      real * 8 bin,xmi,xms,pi
 
15
      parameter (pi=3.14159265358979312d0)
 
16
      include 'dbook.inc'
 
17
      call inihist
 
18
      nwgt_analysis=nwgt
 
19
      if (nwgt_analysis*42.gt.nplots/4) then
 
20
         write (*,*) 'error in analysis_begin: '/
 
21
     &        /'too many histograms, increase NPLOTS to',
 
22
     &        nwgt_analysis*42*4
 
23
         stop 1
 
24
      endif
 
25
      xmi=50.d0
 
26
      xms=130.d0
 
27
      bin=0.8d0
 
28
      do kk=1,nwgt_analysis
 
29
      do j=1,2
 
30
      l=(kk-1)*42+(j-1)*21
 
31
      call bookup(l+ 1,'V pt      '//weights_info(kk)//cc(j)
 
32
     &     ,2.d0,0.d0,200.d0)
 
33
      call bookup(l+ 2,'V pt      '//weights_info(kk)//cc(j)
 
34
     &     ,10.d0,0.d0,1000.d0)
 
35
      call bookup(l+ 3,'V log[pt] '//weights_info(kk)//cc(j)
 
36
     &     ,0.05d0,0.1d0,5.d0)
 
37
      call bookup(l+ 4,'V y       '//weights_info(kk)//cc(j)
 
38
     &     ,0.2d0,-9.d0,9.d0)
 
39
      call bookup(l+ 5,'V eta     '//weights_info(kk)//cc(j)
 
40
     &     ,0.2d0,-9.d0,9.d0)
 
41
      call bookup(l+ 6,'mV        '//weights_info(kk)//cc(j)
 
42
     &     ,bin,xmi,xms)
 
43
c
 
44
      call bookup(l+ 7,'lm pt      '//weights_info(kk)//cc(j)
 
45
     &     ,2.d0,0.d0,200.d0)
 
46
      call bookup(l+ 8,'lm pt      '//weights_info(kk)//cc(j)
 
47
     &     ,10.d0,0.d0,1000.d0)
 
48
      call bookup(l+ 9,'lm log[pt] '//weights_info(kk)//cc(j)
 
49
     &     ,0.05d0,0.1d0,5.d0)
 
50
      call bookup(l+10,'lm eta     '//weights_info(kk)//cc(j)
 
51
     &     ,0.2d0,-9.d0,9.d0)
 
52
      call bookup(l+11,'lp pt      '//weights_info(kk)//cc(j)
 
53
     &     ,2.d0,0.d0,200.d0)
 
54
      call bookup(l+12,'lp pt      '//weights_info(kk)//cc(j)
 
55
     &     ,10.d0,0.d0,1000.d0)
 
56
      call bookup(l+13,'lp log[pt] '//weights_info(kk)//cc(j)
 
57
     &     ,0.05d0,0.1d0,5.d0)
 
58
      call bookup(l+14,'lp eta     '//weights_info(kk)//cc(j)
 
59
     &     ,0.2d0,-9.d0,9.d0)
 
60
c
 
61
      call bookup(l+15,'lmlp delta eta     '//weights_info(kk)//cc(j)
 
62
     $     ,0.2d0,-9.d0,9.d0)
 
63
      call bookup(l+16,'lmlp azimt         '//weights_info(kk)//cc(j)
 
64
     $     ,pi/20.d0,0.d0,pi)
 
65
      call bookup(l+17,'lmlp log[pi-azimt] '//weights_info(kk)//cc(j)
 
66
     $     ,0.05d0,-4.d0,0.1d0)
 
67
      call bookup(l+18,'lmlp inv m         '//weights_info(kk)//cc(j)
 
68
     $     ,bin,xmi,xms)
 
69
      call bookup(l+19,'lmlp pt            '//weights_info(kk)//cc(j)
 
70
     $     ,2.d0,0.d0,200.d0)
 
71
      call bookup(l+20,'lmlp log[pt]       '//weights_info(kk)//cc(j)
 
72
     $     ,0.05d0,0.1d0,5.d0)
 
73
c
 
74
      call bookup(l+21,'total'//weights_info(kk)//cc(j),1.d0,-1.d0,1.d0)
 
75
      enddo
 
76
      enddo
 
77
      return
 
78
      end
 
79
 
 
80
 
 
81
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
82
      subroutine analysis_end(xnorm)
 
83
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
84
      implicit none
 
85
      character*14 ytit
 
86
      double precision xnorm
 
87
      integer i
 
88
      integer kk,l,nwgt_analysis
 
89
      common/c_analysis/nwgt_analysis
 
90
      include 'dbook.inc'
 
91
      call open_topdrawer_file
 
92
      call mclear
 
93
      do i=1,NPLOTS
 
94
         call mopera(i,'+',i,i,xnorm,0.d0)
 
95
         call mfinal(i)
 
96
      enddo
 
97
      ytit='sigma per bin '
 
98
      do kk=1,nwgt_analysis
 
99
      do i=1,2
 
100
      l=(kk-1)*42+(i-1)*21
 
101
C
 
102
      call multitop(l+ 1,3,2,'V pt',' ','LOG')
 
103
      call multitop(l+ 2,3,2,'V pt',' ','LOG')
 
104
      call multitop(l+ 3,3,2,'V log[pt]',' ','LOG')
 
105
      call multitop(l+ 4,3,2,'V y',' ','LOG')
 
106
      call multitop(l+ 5,3,2,'V eta',' ','LOG')
 
107
      call multitop(l+ 6,3,2,'mV',' ','LOG')
 
108
c
 
109
      call multitop(l+ 7,3,2,'lm pt',' ','LOG')
 
110
      call multitop(l+ 8,3,2,'lm pt',' ','LOG')
 
111
      call multitop(l+ 9,3,2,'lm log[pt]',' ','LOG')
 
112
      call multitop(l+10,3,2,'lm eta',' ','LOG')
 
113
      call multitop(l+11,3,2,'lm pt',' ','LOG')
 
114
      call multitop(l+12,3,2,'lm pt',' ','LOG')
 
115
      call multitop(l+13,3,2,'lm log[pt]',' ','LOG')
 
116
      call multitop(l+14,3,2,'lm eta',' ','LOG')
 
117
c
 
118
      call multitop(l+15,3,2,'lmlp deta',' ','LOG')
 
119
      call multitop(l+16,3,2,'lmlp azi',' ','LOG')
 
120
      call multitop(l+17,3,2,'lmlp azi',' ','LOG')
 
121
      call multitop(l+18,3,2,'lmlp inv m',' ','LOG')
 
122
      call multitop(l+19,3,2,'lmlp pt',' ','LOG')
 
123
      call multitop(l+20,3,2,'lmlp pt',' ','LOG')
 
124
c
 
125
      call multitop(l+21,3,2,'total',' ','LOG')
 
126
      enddo
 
127
      enddo
 
128
c
 
129
      call close_topdrawer_file
 
130
      return                
 
131
      end
 
132
 
 
133
 
 
134
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
135
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
 
136
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
137
      implicit none
 
138
      include 'nexternal.inc'
 
139
      integer istatus(nexternal)
 
140
      integer iPDG(nexternal)
 
141
      double precision p(0:4,nexternal)
 
142
      double precision wgts(*)
 
143
      integer ibody
 
144
      double precision wgt,var
 
145
      integer i,kk,l,nwgt_analysis
 
146
      common/c_analysis/nwgt_analysis
 
147
      double precision www,ppl(0:3),pplb(0:3),ppv(0:3),ycut,xmv,ptv,yv
 
148
     $     ,etav,ptl,yl,etal,ptlb,ylb,etalb,ptpair,azi,azinorm,xmll
 
149
     $     ,detallb
 
150
      double precision getrapidity,getinvm,getpseudorap,getdelphi
 
151
      external getrapidity,getinvm,getpseudorap,getdelphi
 
152
      double precision pi
 
153
      parameter (pi=3.14159265358979312d0)
 
154
      if (nexternal.ne.5) then
 
155
         write (*,*) 'error #1 in analysis_fill: '/
 
156
     &        /'only for process "p p > l+ l- [QCD]"'
 
157
         stop 1
 
158
      endif
 
159
      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
 
160
         write (*,*) 'error #2 in analysis_fill: '/
 
161
     &        /'only for process "p p > l+ l- [QCD]"'
 
162
         stop 1
 
163
      endif
 
164
      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
 
165
         write (*,*) 'error #3 in analysis_fill: '/
 
166
     &        /'only for process "p p > l+ l- [QCD]"'
 
167
         stop 1
 
168
      endif
 
169
      if (.not. (abs(ipdg(5)).le.5 .or. ipdg(5).eq.21)) then
 
170
         write (*,*) 'error #4 in analysis_fill: '/
 
171
     &        /'only for process "p p > l+ l- [QCD]"'
 
172
         stop 1
 
173
      endif
 
174
      if (ipdg(3).ne.-11.and.ipdg(3).ne.-13.and.ipdg(3).ne.-15) then
 
175
         write (*,*) 'error #5 in analysis_fill: '/
 
176
     &        /'only for process "p p > l+ l- [QCD]"'
 
177
         stop 1
 
178
      endif
 
179
      if (ipdg(4).ne.11.and.ipdg(4).ne.13.and.ipdg(4).ne.15) then
 
180
         write (*,*) 'error #6 in analysis_fill: '/
 
181
     &        /'only for process "p p > l+ l- [QCD]"'
 
182
         stop 1
 
183
      endif
 
184
 
 
185
      DO i=0,3
 
186
        ppl(i)=p(i,4)
 
187
        pplb(i)=p(i,3)
 
188
        ppv(i)=ppl(i)+pplb(i)
 
189
      ENDDO
 
190
 
 
191
C FILL THE HISTOS
 
192
      YCUT=2.5D0
 
193
C Variables of the vector boson
 
194
      xmv=getinvm(ppv(0),ppv(1),ppv(2),ppv(3))
 
195
      ptv=sqrt(ppv(1)**2+ppv(2)**2)
 
196
      yv=getrapidity(ppv(0),ppv(3))
 
197
      etav=getpseudorap(ppv(0),ppv(1),ppv(2),ppv(3))
 
198
C Variables of the leptons
 
199
      ptl=sqrt(ppl(1)**2+ppl(2)**2)
 
200
      yl=getrapidity(ppl(0),ppl(3))
 
201
      etal=getpseudorap(ppl(0),ppl(1),ppl(2),ppl(3))
 
202
c
 
203
      ptlb=sqrt(pplb(1)**2+pplb(2)**2)
 
204
      ylb=getrapidity(pplb(0),pplb(3))
 
205
      etalb=getpseudorap(pplb(0),pplb(1),pplb(2),pplb(3))
 
206
c
 
207
      ptpair=ptv
 
208
      azi=getdelphi(ppl(1),pplb(1),ppl(2),pplb(2))
 
209
      azinorm=(pi-azi)/pi
 
210
      xmll=xmv
 
211
      detallb=etal-etalb
 
212
c
 
213
      do kk=1,nwgt_analysis
 
214
      www=wgts(kk)
 
215
      l=(kk-1)*42
 
216
      call mfill(l+1,(ptv),(WWW))
 
217
      call mfill(l+2,(ptv),(WWW))
 
218
      if(ptv.gt.0.d0)call mfill(l+3,(log10(ptv)),(WWW))
 
219
      call mfill(l+4,(yv),(WWW))
 
220
      call mfill(l+5,(etav),(WWW))
 
221
      call mfill(l+6,(xmv),(WWW))
 
222
c
 
223
      call mfill(l+7,(ptl),(WWW))
 
224
      call mfill(l+8,(ptl),(WWW))
 
225
      if(ptl.gt.0.d0)call mfill(l+9,(log10(ptl)),(WWW))
 
226
      call mfill(l+10,(etal),(WWW))
 
227
      call mfill(l+11,(ptlb),(WWW))
 
228
      call mfill(l+12,(ptlb),(WWW))
 
229
      if(ptlb.gt.0.d0)call mfill(l+13,(log10(ptlb)),(WWW))
 
230
      call mfill(l+14,(etalb),(WWW))
 
231
c
 
232
      call mfill(l+15,(detallb),(WWW))
 
233
      call mfill(l+16,(azi),(WWW))
 
234
      if(azinorm.gt.0.d0)
 
235
     #  call mfill(l+17,(log10(azinorm)),(WWW))
 
236
      call mfill(l+18,(xmll),(WWW))
 
237
      call mfill(l+19,(ptpair),(WWW))
 
238
      if(ptpair.gt.0)call mfill(l+20,(log10(ptpair)),(WWW))
 
239
      call mfill(l+21,(0d0),(WWW))
 
240
c
 
241
      l=l+21
 
242
 
 
243
      if(abs(etav).lt.ycut)then
 
244
        call mfill(l+1,(ptv),(WWW))
 
245
        call mfill(l+2,(ptv),(WWW))
 
246
        if(ptv.gt.0.d0)call mfill(l+3,(log10(ptv)),(WWW))
 
247
      endif
 
248
      if(ptv.gt.20.d0)then
 
249
        call mfill(l+4,(yv),(WWW))
 
250
        call mfill(l+5,(etav),(WWW))
 
251
      endif
 
252
      if(abs(etav).lt.ycut.and.ptv.gt.20.d0)then
 
253
         call mfill(l+6,(xmv),(WWW))
 
254
         call mfill(l+21,(0d0),(WWW))
 
255
      endif
 
256
c
 
257
      if(abs(etal).lt.ycut)then
 
258
        call mfill(l+7,(ptl),(WWW))
 
259
        call mfill(l+8,(ptl),(WWW))
 
260
        if(ptl.gt.0.d0)call mfill(l+9,(log10(ptl)),(WWW))
 
261
      endif
 
262
      if(ptl.gt.20.d0)call mfill(l+10,(etal),(WWW))
 
263
      if(abs(etalb).lt.ycut)then
 
264
        call mfill(l+11,(ptlb),(WWW))
 
265
        call mfill(l+12,(ptlb),(WWW))
 
266
        if(ptlb.gt.0.d0)call mfill(l+13,(log10(ptlb)),(WWW))
 
267
      endif
 
268
      if(ptlb.gt.20.d0)call mfill(l+14,(etalb),(WWW))
 
269
c
 
270
      if( abs(etal).lt.ycut.and.abs(etalb).lt.ycut .and.
 
271
     #    ptl.gt.20.d0.and.ptlb.gt.20.d0)then
 
272
        call mfill(l+15,(detallb),(WWW))
 
273
        call mfill(l+16,(azi),(WWW))
 
274
        if(azinorm.gt.0.d0)
 
275
     #    call mfill(l+17,(log10(azinorm)),(WWW))
 
276
        call mfill(l+18,(xmll),(WWW))
 
277
        call mfill(l+19,(ptpair),(WWW))
 
278
        if(ptpair.gt.0) 
 
279
     #    call mfill(l+20,(log10(ptpair)),(WWW))
 
280
      endif
 
281
 
 
282
      enddo
 
283
 
 
284
 999  return      
 
285
      end
 
286
 
 
287
 
 
288
      function getrapidity(en,pl)
 
289
      implicit none
 
290
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
291
      parameter (tiny=1.d-8)
 
292
      xplus=en+pl
 
293
      xminus=en-pl
 
294
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
295
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
296
            y=0.5d0*log( xplus/xminus  )
 
297
         else
 
298
            y=sign(1.d0,pl)*1.d8
 
299
         endif
 
300
      else 
 
301
         y=sign(1.d0,pl)*1.d8
 
302
      endif
 
303
      getrapidity=y
 
304
      return
 
305
      end
 
306
 
 
307
 
 
308
      function getpseudorap(en,ptx,pty,pl)
 
309
      implicit none
 
310
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
311
      parameter (tiny=1.d-5)
 
312
c
 
313
      pt=sqrt(ptx**2+pty**2)
 
314
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
315
        eta=sign(1.d0,pl)*1.d8
 
316
      else
 
317
        th=atan2(pt,pl)
 
318
        eta=-log(tan(th/2.d0))
 
319
      endif
 
320
      getpseudorap=eta
 
321
      return
 
322
      end
 
323
 
 
324
 
 
325
      function getinvm(en,ptx,pty,pl)
 
326
      implicit none
 
327
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
328
      parameter (tiny=1.d-5)
 
329
c
 
330
      tmp=en**2-ptx**2-pty**2-pl**2
 
331
      if(tmp.gt.0.d0)then
 
332
        tmp=sqrt(tmp)
 
333
      elseif(tmp.gt.-tiny)then
 
334
        tmp=0.d0
 
335
      else
 
336
        write(*,*)'Attempt to compute a negative mass'
 
337
        stop
 
338
      endif
 
339
      getinvm=tmp
 
340
      return
 
341
      end
 
342
 
 
343
 
 
344
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
345
      implicit none
 
346
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
347
      parameter (tiny=1.d-5)
 
348
c
 
349
      pt1=sqrt(ptx1**2+pty1**2)
 
350
      pt2=sqrt(ptx2**2+pty2**2)
 
351
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
352
        tmp=ptx1*ptx2+pty1*pty2
 
353
        tmp=tmp/(pt1*pt2)
 
354
        if(abs(tmp).gt.1.d0+tiny)then
 
355
          write(*,*)'Cosine larger than 1'
 
356
          stop
 
357
        elseif(abs(tmp).ge.1.d0)then
 
358
          tmp=sign(1.d0,tmp)
 
359
        endif
 
360
        tmp=acos(tmp)
 
361
      else
 
362
        tmp=1.d8
 
363
      endif
 
364
      getdelphi=tmp
 
365
      return
 
366
      end