← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

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