← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2

« back to all changes in this revision

Viewing changes to Template/NLO/FixedOrderAnalysis/analysis_td_pp_V.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/FixedOrderAnalysis/analysis_td_pp_V.f
 
1
c
 
2
c Example analysis for "p p > w+ [QCD]" process.
 
3
c Example analysis for "p p > w- [QCD]" process.
 
4
c Example analysis for "p p > z [QCD]" process.
 
5
c
 
6
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
7
      subroutine analysis_begin(nwgt,weights_info)
 
8
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
9
      implicit none
 
10
      integer nwgt
 
11
      character*(*) weights_info(*)
 
12
      integer j,kk,l,nwgt_analysis
 
13
      common/c_analysis/nwgt_analysis
 
14
      character*5 cc(2)
 
15
      data cc/'     ',' Born'/
 
16
      real * 8 bin,xmi,xms,pi
 
17
      parameter (pi=3.14159265358979312d0)
 
18
      include 'dbook.inc'
 
19
      call inihist
 
20
      nwgt_analysis=nwgt
 
21
      if (nwgt_analysis*10.gt.nplots/4) then
 
22
         write (*,*) 'error in analysis_begin: '/
 
23
     &        /'too many histograms, increase NPLOTS to',
 
24
     &        nwgt_analysis*10*4
 
25
         stop 1
 
26
      endif
 
27
      xmi=40.d0
 
28
      xms=120.d0
 
29
      bin=1.0d0
 
30
      do j=1,2
 
31
      do kk=1,nwgt_analysis
 
32
      l=(kk-1)*10+(j-1)*5
 
33
      call bookup(l+ 1,'V pt     '//weights_info(kk)//cc(j)
 
34
     &     ,2.d0,0.d0,200.d0)
 
35
      call bookup(l+ 2,'V log pt '//weights_info(kk)//cc(j)
 
36
     &     ,0.05d0,0.d0,5.d0)
 
37
      call bookup(l+ 3,'V y      '//weights_info(kk)//cc(j)
 
38
     &     ,0.25d0,-9.d0,9.d0)
 
39
      call bookup(l+ 4,'V eta    '//weights_info(kk)//cc(j)
 
40
     &     ,0.25d0,-9.d0,9.d0)
 
41
      call bookup(l+ 5,'mV       '//weights_info(kk)//cc(j)
 
42
     &     ,bin,xmi,xms)
 
43
      enddo
 
44
      enddo
 
45
      return
 
46
      end
 
47
 
 
48
 
 
49
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
50
      subroutine analysis_end(xnorm)
 
51
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
52
      implicit none
 
53
      character*14 ytit
 
54
      double precision xnorm
 
55
      integer i
 
56
      integer kk,l,nwgt_analysis
 
57
      common/c_analysis/nwgt_analysis
 
58
      include 'dbook.inc'
 
59
      call open_topdrawer_file
 
60
      call mclear
 
61
      do i=1,NPLOTS
 
62
         call mopera(i,'+',i,i,xnorm,0.d0)
 
63
         call mfinal(i)
 
64
      enddo
 
65
      ytit='sigma per bin '
 
66
      do i=1,2
 
67
      do kk=1,nwgt_analysis
 
68
      l=(kk-1)*10+(i-1)*5
 
69
      call multitop(l+ 1,3,2,'V pt     ',' ','LOG')
 
70
      call multitop(l+ 2,3,2,'V log pt ',' ','LOG')
 
71
      call multitop(l+ 3,3,2,'V y      ',' ','LOG')
 
72
      call multitop(l+ 4,3,2,'V eta    ',' ','LOG')
 
73
      call multitop(l+ 5,3,2,'mV       ',' ','LOG')
 
74
      enddo
 
75
      enddo
 
76
      call close_topdrawer_file
 
77
      return                
 
78
      end
 
79
 
 
80
 
 
81
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
82
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
 
83
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
84
      implicit none
 
85
      include 'nexternal.inc'
 
