| 1 | c
|
|---|
| 2 | c
|
|---|
| 3 | c Plotting routines
|
|---|
| 4 | c
|
|---|
| 5 | c
|
|---|
| 6 | subroutine initplot
|
|---|
| 7 | implicit none
|
|---|
| 8 | c Book histograms in this routine. Use mbook or bookup. The entries
|
|---|
| 9 | c of these routines are real*8
|
|---|
| 10 | double precision emax,ebin,etamin,etamax,etabin,xmi,xms,bin
|
|---|
| 11 | double precision pi
|
|---|
| 12 | parameter (pi=3.1415926535897932385d0)
|
|---|
| 13 | integer i,kk
|
|---|
| 14 | include 'run.inc'
|
|---|
| 15 | character*5 cc(2)
|
|---|
| 16 | data cc/' ',' cuts'/
|
|---|
| 17 | c
|
|---|
| 18 | emax=dsqrt(ebeam(1)*ebeam(2))
|
|---|
| 19 | ebin=emax/50.d0
|
|---|
| 20 | etamin=-5.d0
|
|---|
| 21 | etamax=5.d0
|
|---|
| 22 | etabin=0.2d0
|
|---|
| 23 | c resets histograms
|
|---|
| 24 | call inihist
|
|---|
| 25 | c
|
|---|
| 26 | xmi=50.d0
|
|---|
| 27 | xms=130.d0
|
|---|
| 28 | bin=0.8d0
|
|---|
| 29 |
|
|---|
| 30 | do i=1,2
|
|---|
| 31 | kk=(i-1)*50
|
|---|
| 32 | call bookup(kk+ 1,'V pt'//cc(i),2.d0,0.d0,200.d0)
|
|---|
| 33 | call bookup(kk+ 2,'V pt'//cc(i),10.d0,0.d0,1000.d0)
|
|---|
| 34 | call bookup(kk+ 3,'V log[pt]'//cc(i),0.05d0,0.1d0,5.d0)
|
|---|
| 35 | call bookup(kk+ 4,'V y'//cc(i),0.2d0,-9.d0,9.d0)
|
|---|
| 36 | call bookup(kk+ 5,'V eta'//cc(i),0.2d0,-9.d0,9.d0)
|
|---|
| 37 | call bookup(kk+ 6,'mV'//cc(i),bin,xmi,xms)
|
|---|
| 38 | c
|
|---|
| 39 | call bookup(kk+ 7,'l pt'//cc(i),2.d0,0.d0,200.d0)
|
|---|
| 40 | call bookup(kk+ 8,'l pt'//cc(i),10.d0,0.d0,1000.d0)
|
|---|
| 41 | call bookup(kk+ 9,'l log[pt]'//cc(i),0.05d0,0.1d0,5.d0)
|
|---|
| 42 | call bookup(kk+10,'l eta'//cc(i),0.2d0,-9.d0,9.d0)
|
|---|
| 43 | call bookup(kk+11,'lb pt'//cc(i),2.d0,0.d0,200.d0)
|
|---|
| 44 | call bookup(kk+12,'lb pt'//cc(i),10.d0,0.d0,1000.d0)
|
|---|
| 45 | call bookup(kk+13,'lb log[pt]'//cc(i),0.05d0,0.1d0,5.d0)
|
|---|
| 46 | call bookup(kk+14,'lb eta'//cc(i),0.2d0,-9.d0,9.d0)
|
|---|
| 47 | c
|
|---|
| 48 | call bookup(kk+15,'llb delta eta'//cc(i),0.2d0,-9.d0,9.d0)
|
|---|
| 49 | call bookup(kk+16,'llb azimt'//cc(i),pi/20.d0,0.d0,pi)
|
|---|
| 50 | call bookup(kk+17,'llb log[pi-azimt]'//cc(i),0.05d0,-4.d0,0.1d0)
|
|---|
| 51 | call bookup(kk+18,'llb inv m'//cc(i),bin,xmi,xms)
|
|---|
| 52 | call bookup(kk+19,'llb pt'//cc(i),2.d0,0.d0,200.d0)
|
|---|
| 53 | call bookup(kk+20,'llb log[pt]'//cc(i),0.05d0,0.1d0,5.d0)
|
|---|
| 54 | c
|
|---|
| 55 | call bookup(kk+21,'total'//cc(i),1.d0,-1.d0,1.d0)
|
|---|
| 56 | enddo
|
|---|
| 57 | return
|
|---|
| 58 | end
|
|---|
| 59 |
|
|---|
| 60 |
|
|---|
| 61 | subroutine topout
|
|---|
| 62 | implicit none
|
|---|
| 63 | character*14 ytit
|
|---|
| 64 | logical usexinteg,mint
|
|---|
| 65 | common/cusexinteg/usexinteg,mint
|
|---|
| 66 | integer itmax,ncall
|
|---|
| 67 | common/citmax/itmax,ncall
|
|---|
| 68 | real*8 xnorm1,xnorm2
|
|---|
| 69 | logical unwgt
|
|---|
| 70 | double precision evtsgn
|
|---|
| 71 | common /c_unwgt/evtsgn,unwgt
|
|---|
| 72 | integer i,kk
|
|---|
| 73 | include 'dbook.inc'
|
|---|
| 74 | c
|
|---|
| 75 | if (unwgt) then
|
|---|
| 76 | ytit='events per bin'
|
|---|
| 77 | else
|
|---|
| 78 | ytit='sigma per bin '
|
|---|
| 79 | endif
|
|---|
| 80 | xnorm1=1.d0/float(itmax)
|
|---|
| 81 | xnorm2=1.d0/float(ncall*itmax)
|
|---|
| 82 | do i=1,NPLOTS
|
|---|
| 83 | if(usexinteg.and..not.mint) then
|
|---|
| 84 | call mopera(i,'+',i,i,xnorm1,0.d0)
|
|---|
| 85 | elseif(mint) then
|
|---|
| 86 | call mopera(i,'+',i,i,xnorm2,0.d0)
|
|---|
| 87 | endif
|
|---|
| 88 | call mfinal(i)
|
|---|
| 89 | enddo
|
|---|
| 90 | do i=1,2
|
|---|
| 91 | kk=(i-1)*50
|
|---|
| 92 | call multitop(kk+ 1,3,2,'V pt',' ','LOG')
|
|---|
| 93 | call multitop(kk+ 2,3,2,'V pt',' ','LOG')
|
|---|
| 94 | call multitop(kk+ 3,3,2,'V log[pt]',' ','LOG')
|
|---|
| 95 | call multitop(kk+ 4,3,2,'V y',' ','LOG')
|
|---|
| 96 | call multitop(kk+ 5,3,2,'V eta',' ','LOG')
|
|---|
| 97 | call multitop(kk+ 6,3,2,'mV',' ','LOG')
|
|---|
| 98 | enddo
|
|---|
| 99 | do i=1,2
|
|---|
| 100 | kk=(i-1)*50
|
|---|
| 101 | call multitop(kk+ 7,3,2,'l pt',' ','LOG')
|
|---|
| 102 | call multitop(kk+ 8,3,2,'l pt',' ','LOG')
|
|---|
| 103 | call multitop(kk+ 9,3,2,'l log[pt]',' ','LOG')
|
|---|
| 104 | call multitop(kk+10,3,2,'l eta',' ','LOG')
|
|---|
| 105 | call multitop(kk+11,3,2,'l pt',' ','LOG')
|
|---|
| 106 | call multitop(kk+12,3,2,'l pt',' ','LOG')
|
|---|
| 107 | call multitop(kk+13,3,2,'l log[pt]',' ','LOG')
|
|---|
| 108 | call multitop(kk+14,3,2,'l eta',' ','LOG')
|
|---|
| 109 | c
|
|---|
| 110 | call multitop(kk+15,3,2,'llb deta',' ','LOG')
|
|---|
| 111 | call multitop(kk+16,3,2,'llb azi',' ','LOG')
|
|---|
| 112 | call multitop(kk+17,3,2,'llb azi',' ','LOG')
|
|---|
| 113 | call multitop(kk+18,3,2,'llb inv m',' ','LOG')
|
|---|
| 114 | call multitop(kk+19,3,2,'llb pt',' ','LOG')
|
|---|
| 115 | call multitop(kk+20,3,2,'llb pt',' ','LOG')
|
|---|
| 116 | c
|
|---|
| 117 | call multitop(kk+21,3,2,'total',' ','LOG')
|
|---|
| 118 | enddo
|
|---|
| 119 | return
|
|---|
| 120 | end
|
|---|
| 121 |
|
|---|
| 122 |
|
|---|
| 123 | subroutine outfun(pp,ybst_til_tolab,www,itype)
|
|---|
| 124 | C
|
|---|
| 125 | C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
|
|---|
| 126 | C
|
|---|
| 127 | C In MadFKS, the momenta PP given in input to this function are in the
|
|---|
| 128 | C reduced parton c.m. frame. If need be, boost them to the lab frame.
|
|---|
| 129 | C The rapidity of this boost is
|
|---|
| 130 | C
|
|---|
| 131 | C YBST_TIL_TOLAB
|
|---|
| 132 | C
|
|---|
| 133 | C also given in input
|
|---|
| 134 | C
|
|---|
| 135 | C This is the rapidity that enters in the arguments of the sinh() and
|
|---|
| 136 | C cosh() of the boost, in such a way that
|
|---|
| 137 | C ylab = ycm - ybst_til_tolab
|
|---|
| 138 | C where ylab is the rapidity in the lab frame and ycm the rapidity
|
|---|
| 139 | C in the center-of-momentum frame.
|
|---|
| 140 | C
|
|---|
| 141 | C *WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING**WARNING*
|
|---|
| 142 | implicit none
|
|---|
| 143 | include 'nexternal.inc'
|
|---|
| 144 | real*8 pp(0:3,nexternal),ybst_til_tolab,www
|
|---|
| 145 | integer itype
|
|---|
| 146 | real*8 var
|
|---|
| 147 | real*8 ppevs(0:3,nexternal)
|
|---|
| 148 |
|
|---|
| 149 | real*8 getrapidity,getpseudorap,chybst,shybst,chybstmo
|
|---|
| 150 | real*8 xd(1:3)
|
|---|
| 151 | data (xd(i),i=1,3)/0,0,1/
|
|---|
| 152 | real*8 pplab(0:3,nexternal)
|
|---|
| 153 | double precision ppcl(4,nexternal),y(nexternal)
|
|---|
| 154 | double precision pjet(0:3,nexternal)
|
|---|
| 155 | double precision cthjet(nexternal)
|
|---|
| 156 | integer nn,njet,nsub,jet(nexternal)
|
|---|
| 157 | real*8 emax,getcth,cpar,dpar,thrust,dot,shat
|
|---|
| 158 | integer i,j,kk,imax
|
|---|
| 159 |
|
|---|
| 160 | LOGICAL IS_A_J(NEXTERNAL),IS_A_LP(NEXTERNAL),IS_A_LM(NEXTERNAL)
|
|---|
| 161 | COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM
|
|---|
| 162 | c masses
|
|---|
| 163 | double precision pmass(nexternal)
|
|---|
| 164 | double precision pt1, eta1, y1, pt2, eta2, y2, pt3, eta3, y3, ht
|
|---|
| 165 |
|
|---|
| 166 | common/to_mass/pmass
|
|---|
| 167 | double precision ppkt(0:3,nexternal),ecut,djet
|
|---|
| 168 | double precision ycut, palg
|
|---|
| 169 | integer ipartjet(nexternal)
|
|---|
| 170 |
|
|---|
| 171 | double precision pi
|
|---|
| 172 | parameter (pi=3.1415926535897932385d0)
|
|---|
| 173 | double precision tiny
|
|---|
| 174 | parameter (tiny=1d-5)
|
|---|
| 175 | double precision ppl(5),pplb(5),ppv(5),xmv,ptv,yv,thv,etav,ptl,yl
|
|---|
| 176 | $ ,thl,etal,pll,enl,ptlb,ylb,thlb,etalb,pllb,enlb,ptpair,dll
|
|---|
| 177 | $ ,cll,azi,azinorm,xmll,detallb,wmass,wgamma,bwcutoff
|
|---|
| 178 | logical flag
|
|---|
| 179 |
|
|---|
| 180 | c
|
|---|
| 181 | if(itype.eq.11.or.itype.eq.12)then
|
|---|
| 182 | kk=0
|
|---|
| 183 | elseif(itype.eq.20)then
|
|---|
| 184 | kk=50
|
|---|
| 185 | return
|
|---|
| 186 | else
|
|---|
| 187 | write(*,*)'Error in outfun: unknown itype',itype
|
|---|
| 188 | stop
|
|---|
| 189 | endif
|
|---|
| 190 |
|
|---|
| 191 |
|
|---|
| 192 |
|
|---|
| 193 | chybst=cosh(ybst_til_tolab)
|
|---|
| 194 | shybst=sinh(ybst_til_tolab)
|
|---|
| 195 | chybstmo=chybst-1.d0
|
|---|
| 196 | do i=3,nexternal
|
|---|
| 197 | call boostwdir2(chybst,shybst,chybstmo,xd,
|
|---|
| 198 | # pp(0,i),pplab(0,i))
|
|---|
| 199 | enddo
|
|---|
| 200 |
|
|---|
| 201 |
|
|---|
| 202 | DO I=1,4
|
|---|
| 203 | PPL(I)=pplab(mod(I,4),4)
|
|---|
| 204 | PPLB(I)=pplab(mod(I,4),3)
|
|---|
| 205 | ENDDO
|
|---|
| 206 | ppl(5)=0d0
|
|---|
| 207 | pplb(5)=0d0
|
|---|
| 208 |
|
|---|
| 209 | do i=1,4
|
|---|
| 210 | ppv(i)=ppl(i)+pplb(i)
|
|---|
| 211 | enddo
|
|---|
| 212 | ppv(5)=sqrt(ppv(4)**2-ppv(1)**2-ppv(2)**2-ppv(3)**2)
|
|---|
| 213 |
|
|---|
| 214 | C FILL THE HISTOS
|
|---|
| 215 | YCUT=2.5D0
|
|---|
| 216 | C Variables of the vector boson
|
|---|
| 217 | xmv=ppv(5)
|
|---|
| 218 | ptv=sqrt(ppv(1)**2+ppv(2)**2)
|
|---|
| 219 | if(abs(ppv(4)-abs(ppv(3))).gt.tiny)then
|
|---|
| 220 | yv=0.5d0*log( (ppv(4)+ppv(3))/
|
|---|
| 221 | # (ppv(4)-ppv(3)) )
|
|---|
| 222 | else
|
|---|
| 223 | yv=sign(1.d0,ppv(3))*1.d8
|
|---|
| 224 | endif
|
|---|
| 225 | thv = atan2(ptv+tiny,ppv(3))
|
|---|
| 226 | etav= -log(tan(thv/2))
|
|---|
| 227 | C Variables of the leptons
|
|---|
| 228 | ptl=sqrt(ppl(1)**2+ppl(2)**2)
|
|---|
| 229 | if(abs(ppl(4)-abs(ppl(3))).gt.tiny)then
|
|---|
| 230 | yl=0.5d0*log( (ppl(4)+ppl(3))/
|
|---|
| 231 | # (ppl(4)-ppl(3)) )
|
|---|
| 232 | else
|
|---|
| 233 | yl=sign(1.d0,ppl(3))*1.d8
|
|---|
| 234 | endif
|
|---|
| 235 | thl = atan2(ptl+tiny,ppl(3))
|
|---|
| 236 | etal= -log(tan(thl/2))
|
|---|
| 237 | pll = ppl(3)
|
|---|
| 238 | enl = ppl(4)
|
|---|
| 239 | c
|
|---|
| 240 | ptlb=sqrt(pplb(1)**2+pplb(2)**2)
|
|---|
| 241 | if(abs(pplb(4)-abs(pplb(3))).gt.tiny)then
|
|---|
| 242 | ylb=0.5d0*log( (pplb(4)+pplb(3))/
|
|---|
| 243 | # (pplb(4)-pplb(3)) )
|
|---|
| 244 | else
|
|---|
| 245 | ylb=sign(1.d0,pplb(3))*1.d8
|
|---|
| 246 | endif
|
|---|
| 247 | thlb = atan2(ptlb+tiny,pplb(3))
|
|---|
| 248 | etalb= -log(tan(thlb/2))
|
|---|
| 249 | pllb = pplb(3)
|
|---|
| 250 | enlb = pplb(4)
|
|---|
| 251 | c
|
|---|
| 252 | ptpair = dsqrt((ppl(1)+pplb(1))**2+(ppl(2)+pplb(2))**2)
|
|---|
| 253 | dll = ppl(1)*pplb(1)+ppl(2)*pplb(2)
|
|---|
| 254 | cll=0
|
|---|
| 255 | if(ptl.ne.0.and.ptlb.ne.0) cll = dll * (1-tiny)/(ptl*ptlb)
|
|---|
| 256 | if(abs(cll).gt.1) then
|
|---|
| 257 | write(*,*) ' cosine = ',cll ,dll,ptl,ptlb
|
|---|
| 258 | cll = - 1
|
|---|
| 259 | endif
|
|---|
| 260 | azi = (1-tiny)*acos(cll)
|
|---|
| 261 | azinorm = (pi-azi)/pi
|
|---|
| 262 | xmll = dsqrt( ppl(5)**2 + ppl(5)**2 +
|
|---|
| 263 | # 2*(enl*enlb - pll*pllb - dll) )
|
|---|
| 264 | detallb = etal-etalb
|
|---|
| 265 | c
|
|---|
| 266 | kk=0
|
|---|
| 267 | wmass=80.419d0
|
|---|
| 268 | wgamma=2.046d0
|
|---|
| 269 | bwcutoff=15.d0
|
|---|
| 270 |
|
|---|
| 271 | flag=(xmv.ge.wmass-wgamma*bwcutoff.and.
|
|---|
| 272 | & xmv.le.wmass+wgamma*bwcutoff)
|
|---|
| 273 |
|
|---|
| 274 | if(flag)then
|
|---|
| 275 | call mfill(kk+1,ptv,WWW)
|
|---|
| 276 | call mfill(kk+2,ptv,WWW)
|
|---|
| 277 | if(ptv.gt.0.d0)call mfill(kk+3,log10(ptv),WWW)
|
|---|
| 278 | call mfill(kk+4,yv,WWW)
|
|---|
| 279 | call mfill(kk+5,etav,WWW)
|
|---|
| 280 | call mfill(kk+6,xmv,WWW)
|
|---|
| 281 | c
|
|---|
| 282 | call mfill(kk+7,ptl,WWW)
|
|---|
| 283 | call mfill(kk+8,ptl,WWW)
|
|---|
| 284 | if(ptl.gt.0.d0)call mfill(kk+9,log10(ptl),WWW)
|
|---|
| 285 | call mfill(kk+10,etal,WWW)
|
|---|
| 286 | call mfill(kk+11,ptlb,WWW)
|
|---|
| 287 | call mfill(kk+12,ptlb,WWW)
|
|---|
| 288 | if(ptlb.gt.0.d0)call mfill(kk+13,log10(ptlb),WWW)
|
|---|
| 289 | call mfill(kk+14,etalb,WWW)
|
|---|
| 290 | c
|
|---|
| 291 | call mfill(kk+15,detallb,WWW)
|
|---|
| 292 | call mfill(kk+16,azi,WWW)
|
|---|
| 293 | if(azinorm.gt.0.d0)
|
|---|
| 294 | # call mfill(kk+17,log10(azinorm),WWW)
|
|---|
| 295 | call mfill(kk+18,xmll,WWW)
|
|---|
| 296 | call mfill(kk+19,ptpair,WWW)
|
|---|
| 297 | if(ptpair.gt.0)call mfill(kk+20,log10(ptpair),WWW)
|
|---|
| 298 | call mfill(kk+21,0d0,WWW)
|
|---|
| 299 | c
|
|---|
| 300 | kk=50
|
|---|
| 301 | if(abs(etav).lt.ycut)then
|
|---|
| 302 | call mfill(kk+1,ptv,WWW)
|
|---|
| 303 | call mfill(kk+2,ptv,WWW)
|
|---|
| 304 | if(ptv.gt.0.d0)call mfill(kk+3,log10(ptv),WWW)
|
|---|
| 305 | endif
|
|---|
| 306 | if(ptv.gt.20.d0)then
|
|---|
| 307 | call mfill(kk+4,yv,WWW)
|
|---|
| 308 | call mfill(kk+5,etav,WWW)
|
|---|
| 309 | endif
|
|---|
| 310 | if(abs(etav).lt.ycut.and.ptv.gt.20.d0)then
|
|---|
| 311 | call mfill(kk+6,xmv,WWW)
|
|---|
| 312 | call mfill(kk+21,0d0,WWW)
|
|---|
| 313 | endif
|
|---|
| 314 | c
|
|---|
| 315 | if(abs(etal).lt.ycut)then
|
|---|
| 316 | call mfill(kk+7,ptl,WWW)
|
|---|
| 317 | call mfill(kk+8,ptl,WWW)
|
|---|
| 318 | if(ptl.gt.0.d0)call mfill(kk+9,log10(ptl),WWW)
|
|---|
| 319 | endif
|
|---|
| 320 | if(ptl.gt.20.d0)call mfill(kk+10,etal,WWW)
|
|---|
| 321 | if(abs(etalb).lt.ycut)then
|
|---|
| 322 | call mfill(kk+11,ptlb,WWW)
|
|---|
| 323 | call mfill(kk+12,ptlb,WWW)
|
|---|
| 324 | if(ptlb.gt.0.d0)call mfill(kk+13,log10(ptlb),WWW)
|
|---|
| 325 | endif
|
|---|
| 326 | if(ptlb.gt.20.d0)call mfill(kk+14,etalb,WWW)
|
|---|
| 327 | c
|
|---|
| 328 | if( abs(etal).lt.ycut.and.abs(etalb).lt.ycut .and.
|
|---|
| 329 | # ptl.gt.20.d0.and.ptlb.gt.20.d0)then
|
|---|
| 330 | call mfill(kk+15,detallb,WWW)
|
|---|
| 331 | call mfill(kk+16,azi,WWW)
|
|---|
| 332 | if(azinorm.gt.0.d0)
|
|---|
| 333 | # call mfill(kk+17,log10(azinorm),WWW)
|
|---|
| 334 | call mfill(kk+18,xmll,WWW)
|
|---|
| 335 | call mfill(kk+19,ptpair,WWW)
|
|---|
| 336 | if(ptpair.gt.0)
|
|---|
| 337 | # call mfill(kk+20,log10(ptpair),WWW)
|
|---|
| 338 | endif
|
|---|
| 339 | endif
|
|---|
| 340 |
|
|---|
| 341 | 999 return
|
|---|
| 342 | end
|
|---|
| 343 |
|
|---|
| 344 |
|
|---|
| 345 | function getrapidity(en,px,py,pl)
|
|---|
| 346 | implicit none
|
|---|
| 347 | real*8 getrapidity,en,px,py,pl,tiny,xplus,xminus,y
|
|---|
| 348 | parameter (tiny=1.d-8)
|
|---|
| 349 | if (abs(en).lt.abs(pl)) en = dsqrt(px**2+py**2+pl**2)
|
|---|
| 350 | c
|
|---|
| 351 | xplus=en+pl
|
|---|
| 352 | xminus=en-pl
|
|---|
| 353 | if(xplus.gt.tiny.and.xminus.gt.tiny)then
|
|---|
| 354 | if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
|
|---|
| 355 | y=0.5d0*log( xplus/xminus )
|
|---|
| 356 | else
|
|---|
| 357 | y=sign(1.d0,pl)*1.d8
|
|---|
| 358 | endif
|
|---|
| 359 | else
|
|---|
| 360 | y=sign(1.d0,pl)*1.d8
|
|---|
| 361 | endif
|
|---|
| 362 | getrapidity=y
|
|---|
| 363 | return
|
|---|
| 364 | end
|
|---|
| 365 |
|
|---|
| 366 |
|
|---|
| 367 | function getpseudorap(en,ptx,pty,pl)
|
|---|
| 368 | implicit none
|
|---|
| 369 | real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
|
|---|
| 370 | parameter (tiny=1.d-5)
|
|---|
| 371 | c
|
|---|
| 372 | pt=sqrt(ptx**2+pty**2)
|
|---|
| 373 | if(pt.lt.tiny.and.abs(pl).lt.tiny)then
|
|---|
| 374 | eta=sign(1.d0,pl)*1.d8
|
|---|
| 375 | else
|
|---|
| 376 | th=atan2(pt,pl)
|
|---|
| 377 | eta=-log(tan(th/2.d0))
|
|---|
| 378 | endif
|
|---|
| 379 | getpseudorap=eta
|
|---|
| 380 | return
|
|---|
| 381 | end
|
|---|
| 382 |
|
|---|
