← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

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

  • Committer: olivier Mattelaer
  • Date: 2015-03-05 00:14:16 UTC
  • mfrom: (258.1.9 2.3)
  • mto: (258.8.1 2.3)
  • mto: This revision was merged to the branch mainline in revision 259.
  • Revision ID: olivier.mattelaer@uclouvain.be-20150305001416-y9mzeykfzwnl9t0j
partial merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/FixedOrderAnalysis/analysis_HwU_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
 
11
      character*5 cc(2)
 
12
      data cc/'     ','cuts '/
 
13
      real * 8 xmi,xms,pi
 
14
      parameter (pi=3.14159265358979312d0)
 
15
      call HwU_inithist(nwgt,weights_info)
 
16
      xmi=50.d0
 
17
      xms=130.d0
 
18
      do j=1,2
 
19
      l=(j-1)*21
 
20
      call HwU_book(l+ 1,'V pt      '//cc(j) ,100,0.d0,200.d0)
 
21
      call HwU_book(l+ 2,'V pt      '//cc(j) ,100,0.d0,1000.d0)
 
22
      call HwU_book(l+ 3,'V log[pt] '//cc(j) ,98,0.1d0,5.d0)
 
23
      call HwU_book(l+ 4,'V y       '//cc(j) ,90,-9.d0,9.d0)
 
24
      call HwU_book(l+ 5,'V eta     '//cc(j) ,90,-9.d0,9.d0)
 
25
      call HwU_book(l+ 6,'mV        '//cc(j) ,100,xmi,xms)
 
26
c
 
27
      call HwU_book(l+ 7,'lm pt      '//cc(j) ,100,0.d0,200.d0)
 
28
      call HwU_book(l+ 8,'lm pt      '//cc(j) ,100,0.d0,1000.d0)
 
29
      call HwU_book(l+ 9,'lm log[pt] '//cc(j) ,98,0.1d0,5.d0)
 
30
      call HwU_book(l+10,'lm eta     '//cc(j) ,90,-9.d0,9.d0)
 
31
      call HwU_book(l+11,'lp pt      '//cc(j) ,100,0.d0,200.d0)
 
32
      call HwU_book(l+12,'lp pt      '//cc(j) ,100,0.d0,1000.d0)
 
33
      call HwU_book(l+13,'lp log[pt] '//cc(j) ,98,0.1d0,5.d0)
 
34
      call HwU_book(l+14,'lp eta     '//cc(j) ,90,-9.d0,9.d0)
 
35
c
 
36
      call HwU_book(l+15,'lmlp delta eta     '//cc(j) ,90,-9.d0,9.d0)
 
37
      call HwU_book(l+16,'lmlp azimt         '//cc(j) ,20,0.d0,pi)
 
38
      call HwU_book(l+17,'lmlp log[pi-azimt] '//cc(j) ,82,-4.d0,0.1d0)
 
39
      call HwU_book(l+18,'lmlp inv m         '//cc(j) ,100,xmi,xms)
 
40
      call HwU_book(l+19,'lmlp pt            '//cc(j) ,100,0.d0,200.d0)
 
41
      call HwU_book(l+20,'lmlp log[pt]       '//cc(j) ,98,0.1d0,5.d0)
 
42
c
 
43
      call HwU_book(l+21,'total'//cc(j),2,-1.d0,1.d0)
 
44
      enddo
 
45
      return
 
46
      end
 
47
 
 
48
 
 
49
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
50
      subroutine analysis_end(dummy)
 
51
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
52
      implicit none
 
53
      double precision dummy
 
54
      call HwU_write_file
 
55
      return                
 
56
      end
 
57
 
 
58
      
 
59
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
60
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
 
61
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
62
      implicit none
 
63
      include 'nexternal.inc'
 
64
      integer istatus(nexternal)
 
65
      integer iPDG(nexternal)
 
66
      double precision p(0:4,nexternal)
 
67
      double precision wgts(*)
 
68
      integer ibody
 
69
      double precision wgt,var
 
70
      integer i,kk,l,nwgt_analysis
 
71
      common/c_analysis/nwgt_analysis
 
72
      double precision www,ppl(0:3),pplb(0:3),ppv(0:3),ycut,xmv,ptv,yv
 
73
     $     ,etav,ptl,yl,etal,ptlb,ylb,etalb,ptpair,azi,azinorm,xmll
 
74
     $     ,detallb
 
75
      double precision getrapidity,getinvm,getpseudorap,getdelphi
 
76
      external getrapidity,getinvm,getpseudorap,getdelphi
 
77
      double precision pi
 
78
      parameter (pi=3.14159265358979312d0)
 
79
      if (nexternal.ne.5) then
 
80
         write (*,*) 'error #1 in analysis_fill: '/
 
81
     &        /'only for process "p p > l+ l- [QCD]"'
 
82
         stop 1
 
83
      endif
 
84
      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
 
85
         write (*,*) 'error #2 in analysis_fill: '/
 
86
     &        /'only for process "p p > l+ l- [QCD]"'
 
87
         stop 1
 
88
      endif
 
89
      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
 
90
         write (*,*) 'error #3 in analysis_fill: '/
 
91
     &        /'only for process "p p > l+ l- [QCD]"'
 
92
         stop 1
 
93
      endif
 
94
      if (.not. (abs(ipdg(5)).le.5 .or. ipdg(5).eq.21)) then
 
95
         write (*,*) 'error #4 in analysis_fill: '/
 
96
     &        /'only for process "p p > l+ l- [QCD]"'
 
97
         stop 1
 
98
      endif
 
99
      if (ipdg(3).ne.-11.and.ipdg(3).ne.-13.and.ipdg(3).ne.-15) then
 
100
         write (*,*) 'error #5 in analysis_fill: '/
 
101
     &        /'only for process "p p > l+ l- [QCD]"'
 
102
         stop 1
 
103
      endif
 
104
      if (ipdg(4).ne.11.and.ipdg(4).ne.13.and.ipdg(4).ne.15) then
 
105
         write (*,*) 'error #6 in analysis_fill: '/
 
106
     &        /'only for process "p p > l+ l- [QCD]"'
 
107
         stop 1
 
108
      endif
 
109
 
 
110
      DO i=0,3
 
111
        ppl(i)=p(i,4)
 
112
        pplb(i)=p(i,3)
 
113
        ppv(i)=ppl(i)+pplb(i)
 
114
      ENDDO
 
115
 
 
116
C FILL THE HISTOS
 
117
      YCUT=2.5D0
 
118
C Variables of the vector boson
 
119
      xmv=getinvm(ppv(0),ppv(1),ppv(2),ppv(3))
 
120
      ptv=sqrt(ppv(1)**2+ppv(2)**2)
 
121
      yv=getrapidity(ppv(0),ppv(3))
 
122
      etav=getpseudorap(ppv(0),ppv(1),ppv(2),ppv(3))
 
123
C Variables of the leptons
 
124
      ptl=sqrt(ppl(1)**2+ppl(2)**2)
 
125
      yl=getrapidity(ppl(0),ppl(3))
 
126
      etal=getpseudorap(ppl(0),ppl(1),ppl(2),ppl(3))
 
127
c
 
128
      ptlb=sqrt(pplb(1)**2+pplb(2)**2)
 
129
      ylb=getrapidity(pplb(0),pplb(3))
 
130
      etalb=getpseudorap(pplb(0),pplb(1),pplb(2),pplb(3))
 
131
c
 
132
      ptpair=ptv
 
133
      azi=getdelphi(ppl(1),pplb(1),ppl(2),pplb(2))
 
134
      azinorm=(pi-azi)/pi
 
135
      xmll=xmv
 
136
      detallb=etal-etalb
 
137
c
 
138
      l=0
 
139
      call HwU_fill(l+1,(ptv),(WGTS))
 
140
      call HwU_fill(l+2,(ptv),(WGTS))
 
141
      if(ptv.gt.0.d0)call HwU_fill(l+3,(log10(ptv)),(WGTS))
 
142
      call HwU_fill(l+4,(yv),(WGTS))
 
143
      call HwU_fill(l+5,(etav),(WGTS))
 
144
      call HwU_fill(l+6,(xmv),(WGTS))
 
145
c
 
146
      call HwU_fill(l+7,(ptl),(WGTS))
 
147
      call HwU_fill(l+8,(ptl),(WGTS))
 
148
      if(ptl.gt.0.d0)call HwU_fill(l+9,(log10(ptl)),(WGTS))
 
149
      call HwU_fill(l+10,(etal),(WGTS))
 
150
      call HwU_fill(l+11,(ptlb),(WGTS))
 
151
      call HwU_fill(l+12,(ptlb),(WGTS))
 
152
      if(ptlb.gt.0.d0)call HwU_fill(l+13,(log10(ptlb)),(WGTS))
 
153
      call HwU_fill(l+14,(etalb),(WGTS))
 
154
c
 
155
      call HwU_fill(l+15,(detallb),(WGTS))
 
156
      call HwU_fill(l+16,(azi),(WGTS))
 
157
      if(azinorm.gt.0.d0) call HwU_fill(l+17,(log10(azinorm)),(WGTS))
 
158
      call HwU_fill(l+18,(xmll),(WGTS))
 
159
      call HwU_fill(l+19,(ptpair),(WGTS))
 
160
      if(ptpair.gt.0)call HwU_fill(l+20,(log10(ptpair)),(WGTS))
 
161
      call HwU_fill(l+21,(0d0),(WGTS))
 
162
c
 
163
      l=l+21
 
164
 
 
165
      if(abs(etav).lt.ycut)then
 
166
        call HwU_fill(l+1,(ptv),(WGTS))
 
167
        call HwU_fill(l+2,(ptv),(WGTS))
 
168
        if(ptv.gt.0.d0)call HwU_fill(l+3,(log10(ptv)),(WGTS))
 
169
      endif
 
170
      if(ptv.gt.20.d0)then
 
171
        call HwU_fill(l+4,(yv),(WGTS))
 
172
        call HwU_fill(l+5,(etav),(WGTS))
 
173
      endif
 
174
      if(abs(etav).lt.ycut.and.ptv.gt.20.d0)then
 
175
         call HwU_fill(l+6,(xmv),(WGTS))
 
176
         call HwU_fill(l+21,(0d0),(WGTS))
 
177
      endif
 
178
c
 
179
      if(abs(etal).lt.ycut)then
 
180
        call HwU_fill(l+7,(ptl),(WGTS))
 
181
        call HwU_fill(l+8,(ptl),(WGTS))
 
182
        if(ptl.gt.0.d0)call HwU_fill(l+9,(log10(ptl)),(WGTS))
 
183
      endif
 
184
      if(ptl.gt.20.d0)call HwU_fill(l+10,(etal),(WGTS))
 
185
      if(abs(etalb).lt.ycut)then
 
186
        call HwU_fill(l+11,(ptlb),(WGTS))
 
187
        call HwU_fill(l+12,(ptlb),(WGTS))
 
188
        if(ptlb.gt.0.d0)call HwU_fill(l+13,(log10(ptlb)),(WGTS))
 
189
      endif
 
190
      if(ptlb.gt.20.d0)call HwU_fill(l+14,(etalb),(WGTS))
 
191
c
 
192
      if( abs(etal).lt.ycut.and.abs(etalb).lt.ycut .and.
 
193
     &     ptl.gt.20.d0.and.ptlb.gt.20.d0)then
 
194
        call HwU_fill(l+15,(detallb),(WGTS))
 
195
        call HwU_fill(l+16,(azi),(WGTS))
 
196
        if(azinorm.gt.0.d0) call HwU_fill(l+17,(log10(azinorm)),(WGTS))
 
197
        call HwU_fill(l+18,(xmll),(WGTS))
 
198
        call HwU_fill(l+19,(ptpair),(WGTS))
 
199
        if(ptpair.gt.0) call HwU_fill(l+20,(log10(ptpair)),(WGTS))
 
200
      endif
 
201
 
 
202
 999  return      
 
203
      end
 
204
 
 
205
 
 
206
      function getrapidity(en,pl)
 
207
      implicit none
 
208
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
209
      parameter (tiny=1.d-8)
 
210
      xplus=en+pl
 
211
      xminus=en-pl
 
212
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
213
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
214
            y=0.5d0*log( xplus/xminus  )
 
215
         else
 
216
            y=sign(1.d0,pl)*1.d8
 
217
         endif
 
218
      else 
 
219
         y=sign(1.d0,pl)*1.d8
 
220
      endif
 
221
      getrapidity=y
 
222
      return
 
223
      end
 
224
 
 
225
 
 
226
      function getpseudorap(en,ptx,pty,pl)
 
227
      implicit none
 
228
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
229
      parameter (tiny=1.d-5)
 
230
c
 
231
      pt=sqrt(ptx**2+pty**2)
 
232
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
233
        eta=sign(1.d0,pl)*1.d8
 
234
      else
 
235
        th=atan2(pt,pl)
 
236
        eta=-log(tan(th/2.d0))
 
237
      endif
 
238
      getpseudorap=eta
 
239
      return
 
240
      end
 
241
 
 
242
 
 
243
      function getinvm(en,ptx,pty,pl)
 
244
      implicit none
 
245
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
246
      parameter (tiny=1.d-5)
 
247
c
 
248
      tmp=en**2-ptx**2-pty**2-pl**2
 
249
      if(tmp.gt.0.d0)then
 
250
        tmp=sqrt(tmp)
 
251
      elseif(tmp.gt.-tiny)then
 
252
        tmp=0.d0
 
253
      else
 
254
        write(*,*)'Attempt to compute a negative mass'
 
255
        stop
 
256
      endif
 
257
      getinvm=tmp
 
258
      return
 
259
      end
 
260
 
 
261
 
 
262
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
263
      implicit none
 
264
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
265
      parameter (tiny=1.d-5)
 
266
c
 
267
      pt1=sqrt(ptx1**2+pty1**2)
 
268
      pt2=sqrt(ptx2**2+pty2**2)
 
269
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
270
        tmp=ptx1*ptx2+pty1*pty2
 
271
        tmp=tmp/(pt1*pt2)
 
272
        if(abs(tmp).gt.1.d0+tiny)then
 
273
          write(*,*)'Cosine larger than 1'
 
274
          stop
 
275
        elseif(abs(tmp).ge.1.d0)then
 
276
          tmp=sign(1.d0,tmp)
 
277
        endif
 
278
        tmp=acos(tmp)
 
279
      else
 
280
        tmp=1.d8
 
281
      endif
 
282
      getdelphi=tmp
 
283
      return
 
284
      end