86
      integer istatus(nexternal)
 
87
      integer iPDG(nexternal)
 
88
      double precision p(0:4,nexternal)
 
89
      double precision wgts(*)
 
90
      integer ibody
 
91
      double precision wgt,var
 
92
      integer i,kk,l,nwgt_analysis
 
93
      common/c_analysis/nwgt_analysis
 
94
      double precision www,pv(0:3),xmv,ptv,yv,etav
 
95
      double precision getrapidity,getpseudorap,dot
 
96
      external getrapidity,getpseudorap,dot
 
97
      if (nexternal.ne.4) then
 
98
         write (*,*) 'error #1 in analysis_fill: '/
 
99
     &        /'only for process "p p > V [QCD]"'
 
100
         stop 1
 
101
      endif
 
102
      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
 
103
         write (*,*) 'error #2 in analysis_fill: '/
 
104
     &        /'only for process "p p > V [QCD]"'
 
105
         stop 1
 
106
      endif
 
107
      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
 
108
         write (*,*) 'error #3 in analysis_fill: '/
 
109
     &        /'only for process "p p > V [QCD]"'
 
110
         stop 1
 
111
      endif
 
112
      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
 
113
         write (*,*) 'error #4 in analysis_fill: '/
 
114
     &        /'only for process "p p > V [QCD]"'
 
115
         stop 1
 
116
      endif
 
117
      if (abs(ipdg(3)).ne.24.and.ipdg(3).ne.23) then
 
118
         write (*,*) 'error #5 in analysis_fill: '/
 
119
     &        /'only for process "p p > V [QCD]"'
 
120
         stop 1
 
121
      endif
 
122
C
 
123
      do i=0,3
 
124
        pv(i)=p(i,3)
 
125
      enddo
 
126
      xmv=sqrt(max(dot(pv,pv),0d0))
 
127
      ptv=sqrt(max(pv(1)**2+pv(2)**2,0d0))
 
128
      yv=getrapidity(pv(0),pv(3))
 
129
      etav=getpseudorap(pv(0),pv(1),pv(2),pv(3))
 
130
C
 
131
      do i=1,2
 
132
         do kk=1,nwgt_analysis
 
133
            www=wgts(kk)
 
134
            l=(kk-1)*10+(i-1)*5
 
135
            if (ibody.ne.3 .and.i.eq.2) cycle
 
136
            call mfill(l+1,ptv,WWW)
 
137
            if(ptv.gt.0) call mfill(l+2,log10(ptv),WWW)
 
138
            call mfill(l+3,yv,WWW)
 
139
            call mfill(l+4,etav,WWW)
 
140
            call mfill(l+5,xmv,WWW)
 
141
         enddo
 
142
      enddo
 
143
C
 
144
 999  return      
 
145
      end
 
146
 
 
147
 
 
148
 
 
149
      function getrapidity(en,pl)
 
150
      implicit none
 
151
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
152
      parameter (tiny=1.d-8)
 
153
      xplus=en+pl
 
154
      xminus=en-pl
 
155
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
156
         if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny)then
 
157
            y=0.5d0*log( xplus/xminus  )
 
158
         else
 
159
            y=sign(1.d0,pl)*1.d8
 
160
         endif
 
161
      else 
 
162
         y=sign(1.d0,pl)*1.d8
 
163
      endif
 
164
      getrapidity=y
 
165
      return
 
166
      end
 
167
 
 
168
 
 
169
      function getpseudorap(en,ptx,pty,pl)
 
170
      implicit none
 
171
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
172
      parameter (tiny=1.d-5)
 
173
c
 
174
      pt=sqrt(ptx**2+pty**2)
 
175
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
176
        eta=sign(1.d0,pl)*1.d8
 
177
      else
 
178
        th=atan2(pt,pl)
 
179
        eta=-log(tan(th/2.d0))
 
180
      endif
 
181
      getpseudorap=eta
 
182
      return
 
183
      end