← Back to branch summary

~madteam/mg5amcnlo/series2.0

« back to all changes in this revision

Viewing changes to Template/NLO/FixedOrderAnalysis/analysis_HwU_pp_hjj.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_hjj.f
 
1
c
 
2
c Example analysis for "p p > h j j [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*8 cc(2)
 
12
      data cc/'        ','vbfcuts '/
 
13
      double precision pi
 
14
      PARAMETER (PI=3.14159265358979312D0)
 
15
      include 'dbook.inc'
 
16
      real*8 vetomin, vetomax
 
17
      integer nbinveto
 
18
      common /to_veto_hist/vetomin,vetomax,nbinveto
 
19
c      
 
20
      call HwU_inithist(nwgt,weights_info)
 
21
      vetomin = 0d0
 
22
      vetomax = 100d0
 
23
      nbinveto = 50
 
24
      do i=1,2
 
25
         l=(i-1)*54
 
26
         call HwU_book(l+  1,'total ' //cc(i),5,0.5d0,5.5d0)
 
27
 
 
28
         call HwU_book(l+  2,'Higgs pT ' //cc(i),50,0.d0,400.d0)
 
29
         call HwU_book(l+  3,'Higgs pT ' //cc(i),50,0.d0,800.d0)
 
30
         call HwU_book(l+  4,'Higgs logpT ' //cc(i),50,0.d0,4.d0)
 
31
         call HwU_book(l+  5,'Higgs eta ' //cc(i),50,-6.d0,6.d0)
 
32
         call HwU_book(l+  6,'Higgs y ' //cc(i),50,-6.d0,6.d0)
 
33
 
 
34
         call HwU_book(l+  7,'j1 pT ' //cc(i),50,0.d0,400.d0)
 
35
         call HwU_book(l+  8,'j1 pT ' //cc(i),50,0.d0,800.d0)
 
36
         call HwU_book(l+  9,'j1 logpT ' //cc(i),50,0.d0,4.d0)
 
37
         call HwU_book(l+ 10,'j1 eta ' //cc(i),50,-6.d0,6.d0)
 
38
         call HwU_book(l+ 11,'j1 y ' //cc(i),50,-6.d0,6.d0)
 
39
 
 
40
         call HwU_book(l+ 12,'j2 pT ' //cc(i),50,0.d0,400.d0)
 
41
         call HwU_book(l+ 13,'j2 pT ' //cc(i),50,0.d0,800.d0)
 
42
         call HwU_book(l+ 14,'j2 logpT ' //cc(i),50,0.d0,4.d0)
 
43
         call HwU_book(l+ 15,'j2 eta ' //cc(i),50,-6.d0,6.d0)
 
44
         call HwU_book(l+ 16,'j2 y ' //cc(i),50,-6.d0,6.d0)
 
45
 
 
46
         call HwU_book(l+ 17,'j3 pT ' //cc(i),50,0.d0,400.d0)
 
47
         call HwU_book(l+ 18,'j3 pT ' //cc(i),50,0.d0,800.d0)
 
48
         call HwU_book(l+ 19,'j3 logpT ' //cc(i),50,0.d0,4.d0)
 
49
         call HwU_book(l+ 20,'j3 eta ' //cc(i),50,-6.d0,6.d0)
 
50
         call HwU_book(l+ 21,'j3 y ' //cc(i),50,-6.d0,6.d0)
 
51
 
 
52
         call HwU_book(l+ 22,'H+j1 pT ' //cc(i),50,0.d0,400.d0)
 
53
         call HwU_book(l+ 23,'H+j1 pT ' //cc(i),50,0.d0,800.d0)
 
54
         call HwU_book(l+ 24,'H+j1 logpT ' //cc(i),50,0.d0,4.d0)
 
55
         call HwU_book(l+ 25,'H+j1 eta ' //cc(i),50,-6.d0,6.d0)
 
56
         call HwU_book(l+ 26,'H+j1 y ' //cc(i),50,-6.d0,6.d0)
 
57
 
 
58
         call HwU_book(l+ 27,'j1+j2 pT ' //cc(i),50,0.d0,400.d0)
 
59
         call HwU_book(l+ 28,'j1+j2 pT ' //cc(i),50,0.d0,800.d0)
 
60
         call HwU_book(l+ 29,'j1+j2 logpT ' //cc(i),50,0.d0,4.d0)
 
61
         call HwU_book(l+ 30,'j1+j2 eta ' //cc(i),50,-6.d0,6.d0)
 
62
         call HwU_book(l+ 31,'j1+j2 y ' //cc(i),50,-6.d0,6.d0)
 
63
 
 
64
         call HwU_book(l+ 32,'syst pT ' //cc(i),50,0.d0,400.d0)
 
65
         call HwU_book(l+ 33,'syst pT ' //cc(i),50,0.d0,800.d0)
 
66
         call HwU_book(l+ 34,'syst logpT ' //cc(i),50,0.d0,4.d0)
 
67
         call HwU_book(l+ 35,'syst eta ' //cc(i),50,-10.d0,10.d0)
 
68
         call HwU_book(l+ 36,'syst y ' //cc(i),50,-6.d0,6.d0)
 
69
 
 
70
         call HwU_book(l+ 37,'Dphi H-j1 ' //cc(i),50,0d0,pi)
 
71
         call HwU_book(l+ 38,'Dphi H-j2 ' //cc(i),50,0d0,pi)
 
72
         call HwU_book(l+ 39,'Dphi j1-j2 ' //cc(i),50,0d0,pi)
 
73
         
 
74
         call HwU_book(l+ 40,'DR H-j1 ' //cc(i),50,0d0,10.d0)
 
75
         call HwU_book(l+ 41,'DR H-j2 ' //cc(i),50,0d0,10.d0)
 
76
         call HwU_book(l+ 42,'DR j1-j2 ' //cc(i),50,0d0,10.d0)
 
77
 
 
78
         call HwU_book(l+ 43,'mj1j2 ' //cc(i),50,0d0,3000.d0)
 
79
 
 
80
c Nason-Oleari plots (hep-ph/0911.5299)
 
81
         call HwU_book(l+ 44,'|yj1-yj2| ' //cc(i),25,0.d0,10.d0)
 
82
         call HwU_book(l+ 45,'yj3_rel ' //cc(i),50,-6.d0,6.d0)
 
83
         call HwU_book(l+ 46,'njets ' //cc(i),100,-0.5d0,9.5d0)
 
84
         call HwU_book(l+ 47,'ptrel_j1 ' //cc(i),50,0.d0,200.d0)
 
85
         call HwU_book(l+ 48,'ptrel_j2 ' //cc(i),50,0.d0,200.d0)
 
86
         call HwU_book(l+ 49,'P-veto ' //cc(i), dble(nbinveto),vetomin
 
87
     &        ,vetomax)
 
88
         call HwU_book(l+ 50,'jveto pT ' //cc(i), dble(nbinveto),vetomin
 
89
     &        ,vetomax)
 
90
         call HwU_book(l+ 51,'jveto pT ' //cc(i), dble(nbinveto),
 
91
     &        vetomin,2d0*vetomax)
 
92
         call HwU_book(l+ 52,'jveto logpT ' //cc(i),50,0.d0,4.d0)
 
93
         call HwU_book(l+ 53,'jveto eta ' //cc(i),50,-6.d0,6.d0)
 
94
         call HwU_book(l+ 54,'jveto y ' //cc(i),50,-6.d0,6.d0)
 
95
 
 
96
      enddo
 
97
      return
 
98
      end
 
99
 
 
100
 
 
101
 
 
102
 
 
103
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
104
      subroutine analysis_end(dummy)
 
105
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
106
      implicit none
 
107
      double precision dummy
 
108
      call HwU_write_file
 
109
      return                
 
110
      end
 
111
 
 
112
 
 
113
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
114
      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
 
115
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
116
      implicit none
 
117
      include 'nexternal.inc'
 
118
      integer istatus(nexternal)
 
119
      integer iPDG(nexternal)
 
120
      double precision p(0:4,nexternal)
 
121
      double precision wgts(*)
 
122
      integer ibody
 
123
      double precision wgt,var
 
124
      integer i,j,k,kk,l
 
125
      double precision www,pQCD(0:3,nexternal),palg,rfj,sycut,yjmax
 
126
     $     ,pjet(0:3,nexternal),ptjet(nexternal),yjet(nexternal)
 
127
     $     ,etajet(nexternal),ptj_tag,deltay12,mj1j2min,ph(0:3),pj1(0:3)
 
128
     $     ,pj2(0:3),pj3(0:3),pjj(0:3),pHj(0:3),psyst(0:3),pjveto(0:3)
 
129
     $     ,ptH,etaH,yH,njdble,ptj1,etaj1,yj1,ptHj,etaHj,yHj,DphiHj1
 
130
     $     ,DRHj1,ptj2,etaj2,yj2,ptjj,etajj,yjj,ptsyst,etasyst,ysyst
 
131
     $     ,DphiHj2,Dphij1j2,DRHj2,DRj1j2,mj1j2,Dyj1j2,ptj3,etaj3,yj3
 
132
     $     ,yj3rel,chy1,shy1,chy1mo,chy2,shy2,chy2mo,ptrel_j1,ptrel_j2
 
133
     $     ,ppboost(0:3,nexternal),prel_j1(0:3),prel_j2(0:3)
 
134
     $     ,pj1boost(0:3),pj2boost(0:3),pt_veto,xsecup2
 
135
      logical pass_tag_cuts,flag
 
136
      integer nQCD,jet(nexternal),ij1y,ij2y,ij3y,njet,njety,ijveto
 
137
     $     ,ijvetoy,ij1,ij2,ij3
 
138
      double precision getptv4,getrapidityv4,getpseudorapv4,getdelphiv4
 
139
     $     ,getdrv4,getinvmv4,getmod
 
140
      external getptv4,getrapidityv4,getpseudorapv4,getdelphiv4,getdrv4
 
141
     $     ,getinvmv4,getmod
 
142
      real*8 vetomin, vetomax
 
143
      integer nbinveto
 
144
      common /to_veto_hist/vetomin,vetomax,nbinveto
 
145
      real*8 xd(1:3)
 
146
      data (xd(i),i=1,3)/0d0,0d0,1d0/
 
147
      if (nexternal.ne.6) then
 
148
         write (*,*) 'error #1 in analysis_fill: '/
 
149
     &        /'only for process "p p > h j j [QCD]"'
 
150
         stop 1
 
151
      endif
 
152
      if (.not. (abs(ipdg(1)).le.5 .or. ipdg(1).eq.21)) then
 
153
         write (*,*) 'error #2 in analysis_fill: '/
 
154
     &        /'only for process "p p > h j j [QCD]"'
 
155
         stop 1
 
156
      endif
 
157
      if (.not. (abs(ipdg(2)).le.5 .or. ipdg(2).eq.21)) then
 
158
         write (*,*) 'error #3 in analysis_fill: '/
 
159
     &        /'only for process "p p > h j j [QCD]"'
 
160
         stop 1
 
161
      endif
 
162
      if (.not. (abs(ipdg(4)).le.5 .or. ipdg(4).eq.21)) then
 
163
         write (*,*) 'error #4 in analysis_fill: '/
 
164
     &        /'only for process "p p > h j j [QCD]"'
 
165
         stop 1
 
166
      endif
 
167
      if (.not. (abs(ipdg(5)).le.5 .or. ipdg(5).eq.21)) then
 
168
         write (*,*) 'error #5 in analysis_fill: '/
 
169
     &        /'only for process "p p > h j j [QCD]"'
 
170
         stop 1
 
171
      endif
 
172
      if (.not. (abs(ipdg(6)).le.5 .or. ipdg(6).eq.21)) then
 
173
         write (*,*) 'error #6 in analysis_fill: '/
 
174
     &        /'only for process "p p > h j j [QCD]"'
 
175
         stop 1
 
176
      endif
 
177
      if (ipdg(3).ne.25) then
 
178
         write (*,*) 'error #7 in analysis_fill: '/
 
179
     &        /'only for process "p p > h j j [QCD]"'
 
180
         stop 1
 
181
      endif
 
182
c
 
183
 
 
184
c Put all (light) QCD partons in momentum array for jet clustering.
 
185
      nQCD=0
 
186
      do j=nincoming+1,nexternal
 
187
         if (abs(ipdg(j)).le.5 .or. ipdg(j).eq.21) then
 
188
            nQCD=nQCD+1
 
189
            do i=0,3
 
190
               pQCD(i,nQCD)=p(i,j) 
 
191
            enddo
 
192
         endif
 
193
      enddo
 
194
      
 
195
C---CLUSTER THE EVENT
 
196
      palg  = -1.d0
 
197
      rfj   = 0.4d0
 
198
      sycut = 20d0
 
199
      yjmax = 4.5d0
 
200
      do i=1,nexternal
 
201
         do j=0,3
 
202
            pjet(j,i)=0d0
 
203
         enddo
 
204
         ptjet(i)=0d0
 
205
         yjet(i)=0d0
 
206
         etajet(i)=0d0
 
207
         jet(i)=0
 
208
      enddo
 
209
      ij1y=0
 
210
      ij2y=0
 
211
      ij3y=0
 
212
      njet=0
 
213
      njety=0
 
214
      ijveto = 0
 
215
      ijvetoy = 0
 
216
 
 
217
 
 
218
c******************************************************************************
 
219
c     call FASTJET to get all the jets
 
220
c
 
221
c     INPUT:
 
222
c     input momenta:               pQCD(0:3,nexternal), energy is 0th component
 
223
c     number of input momenta:     nQCD
 
224
c     radius parameter:            rfj
 
225
c     minumum jet pt:              sycut
 
226
c     jet algorithm:               palg, 1.0=kt, 0.0=C/A, -1.0 = anti-kt
 
227
c
 
228
c     OUTPUT:
 
229
c     jet momenta:                             pjet(0:3,nexternal), E is 0th cmpnt
 
230
c     the number of jets (with pt > SYCUT):    njet
 
231
c     the jet for a given particle 'i':        jet(i),   note that this is
 
232
c     the particle in pQCD, which doesn't necessarily correspond to the particle
 
233
c     label in the process
 
234
c
 
235
      call amcatnlo_fastjetppgenkt(pQCD,nQCD,rfj,sycut,palg,pjet,njet
 
236
     $     ,jet)
 
237
c
 
238
c******************************************************************************
 
239
      do i=1,njet
 
240
         ptjet(i)=getptv4(pjet(0,i))
 
241
         if(i.gt.1)then
 
242
            if (ptjet(i).gt.ptjet(i-1)) then
 
243
               write (*,*) "Error 1: jets should be ordered in pt"
 
244
               stop
 
245
            endif
 
246
         endif
 
247
         yjet(i)=getrapidityv4(pjet(0,i))
 
248
         etajet(i)=getpseudorapv4(pjet(0,i))
 
249
c look for veto jet without y cuts
 
250
         if (i.gt.2.and.yjet(i).gt.min(yjet(1),yjet(2)).and.
 
251
     &        yjet(i).lt.max(yjet(1),yjet(2)).and.ijveto.eq.0) ijveto=i
 
252
 
 
253
C now look for jets within the rapidity cuts
 
254
         if (dabs(yjet(i)).lt.yjmax) then
 
255
            njety=njety+1
 
256
            if (ij1y.eq.0) then
 
257
               ij1y = i
 
258
            else if (ij2y.eq.0) then
 
259
               ij2y = i
 
260
            else if (ij3y.eq.0) then
 
261
               ij3y = i
 
262
            endif
 
263
c look for veto jet with y cuts
 
264
            if (ij3y.gt.0.and.
 
265
     &           yjet(i).gt.min(yjet(ij1y),yjet(ij2y)).and.
 
266
     &           yjet(i).lt.max(yjet(ij1y),yjet(ij2y)).and.ijvetoy.eq.0) 
 
267
     &           ijvetoy = i
 
268
         endif
 
269
      enddo
 
270
 
 
271
c Nason-Oleari cuts (hep-ph/0911.5299)
 
272
      ptj_tag  = 20d0
 
273
      deltay12 = 4.d0
 
274
      mj1j2min = 600d0
 
275
 
 
276
c this is the loop for w-o / w vbf cuts
 
277
      do i=1,2
 
278
      if(i.eq.1) then 
 
279
         ij1 = 1
 
280
         ij2 = 2
 
281
         ij3 = 3
 
282
      endif
 
283
      if(i.eq.2) then
 
284
         njet = njety
 
285
         ijveto = ijvetoy
 
286
         ij1 = ij1y
 
287
         ij2 = ij2y
 
288
         ij3 = ij3y
 
289
      endif
 
290
 
 
291
c Load momenta
 
292
      xsecup2=1d0
 
293
      do k=0,3
 
294
         pH(k)   =p(k,3)
 
295
         pj1(k)  =pjet(k,ij1)
 
296
         pj2(k)  =pjet(k,ij2)
 
297
         pj3(k)  =pjet(k,ij3)
 
298
         pjj(k)  =pjet(k,ij1)+pjet(k,ij2)
 
299
         pHj(k)  =pjet(k,ij1)+pH(k)
 
300
         psyst(k)=pjet(k,ij1)+pjet(k,ij2)+pH(k)
 
301
         pjveto(k)=pjet(k,ijveto)
 
302
      enddo
 
303
 
 
304
c Define observables
 
305
c Higgs
 
306
         ptH     = getptv4(pH)
 
307
         etaH    = getpseudorapv4(pH)
 
308
         yH      = getrapidityv4(pH)
 
309
         njdble  = dble(njet)
 
310
c At least one jet
 
311
      if(njet.ge.1)then
 
312
        ptj1    = getptv4(pj1)
 
313
        etaj1   = getpseudorapv4(pj1)
 
314
        yj1     = getrapidityv4(pj1)
 
315
        ptHj    = getptv4(pHj)
 
316
        etaHj   = getpseudorapv4(pHj)
 
317
        yHj     = getrapidityv4(pHj)
 
318
        DphiHj1 = getdelphiv4(pH,pj1)
 
319
        DRHj1   = getdrv4(pH,pj1)
 
320
      endif
 
321
c At least two jets
 
322
      if(njet.ge.2)then
 
323
        ptj2    = getptv4(pj2)
 
324
        etaj2   = getpseudorapv4(pj2)
 
325
        yj2     = getrapidityv4(pj2)
 
326
        ptjj    = getptv4(pjj)
 
327
        etajj   = getpseudorapv4(pjj)
 
328
        yjj     = getrapidityv4(pjj)
 
329
        ptsyst  = getptv4(psyst)
 
330
        etasyst = getpseudorapv4(psyst)
 
331
        ysyst   = getrapidityv4(psyst)
 
332
        DphiHj2 = getdelphiv4(pH,pj2)
 
333
        Dphij1j2= getdelphiv4(pj1,pj2)
 
334
        DRHj2   = getdrv4(pH,pj2)
 
335
        DRj1j2  = getdrv4(pj1,pj2)
 
336
        mj1j2   = getinvmv4(pjj)
 
337
        Dyj1j2  = abs(yj1-yj2)
 
338
      endif
 
339
c At least three jets
 
340
      if(njet.ge.3)then
 
341
        ptj3    = getptv4(pj3)
 
342
        etaj3   = getpseudorapv4(pj3)
 
343
        yj3     = getrapidityv4(pj3)
 
344
        yj3rel  = yj3-(yj1+yj2)/2d0
 
345
      endif
 
346
c
 
347
      chy1=cosh(yj1)
 
348
      shy1=sinh(yj1)
 
349
      chy1mo=chy1-1.d0
 
350
      chy2=cosh(yj2)
 
351
      shy2=sinh(yj2)
 
352
      chy2mo=chy2-1.d0
 
353
 
 
354
      call boostwdir2(chy1,shy1,chy1mo,xd,pj1,pj1boost)
 
355
      call boostwdir2(chy2,shy2,chy2mo,xd,pj2,pj2boost)
 
356
      ptrel_j1=0d0
 
357
      ptrel_j2=0d0
 
358
 
 
359
      pass_tag_cuts = njety.ge.2 .and.
 
360
     &                ptj1.ge.ptj_tag .and.
 
361
     &                ptj2.ge.ptj_tag .and.
 
362
     &                abs(yj1-yj2).ge.deltay12 .and.
 
363
     &                yj1*yj2.le.0d0 .and.
 
364
     &                mj1j2.ge.mj1j2min 
 
365
 
 
366
      if(i.eq.1) then 
 
367
         flag=.true.
 
368
      endif
 
369
 
 
370
      if(i.eq.2) then
 
371
         flag=pass_tag_cuts
 
372
      endif
 
373
 
 
374
 
 
375
      do j=1,nQCD
 
376
         if(njet.ge.1.and.jet(j).eq.1)then
 
377
           call boostwdir2(chy1,shy1,chy1mo,xd,pQCD(0,j),ppboost(0,j))
 
378
           call getwedge(ppboost(0,j),pj1boost,prel_j1)
 
379
           ptrel_j1=ptrel_j1+getmod(prel_j1)/getmod(pj1boost)
 
380
         elseif(njet.ge.2.and.jet(j).eq.2)then
 
381
           call boostwdir2(chy2,shy2,chy2mo,xd,pQCD(0,j),ppboost(0,j))
 
382
           call getwedge(ppboost(0,j),pj2boost,prel_j2)
 
383
           ptrel_j2=ptrel_j2+getmod(prel_j2)/getmod(pj2boost)
 
384
         endif
 
385
      enddo
 
386
 
 
387
      l=(i-1)*54
 
388
      if(flag)then
 
389
         call HwU_fill(l+  1,1d0,wgts)
 
390
         call HwU_fill(l+  2,ptH,wgts)
 
391
         call HwU_fill(l+  3,ptH,wgts)
 
392
         if(ptH.gt.0d0) call HwU_fill(l+  4,log10(ptH),wgts)
 
393
         call HwU_fill(l+  5,etaH,wgts)
 
394
         call HwU_fill(l+  6,yH,wgts)
 
395
         call HwU_fill(l+ 46,njdble,wgts)
 
396
 
 
397
         if(njet.ge.1)then
 
398
            call HwU_fill(l+  7,ptj1,wgts)
 
399
            call HwU_fill(l+  8,ptj1,wgts)
 
400
            if (ptj1.gt.0d0) call HwU_fill(l+  9,log10(ptj1),wgts)
 
401
            call HwU_fill(l+ 10,etaj1,wgts)
 
402
            call HwU_fill(l+ 11,yj1,wgts)
 
403
            call HwU_fill(l+ 22,ptHj,wgts)
 
404
            call HwU_fill(l+ 23,ptHj,wgts)
 
405
            if(ptHj.gt.0d0) call HwU_fill(l+ 24,log10(ptHj),wgts)
 
406
            call HwU_fill(l+ 25,etaHj,wgts)
 
407
            call HwU_fill(l+ 26,yHj,wgts)
 
408
            call HwU_fill(l+ 37,DphiHj1,wgts)
 
409
            call HwU_fill(l+ 40,DRHj1,wgts)
 
410
            call HwU_fill(l+ 47,ptrel_j1,wgts)
 
411
         endif 
 
412
         
 
413
         if(njet.ge.2)then
 
414
            call HwU_fill(l+ 12,ptj2,wgts)
 
415
            call HwU_fill(l+ 13,ptj2,wgts)
 
416
            if (ptj2.gt.0d0) call HwU_fill(l+ 14,log10(ptj2),wgts)
 
417
            call HwU_fill(l+ 15,etaj2,wgts)
 
418
            call HwU_fill(l+ 16,yj2,wgts)
 
419
            call HwU_fill(l+ 27,ptjj,wgts)
 
420
            call HwU_fill(l+ 28,ptjj,wgts)
 
421
            if(ptjj.gt.0d0) call HwU_fill(l+ 29,log10(ptjj),wgts)
 
422
            call HwU_fill(l+ 30,etajj,wgts)
 
423
            call HwU_fill(l+ 31,yjj,wgts)
 
424
            call HwU_fill(l+ 32,ptsyst,wgts)
 
425
            call HwU_fill(l+ 33,ptsyst,wgts)
 
426
            if(ptsyst.gt.0d0) call HwU_fill(l+ 34,log10(ptsyst),wgts)
 
427
            call HwU_fill(l+ 35,etasyst,wgts)
 
428
            call HwU_fill(l+ 36,ysyst,wgts)
 
429
            call HwU_fill(l+ 38,DphiHj2,wgts)
 
430
            call HwU_fill(l+ 39,Dphij1j2,wgts)
 
431
            call HwU_fill(l+ 41,DRHj2,wgts)
 
432
            call HwU_fill(l+ 42,DRj1j2,wgts)
 
433
            call HwU_fill(l+ 43,mj1j2,wgts)
 
434
            call HwU_fill(l+ 44,Dyj1j2,wgts)
 
435
            call HwU_fill(l+ 48,ptrel_j2,wgts)
 
436
         endif
 
437
         
 
438
         if(njet.ge.3)then
 
439
            call HwU_fill(l+ 17,ptj3,wgts)
 
440
            call HwU_fill(l+ 18,ptj3,wgts)
 
441
            if(ptj3.gt.0d0) call HwU_fill(l+ 19,log10(ptj3),wgts)
 
442
            call HwU_fill(l+ 20,etaj3,wgts)
 
443
            call HwU_fill(l+ 21,yj3,wgts)
 
444
            call HwU_fill(l+ 45,yj3rel,wgts)
 
445
         endif
 
446
         if (ijveto.gt.0) then
 
447
            pt_veto = getptv4(pjveto)
 
448
            do k=1,nbinveto
 
449
               if (pt_veto.gt. (vetomin+(vetomax-vetomin)*dble(k-1)
 
450
     $              /dble(nbinveto))) then
 
451
                  call HwU_fill(l+49, (vetomax-vetomin)* dble(k)
 
452
     $                 /dble(nbinveto)*0.99, wgts/xsecup2)
 
453
               endif
 
454
            enddo
 
455
            call HwU_fill(l+50,pt_veto,wgts)
 
456
            call HwU_fill(l+51,pt_veto,wgts)
 
457
            if (pt_veto.gt.0d0) call HwU_fill(l+52,dlog10(pt_veto),wgts)
 
458
            call HwU_fill(l+53,getpseudorapv4(pjveto),wgts)
 
459
            call HwU_fill(l+54,getrapidityv4(pjveto),wgts)
 
460
         endif
 
461
      endif
 
462
      enddo
 
463
 
 
464
 999  return      
 
465
      end
 
466
 
 
467
 
 
468
      function getrapidity(en,pl)
 
469
      implicit none
 
470
      real*8 getrapidity,en,pl,tiny,xplus,xminus,y
 
471
      parameter (tiny=1.d-8)
 
472
c
 
473
      xplus=en+pl
 
474
      xminus=en-pl
 
475
      if(xplus.gt.tiny.and.xminus.gt.tiny)then
 
476
        if( (xplus/xminus).gt.tiny.and.(xminus/xplus).gt.tiny )then
 
477
          y=0.5d0*log( xplus/xminus )
 
478
        else
 
479
          y=sign(1.d0,pl)*1.d8
 
480
        endif
 
481
      else
 
482
        y=sign(1.d0,pl)*1.d8
 
483
      endif
 
484
      getrapidity=y
 
485
      return
 
486
      end
 
487
 
 
488
 
 
489
      function getpseudorap(en,ptx,pty,pl)
 
490
      implicit none
 
491
      real*8 getpseudorap,en,ptx,pty,pl,tiny,pt,eta,th
 
492
      parameter (tiny=1.d-5)
 
493
c
 
494
      pt=sqrt(ptx**2+pty**2)
 
495
      if(pt.lt.tiny.and.abs(pl).lt.tiny)then
 
496
        eta=sign(1.d0,pl)*1.d8
 
497
      else
 
498
        th=atan2(pt,pl)
 
499
        eta=-log(tan(th/2.d0))
 
500
      endif
 
501
      getpseudorap=eta
 
502
      return
 
503
      end
 
504
 
 
505
 
 
506
      function getinvm(en,ptx,pty,pl)
 
507
      implicit none
 
508
      real*8 getinvm,en,ptx,pty,pl,tiny,tmp
 
509
      parameter (tiny=1.d-5)
 
510
c
 
511
      tmp=en**2-ptx**2-pty**2-pl**2
 
512
      if(tmp.gt.0.d0)then
 
513
        tmp=sqrt(tmp)
 
514
      elseif(tmp.gt.-tiny)then
 
515
        tmp=0.d0
 
516
      else
 
517
        write(*,*)'Attempt to compute a negative mass'
 
518
        stop
 
519
      endif
 
520
      getinvm=tmp
 
521
      return
 
522
      end
 
523
 
 
524
 
 
525
      function getdelphi(ptx1,pty1,ptx2,pty2)
 
526
      implicit none
 
527
      real*8 getdelphi,ptx1,pty1,ptx2,pty2,tiny,pt1,pt2,tmp
 
528
      parameter (tiny=1.d-5)
 
529
c
 
530
      pt1=sqrt(ptx1**2+pty1**2)
 
531
      pt2=sqrt(ptx2**2+pty2**2)
 
532
      if(pt1.ne.0.d0.and.pt2.ne.0.d0)then
 
533
        tmp=ptx1*ptx2+pty1*pty2
 
534
        tmp=tmp/(pt1*pt2)
 
535
        if(abs(tmp).gt.1.d0+tiny)then
 
536
          write(*,*)'Cosine larger than 1'
 
537
          stop
 
538
        elseif(abs(tmp).ge.1.d0)then
 
539
          tmp=sign(1.d0,tmp)
 
540
        endif
 
541
        tmp=acos(tmp)
 
542
      else
 
543
        tmp=1.d8
 
544
      endif
 
545
      getdelphi=tmp
 
546
      return
 
547
      end
 
548
 
 
549
 
 
550
      function getdr(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
 
551
      implicit none
 
552
      real*8 getdr,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
 
553
     # getpseudorap,getdelphi
 
554
c
 
555
      deta=getpseudorap(en1,ptx1,pty1,pl1)-
 
556
     #     getpseudorap(en2,ptx2,pty2,pl2)
 
557
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
 
558
      getdr=sqrt(dphi**2+deta**2)
 
559
      return
 
560
      end
 
561
 
 
562
 
 
563
      function getdry(en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2)
 
564
      implicit none
 
565
      real*8 getdry,en1,ptx1,pty1,pl1,en2,ptx2,pty2,pl2,deta,dphi,
 
566
     # getrapidity,getdelphi
 
567
c
 
568
      deta=getrapidity(en1,pl1)-
 
569
     #     getrapidity(en2,pl2)
 
570
      dphi=getdelphi(ptx1,pty1,ptx2,pty2)
 
571
      getdry=sqrt(dphi**2+deta**2)
 
572
      return
 
573
      end
 
574
 
 
575
 
 
576
      function getptv(p)
 
577
      implicit none
 
578
      real*8 getptv,p(5)
 
579
c
 
580
      getptv=sqrt(p(1)**2+p(2)**2)
 
581
      return
 
582
      end
 
583
 
 
584
 
 
585
      function getpseudorapv(p)
 
586
      implicit none
 
587
      real*8 getpseudorapv,p(5)
 
588
      real*8 getpseudorap
 
589
c
 
590
      getpseudorapv=getpseudorap(p(4),p(1),p(2),p(3))
 
591
      return
 
592
      end
 
593
 
 
594
 
 
595
      function getrapidityv(p)
 
596
      implicit none
 
597
      real*8 getrapidityv,p(5)
 
598
      real*8 getrapidity
 
599
c
 
600
      getrapidityv=getrapidity(p(4),p(3))
 
601
      return
 
602
      end
 
603
 
 
604
 
 
605
      function getdrv(p1,p2)
 
606
      implicit none
 
607
      real*8 getdrv,p1(5),p2(5)
 
608
      real*8 getdr
 
609
c
 
610
      getdrv=getdr(p1(4),p1(1),p1(2),p1(3),
 
611
     #             p2(4),p2(1),p2(2),p2(3))
 
612
      return
 
613
      end
 
614
 
 
615
 
 
616
      function getinvmv(p)
 
617
      implicit none
 
618
      real*8 getinvmv,p(5)
 
619
      real*8 getinvm
 
620
c
 
621
      getinvmv=getinvm(p(4),p(1),p(2),p(3))
 
622
      return
 
623
      end
 
624
 
 
625
 
 
626
      function getdelphiv(p1,p2)
 
627
      implicit none
 
628
      real*8 getdelphiv,p1(5),p2(5)
 
629
      real*8 getdelphi
 
630
c
 
631
      getdelphiv=getdelphi(p1(1),p1(2),
 
632
     #                     p2(1),p2(2))
 
633
      return
 
634
      end
 
635
 
 
636
 
 
637
      function getptv4(p)
 
638
      implicit none
 
639
      real*8 getptv4,p(0:3)
 
640
c
 
641
      getptv4=sqrt(p(1)**2+p(2)**2)
 
642
      return
 
643
      end
 
644
 
 
645
 
 
646
      function getpseudorapv4(p)
 
647
      implicit none
 
648
      real*8 getpseudorapv4,p(0:3)
 
649
      real*8 getpseudorap
 
650
c
 
651
      getpseudorapv4=getpseudorap(p(0),p(1),p(2),p(3))
 
652
      return
 
653
      end
 
654
 
 
655
 
 
656
      function getrapidityv4(p)
 
657
      implicit none
 
658
      real*8 getrapidityv4,p(0:3)
 
659
      real*8 getrapidity
 
660
c
 
661
      getrapidityv4=getrapidity(p(0),p(3))
 
662
      return
 
663
      end
 
664
 
 
665
 
 
666
      function getdrv4(p1,p2)
 
667
      implicit none
 
668
      real*8 getdrv4,p1(0:3),p2(0:3)
 
669
      real*8 getdr
 
670
c
 
671
      getdrv4=getdr(p1(0),p1(1),p1(2),p1(3),
 
672
     #              p2(0),p2(1),p2(2),p2(3))
 
673
      return
 
674
      end
 
675
 
 
676
 
 
677
      function getinvmv4(p)
 
678
      implicit none
 
679
      real*8 getinvmv4,p(0:3)
 
680
      real*8 getinvm
 
681
c
 
682
      getinvmv4=getinvm(p(0),p(1),p(2),p(3))
 
683
      return
 
684
      end
 
685
 
 
686
 
 
687
      function getdelphiv4(p1,p2)
 
688
      implicit none
 
689
      real*8 getdelphiv4,p1(0:3),p2(0:3)
 
690
      real*8 getdelphi
 
691
c
 
692
      getdelphiv4=getdelphi(p1(1),p1(2),
 
693
     #                      p2(1),p2(2))
 
694
      return
 
695
      end
 
696
 
 
697
 
 
698
      function getcosv4(q1,q2)
 
699
      implicit none
 
700
      real*8 getcosv4,q1(0:3),q2(0:3)
 
701
      real*8 xnorm1,xnorm2,tmp
 
702
c
 
703
      if(q1(0).lt.0.d0.or.q2(0).lt.0.d0)then
 
704
        getcosv4=-1.d10
 
705
        return
 
706
      endif
 
707
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
 
708
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
 
709
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
 
710
        tmp=-1.d10
 
711
      else
 
712
        tmp=q1(1)*q2(1)+q1(2)*q2(2)+q1(3)*q2(3)
 
713
        tmp=tmp/(xnorm1*xnorm2)
 
714
        if(abs(tmp).gt.1.d0.and.abs(tmp).le.1.001d0)then
 
715
          tmp=sign(1.d0,tmp)
 
716
        elseif(abs(tmp).gt.1.001d0)then
 
717
          write(*,*)'Error in getcosv4',tmp
 
718
          stop
 
719
        endif
 
720
      endif
 
721
      getcosv4=tmp
 
722
      return
 
723
      end
 
724
 
 
725
 
 
726
 
 
727
      function getmod(p)
 
728
      implicit none
 
729
      double precision p(0:3),getmod
 
730
 
 
731
      getmod=sqrt(p(1)**2+p(2)**2+p(3)**2)
 
732
 
 
733
      return
 
734
      end
 
735
 
 
736
 
 
737
 
 
738
      subroutine getperpenv4(q1,q2,qperp)
 
739
c Normal to the plane defined by \vec{q1},\vec{q2}
 
740
      implicit none
 
741
      real*8 q1(0:3),q2(0:3),qperp(0:3)
 
742
      real*8 xnorm1,xnorm2
 
743
      integer i
 
744
c
 
745
      xnorm1=sqrt(q1(1)**2+q1(2)**2+q1(3)**2)
 
746
      xnorm2=sqrt(q2(1)**2+q2(2)**2+q2(3)**2)
 
747
      if(xnorm1.lt.1.d-6.or.xnorm2.lt.1.d-6)then
 
748
        do i=1,4
 
749
          qperp(i)=-1.d10
 
750
        enddo
 
751
      else
 
752
        qperp(1)=q1(2)*q2(3)-q1(3)*q2(2)
 
753
        qperp(2)=q1(3)*q2(1)-q1(1)*q2(3)
 
754
        qperp(3)=q1(1)*q2(2)-q1(2)*q2(1)
 
755
        do i=1,3
 
756
          qperp(i)=qperp(i)/(xnorm1*xnorm2)
 
757
        enddo
 
758
        qperp(0)=1.d0
 
759
      endif
 
760
      return
 
761
      end
 
762
 
 
763
 
 
764
 
 
765
 
 
766
      subroutine getwedge(p1,p2,pout)
 
767
      implicit none
 
768
      real*8 p1(0:3),p2(0:3),pout(0:3)
 
769
 
 
770
      pout(1)=p1(2)*p2(3)-p1(3)*p2(2)
 
771
      pout(2)=p1(3)*p2(1)-p1(1)*p2(3)
 
772
      pout(3)=p1(1)*p2(2)-p1(2)*p2(1)
 
773
      pout(0)=0d0
 
774
 
 
775
      return
 
776
      end
 
777