← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_407857/PROC_407857/SubProcesses/unwgt.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_407857/PROC_407857/SubProcesses/unwgt.f
 
1
      subroutine Zoom_Event(wgt,p)
 
2
C**************************************************************************
 
3
C     Determines if region needs to be investigated in case of large
 
4
c     weight events.
 
5
C**************************************************************************
 
6
      IMPLICIT NONE
 
7
c
 
8
c     Constant
 
9
c
 
10
      integer    max_zoom
 
11
      parameter (max_zoom=2000)
 
12
      include 'genps.inc'
 
13
      include 'nexternal.inc'
 
14
 
 
15
c
 
16
c     Arguments
 
17
c
 
18
      double precision wgt, p(0:3,nexternal)
 
19
c
 
20
c     Local
 
21
c
 
22
      double precision xstore(2),gstore,qstore(2)
 
23
      double precision trunc_wgt, xsum, wstore,pstore(0:3,nexternal)
 
24
      integer ix, i,j
 
25
 
 
26
C     
 
27
C     GLOBAL
 
28
C
 
29
      double precision twgt, maxwgt,swgt(maxevents)
 
30
      integer                             lun, nw, itmin
 
31
      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
32
      integer nzoom
 
33
      double precision  tx(1:3,maxinvar)
 
34
      common/to_xpoints/tx, nzoom
 
35
      double precision xzoomfact
 
36
      common/to_zoom/  xzoomfact
 
37
      include 'run.inc'
 
38
      include 'coupl.inc'
 
39
c
 
40
c     DATA
 
41
c
 
42
      data trunc_wgt /-1d0/
 
43
      data xsum/0d0/
 
44
      data wstore /0d0/
 
45
      save ix, pstore, wstore, xstore, gstore, qstore
 
46
c-----
 
47
c  Begin Code
 
48
c-----
 
49
      if (trunc_Wgt .lt. 0d0 .and. twgt .gt. 0d0) then
 
50
         write(*,*) 'Selecting zoom level', twgt*500, wgt
 
51
      endif
 
52
      if (twgt .lt. 0d0) then
 
53
         write(*,*) 'Resetting zoom iteration', twgt
 
54
         twgt = -twgt
 
55
         trunc_wgt = twgt * 500d0
 
56
      endif
 
57
      if (nw .eq. 0) then
 
58
         trunc_wgt = twgt * 500d0
 
59
      endif
 
60
      trunc_wgt=max(trunc_wgt, twgt*500d0)
 
61
      if (nzoom .eq. 0 .and. trunc_wgt .gt. 0d0 ) then
 
62
         if (wgt .gt. trunc_wgt) then
 
63
            write(*,*) 'Zooming on large event ',wgt / trunc_wgt
 
64
            wstore=wgt
 
65
            do i=1,nexternal
 
66
               do j=0,3
 
67
                  pstore(j,i) = p(j,i)
 
68
               enddo
 
69
            enddo
 
70
            do i=1,2
 
71
               xstore(i)=xbk(i)
 
72
               qstore(i)=q2fact(i)
 
73
            enddo
 
74
            gstore=g
 
75
            xsum = wgt
 
76
            nzoom = max_zoom
 
77
            wgt=0d0
 
78
            ix = 1
 
79
         endif
 
80
      elseif (trunc_wgt .gt. 0 .and. wgt .gt. 0d0) then
 
81
         xsum = xsum + wgt
 
82
         if (nzoom .gt. 1) wgt = 0d0
 
83
         ix = ix + 1
 
84
      endif
 
85
      if (xsum .ne. 0d0 .and. nzoom .le. 1) then
 
86
         if (wgt .gt. 0d0) then
 
87
c            xzoomfact = xsum/real(max_zoom) / wgt !Store avg wgt
 
88
            xzoomfact = wstore / wgt  !Store large wgt
 
89
         else
 
90
            xzoomfact = -xsum/real(max_zoom)
 
91
         endif
 
92
         wgt = max(xsum/real(max_zoom),trunc_wgt)  !Don't make smaller then truncated wgt
 
93
         do i=1,nexternal
 
94
            do j=0,3
 
95
               p(j,i) = pstore(j,i)
 
96
            enddo
 
97
         enddo
 
98
         do i=1,2
 
99
            xbk(i)=xstore(i)
 
100
            q2fact(i)=qstore(i)
 
101
         enddo
 
102
         g=gstore
 
103
         write(*,'(a,2e15.3,2f15.3, i8)') 'Stored wgt ',
 
104
     $            wgt/trunc_wgt, wstore, wstore/wgt, real(ix)/max_zoom, ix
 
105
         trunc_wgt = max(trunc_wgt, wgt)
 
106
         xsum = 0d0
 
107
         nzoom = 0
 
108
      endif
 
109
      end
 
110
 
 
111
      subroutine clear_events
 
112
c-------------------------------------------------------------------
 
113
c     delete all events thus far, start from scratch
 
114
c------------------------------------------------------------------
 
115
      implicit none
 
116
c
 
117
c     Constants
 
118
c
 
119
      include 'genps.inc'
 
120
      include 'nexternal.inc'
 
121
c
 
122
c     Global
 
123
c
 
124
      integer iseed, nover, nstore
 
125
C     
 
126
C     GLOBAL
 
127
C
 
128
      double precision twgt, maxwgt,swgt(maxevents)
 
129
      integer                             lun, nw, itmin
 
130
      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
131
 
 
132
c-----
 
133
c  Begin Code
 
134
c-----
 
135
c      write(*,*) 'storing Events'
 
136
      call store_events
 
137
      rewind(lun)
 
138
      nw = 0
 
139
      maxwgt = 0d0
 
140
      end
 
141
 
 
142
      SUBROUTINE unwgt(px,wgt,numproc)
 
143
C**************************************************************************
 
144
C     Determines if event should be written based on its weight
 
145
C**************************************************************************
 
146
      IMPLICIT NONE
 
147
c
 
148
c     Constants
 
149
c
 
150
      include 'genps.inc'
 
151
      include 'nexternal.inc'
 
152
c
 
153
c     Arguments
 
154
c
 
155
      double precision px(0:3,nexternal),wgt
 
156
      integer numproc
 
157
c
 
158
c     Local
 
159
c
 
160
      integer idum, i,j
 
161
      double precision uwgt,yran,fudge, p(0:3,nexternal), xwgt
 
162
C     
 
163
C     GLOBAL
 
164
C
 
165
      double precision twgt, maxwgt,swgt(maxevents)
 
166
      integer                             lun, nw, itmin
 
167
      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
168
 
 
169
      double precision    matrix
 
170
      common/to_maxmatrix/matrix
 
171
 
 
172
      logical               zooming
 
173
      common /to_zoomchoice/zooming
 
174
 
 
175
c
 
176
c     External
 
177
c
 
178
      real xran1
 
179
      external xran1
 
180
c
 
181
c     Data
 
182
c
 
183
      data idum/-1/
 
184
      data yran/1d0/
 
185
      data fudge/10d0/
 
186
C-----
 
187
C  BEGIN CODE
 
188
C-----
 
189
      if (twgt .ge. 0d0) then
 
190
         do i=1,nexternal
 
191
            do j=0,3
 
192
               p(j,i)=px(j,i)
 
193
            enddo
 
194
         enddo
 
195
         xwgt = abs(wgt)
 
196
         if (zooming) call zoom_event(xwgt,P)
 
197
         if (xwgt .eq. 0d0) return
 
198
         yran = xran1(idum)
 
199
         if (xwgt .gt. twgt*fudge*yran) then
 
200
            uwgt = max(xwgt,twgt*fudge)
 
201
c           Set sign of uwgt to sign of wgt
 
202
            uwgt = dsign(uwgt,wgt)
 
203
            if (twgt .gt. 0) uwgt=uwgt/twgt/fudge
 
204
c            call write_event(p,uwgt)
 
205
c            write(29,'(2e15.5)') matrix,wgt
 
206
c $B$ S-COMMENT_C $B$
 
207
            call write_leshouche(p,uwgt,numproc)
 
208
         elseif (xwgt .gt. 0d0 .and. nw .lt. 5) then
 
209
            call write_leshouche(p,wgt/twgt*1d-6,numproc)
 
210
c $E$ S-COMMENT_C $E$
 
211
         endif
 
212
         maxwgt=max(maxwgt,xwgt)
 
213
      endif
 
214
      end
 
215
 
 
216
      subroutine store_events()
 
217
C**************************************************************************
 
218
C     Takes events from scratch file (lun) and writes them to a permanent
 
219
c     file  events.dat
 
220
C**************************************************************************
 
221
      IMPLICIT NONE
 
222
c
 
223
c     Constants
 
224
c
 
225
      include 'genps.inc'
 
226
      include 'nexternal.inc'
 
227
      double precision trunc_max
 
228
      parameter       (trunc_max = 0.01)  !Maximum % cross section to truncate
 
229
c
 
230
c     Arguments
 
231
c
 
232
c
 
233
c     Local
 
234
c
 
235
      integer i, lunw, ic(7,2*nexternal-3), n, j
 
236
      logical done
 
237
      double precision wgt,p(0:4,2*nexternal-3)
 
238
      double precision xsec,xsecabs,xerr,xscale,xtot
 
239
      double precision xsum, xover
 
240
      double precision target_wgt,orig_Wgt(maxevents)
 
241
      logical store_event(maxevents)
 
242
      integer iseed, nover, nstore
 
243
      double precision scale,aqcd,aqed
 
244
      integer ievent
 
245
      character*300 buff
 
246
C     
 
247
C     GLOBAL
 
248
C
 
249
      double precision twgt, maxwgt,swgt(maxevents)
 
250
      integer                             lun, nw, itmin
 
251
      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
252
 
 
253
      integer                   neventswritten
 
254
      common /to_eventswritten/ neventswritten
 
255
c      save neventswritten
 
256
 
 
257
      integer ngroup
 
258
      common/to_group/ngroup
 
259
 
 
260
c
 
261
c     external
 
262
c
 
263
      real xran1
 
264
 
 
265
      data iseed/0/
 
266
      data neventswritten/0/
 
267
C-----
 
268
C  BEGIN CODE
 
269
C-----
 
270
c
 
271
c     First scale all of the events to the total cross section
 
272
c
 
273
      if (nw .le. 0) return
 
274
      call sample_result(xsecabs,xsec,xerr,itmin)
 
275
      if (xsecabs .le. 0) return   !Fix by TS 12/3/2010
 
276
      xtot=0
 
277
      call dsort(nw, swgt)
 
278
      do i=1,nw
 
279
         xtot=xtot+dabs(swgt(i))
 
280
      enddo
 
281
c
 
282
c     Determine minimum target weight given truncation parameter
 
283
c
 
284
      xsum = 0d0
 
285
      i = nw
 
286
      do while (xsum-dabs(swgt(i))*(nw-i) .lt. xtot*trunc_max .and. i .gt. 2)
 
287
         xsum = xsum + dabs(swgt(i))
 
288
         i = i-1
 
289
      enddo
 
290
      if (i .lt. nw) i=i+1
 
291
      target_wgt = dabs(swgt(i))
 
292
c
 
293
c     Select events which will be written
 
294
c
 
295
      xsum = 0d0
 
296
      nstore = 0
 
297
      rewind(lun)
 
298
      done = .false. 
 
299
      do i=1,nw
 
300
         if (.not. done) then
 
301
            call read_event(lun,P,wgt,n,ic,ievent,scale,aqcd,aqed,buff,done)
 
302
         else
 
303
            wgt = 0d0
 
304
         endif
 
305
         if (dabs(wgt) .gt. target_wgt*xran1(iseed)) then
 
306
            xsum=xsum+max(dabs(wgt),target_Wgt)
 
307
            store_event(i)=.true.
 
308
            nstore=nstore+1
 
309
         else
 
310
            store_event(i) = .false.
 
311
         endif
 
312
      enddo
 
313
      xscale = xsecabs/xsum
 
314
      target_wgt = target_wgt*xscale
 
315
      rewind(lun)
 
316
c     JA 8/17/2011 Don't check for previously stored events
 
317
c      if (nstore .le. neventswritten) then
 
318
c         write(*,*) 'No improvement in events',nstore, neventswritten
 
319
c         return
 
320
c      endif
 
321
      lunw = 25
 
322
      open(unit = lunw, file='events.lhe', status='unknown')
 
323
      done = .false.
 
324
      i=0      
 
325
      xtot = 0
 
326
      xover = 0
 
327
      nover = 0
 
328
      do j=1,nw
 
329
         if (.not. done) then
 
330
            call read_event(lun,P,wgt,n,ic,ievent,scale,aqcd,aqed,buff,done)
 
331
         else
 
332
            write(*,*) 'Error done early',j,nw
 
333
         endif
 
334
         if (store_event(j) .and. .not. done) then
 
335
            wgt=wgt*xscale
 
336
            wgt = dsign(max(dabs(wgt), target_wgt),wgt)
 
337
            if (dabs(wgt) .gt. target_wgt) then
 
338
               xover = xover + dabs(wgt) - target_wgt
 
339
               nover = nover+1
 
340
            endif
 
341
            xtot = xtot + dabs(wgt)
 
342
            i=i+1
 
343
            call write_Event(lunw,p,wgt,n,ic,ngroup,scale,aqcd,aqed,buff)
 
344
         endif
 
345
      enddo
 
346
      write(*,*) 'Found ',nw,' events.'
 
347
      write(*,*) 'Wrote ',i ,' events.'
 
348
      write(*,*) 'Actual xsec ',xsec
 
349
      write(*,*) 'Correct abs xsec ',xsecabs
 
350
      write(*,*) 'Event xsec ', xtot
 
351
      write(*,*) 'Events wgts > 1: ', nover
 
352
      write(*,*) '% Cross section > 1: ',xover, xover/xtot*100.
 
353
      neventswritten = i
 
354
 99   close(lunw)
 
355
c      close(lun)
 
356
      end
 
357
 
 
358
      SUBROUTINE write_leshouche(p,wgt,numproc)
 
359
C**************************************************************************
 
360
C     Writes out information for event
 
361
C**************************************************************************
 
362
      IMPLICIT NONE
 
363
c
 
364
c     Constants
 
365
c
 
366
      double precision zero
 
367
      parameter       (ZERO = 0d0)
 
368
      include 'genps.inc'
 
369
      include 'nexternal.inc'
 
370
      include 'maxamps.inc'
 
371
      include 'message.inc'
 
372
      include 'cluster.inc'
 
373
      include 'run.inc'
 
374
c
 
375
c     Arguments
 
376
c
 
377
      double precision p(0:3,nexternal),wgt
 
378
      integer numproc
 
379
c
 
380
c     Local
 
381
c
 
382
      integer i,j,k
 
383
      double precision xtarget,targetamp(maxflow)
 
384
      integer ic, nc, jpart(7,-nexternal+3:2*nexternal-3)
 
385
      integer ida(2),ito(-nexternal+3:nexternal),ns,nres,ires,icloop
 
386
      integer iseed
 
387
      double precision pboost(0:3),pb(0:4,-nexternal+3:2*nexternal-3),eta
 
388
      double precision ptcltmp(nexternal), pdum(0:3)
 
389
 
 
390
      integer idup(nexternal,maxproc,maxsproc)
 
391
      integer mothup(2,nexternal)
 
392
      integer icolup(2,nexternal,maxflow,maxsproc)
 
393
 
 
394
      integer isym(nexternal,99), nsym, jsym
 
395
 
 
396
      double precision sscale,aaqcd,aaqed
 
397
      integer ievent,npart
 
398
      logical flip
 
399
 
 
400
      real ran1
 
401
      external ran1
 
402
 
 
403
      character*300 buff
 
404
 
 
405
C     
 
406
C     GLOBAL
 
407
C
 
408
      double precision twgt, maxwgt,swgt(maxevents)
 
409
      integer                             lun, nw, itmin
 
410
      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
411
 
 
412
      integer          IPSEL
 
413
      COMMON /SubProc/ IPSEL
 
414
 
 
415
      Double Precision amp2(maxamps), jamp2(0:maxflow)
 
416
      common/to_amps/  amp2,       jamp2
 
417
 
 
418
      character*101       hel_buf
 
419
      common/to_helicity/hel_buf
 
420
 
 
421
      integer           mincfig, maxcfig
 
422
      common/to_configs/mincfig, maxcfig
 
423
 
 
424
      integer ngroup
 
425
      common/to_group/ngroup
 
426
 
 
427
      double precision stot,m1,m2
 
428
      common/to_stot/stot,m1,m2
 
429
c
 
430
c     Data
 
431
c
 
432
      include 'leshouche.inc'
 
433
      data iseed /0/
 
434
 
 
435
      double precision pmass(nexternal)
 
436
      common/to_mass/  pmass
 
437
 
 
438
c      integer ncols,ncolflow(maxamps),ncolalt(maxamps)
 
439
c      common/to_colstats/ncols,ncolflow,ncolalt,ic
 
440
c      data ncolflow/maxamps*0/
 
441
c      data ncolalt/maxamps*0/
 
442
 
 
443
      include 'symswap.inc'
 
444
      include 'coupl.inc'
 
445
C-----
 
446
C  BEGIN CODE
 
447
C-----
 
448
      
 
449
      if (nw .ge. maxevents) return
 
450
 
 
451
c     Store weight for event
 
452
      nw = nw+1
 
453
      swgt(nw)=wgt
 
454
c
 
455
c     In case of identical particles symmetry, choose assignment
 
456
c
 
457
      xtarget = ran1(iseed)*nsym
 
458
      jsym = 1
 
459
      do while (xtarget .gt. jsym .and. jsym .lt. nsym)
 
460
         jsym = jsym+1
 
461
      enddo
 
462
c
 
463
c     Fill jpart color and particle info
 
464
c
 
465
      do i=1,nexternal
 
466
         jpart(1,isym(i,jsym)) = idup(i,ipsel,numproc)
 
467
         jpart(2,isym(i,jsym)) = mothup(1,i)
 
468
         jpart(3,isym(i,jsym)) = mothup(2,i)
 
469
c        Color info is filled in mothup
 
470
         jpart(4,isym(i,jsym)) = 0
 
471
         jpart(5,isym(i,jsym)) = 0
 
472
         jpart(6,isym(i,jsym)) = 1
 
473
      enddo
 
474
      do i=1,nincoming
 
475
         jpart(6,isym(i,jsym))=-1
 
476
      enddo
 
477
 
 
478
c   Set helicities
 
479
c      write(*,*) 'Getting helicity',hel_buf(1:50)
 
480
      read(hel_buf,'(20i5)') (jpart(7,isym(i, jsym)),i=1,nexternal)
 
481
c      write(*,*) 'ihel',jpart(7,1),jpart(7,2)
 
482
 
 
483
c   Fix ordering of ptclus
 
484
      do i=1,nexternal
 
485
        ptcltmp(isym(i,jsym)) = ptclus(i)
 
486
      enddo
 
487
      do i=1,nexternal
 
488
        ptclus(i) = ptcltmp(i)
 
489
      enddo
 
490
 
 
491
c     Check if we have flipped particle 1 and 2, and flip back
 
492
      flip = .false.
 
493
      if (p(3,1).lt.0) then
 
494
         do j=0,3
 
495
            pdum(j)=p(j,1)
 
496
            p(j,1)=p(j,2)
 
497
            p(j,2)=pdum(j)
 
498
         enddo
 
499
         flip = .true.
 
500
      endif
 
501
 
 
502
c
 
503
c     Boost momentum to lab frame
 
504
c
 
505
      pboost(0)=1d0
 
506
      pboost(1)=0d0
 
507
      pboost(2)=0d0
 
508
      pboost(3)=0d0
 
509
      if (nincoming.eq.2) then
 
510
         if (xbk(1) .gt. 0d0 .and. xbk(1) .lt. 1d0 .and.
 
511
     $       xbk(2) .gt. 0d0 .and. xbk(2) .lt. 1d0) then
 
512
            eta = sqrt(xbk(1)*ebeam(1)/
 
513
     $                 (xbk(2)*ebeam(2)))
 
514
            pboost(0)=p(0,1)*(eta + 1d0/eta)
 
515
            pboost(3)=p(0,1)*(eta - 1d0/eta)
 
516
         else if (xbk(1) .gt. 0d0 .and. xbk(1) .lt. 1d0 .and.
 
517
     $       xbk(2) .eq. 1d0) then
 
518
            pboost(0)=xbk(1)*ebeam(1)+ebeam(2)
 
519
            pboost(3)=xbk(1)*ebeam(1)-
 
520
     $           sqrt(ebeam(2)**2-m2**2)
 
521
         else if (xbk(1) .eq. 1d0 .and.
 
522
     $       xbk(2) .gt. 0d0 .and. xbk(2) .lt. 1d0) then
 
523
            pboost(0)=ebeam(1)+xbk(2)*ebeam(2)
 
524
            pboost(3)=sqrt(ebeam(1)**2-m1**2)-
 
525
     $           xbk(2)*ebeam(2)
 
526
         else if (xbk(1) .eq. 1d0 .and. xbk(2) .eq. 1d0) then
 
527
            pboost(0)=ebeam(1)+ebeam(2)
 
528
            pboost(3)=sqrt(ebeam(1)**2-m1**2)-
 
529
     $           sqrt(ebeam(2)**2-m2**2)
 
530
         else
 
531
            write(*,*) 'Warning bad x1 or x2 in write_leshouche',
 
532
     $           xbk(1),xbk(2)
 
533
         endif
 
534
      endif
 
535
      do j=1,nexternal
 
536
         call boostx(p(0,j),pboost,pb(0,isym(j,jsym)))
 
537
c      Add mass information in pb(4)
 
538
         pb(4,isym(j,jsym))=pmass(j)
 
539
      enddo
 
540
 
 
541
c
 
542
c     Add info on resonant mothers
 
543
c
 
544
      call addmothers(ipsel,jpart,pb,isym,jsym,sscale,aaqcd,aaqed,buff,
 
545
     $                npart,numproc)
 
546
 
 
547
c     Need to flip after addmothers, since color might get overwritten
 
548
      if (flip) then
 
549
         do i=1,7
 
550
            j=jpart(i,1)
 
551
            jpart(i,1)=jpart(i,2)
 
552
            jpart(i,2)=j
 
553
         enddo
 
554
         ptcltmp(1)=ptclus(1)
 
555
         ptclus(1)=ptclus(2)
 
556
         ptclus(2)=ptcltmp(1)
 
557
      endif
 
558
 
 
559
c
 
560
c     Write events to lun
 
561
c
 
562
c      write(*,*) 'Writing event'
 
563
      if(q2fact(1).gt.0.and.q2fact(2).gt.0)then
 
564
         sscale = sqrt(max(q2fact(1),q2fact(2)))
 
565
      else if(q2fact(1).gt.0)then
 
566
         sscale = sqrt(q2fact(1))
 
567
      else if(q2fact(2).gt.0)then
 
568
         sscale = sqrt(q2fact(2))
 
569
      else
 
570
         sscale = 0d0
 
571
      endif
 
572
      aaqcd = g*g/4d0/3.1415926d0
 
573
      aaqed = gal(1)*gal(1)/4d0/3.1415926d0
 
574
 
 
575
      if (btest(mlevel,3)) then
 
576
        write(*,*)' write_leshouche: SCALUP to: ',sscale
 
577
      endif
 
578
      
 
579
      
 
580
      call write_event(lun,pb(0,1),wgt,npart,jpart(1,1),ngroup,
 
581
     &   sscale,aaqcd,aaqed,buff)
 
582
      if(btest(mlevel,1))
 
583
     &   call write_event(6,pb(0,1),wgt,npart,jpart(1,1),ngroup,
 
584
     &   sscale,aaqcd,aaqed,buff)
 
585
 
 
586
      end
 
587
      
 
588
      integer function n_unwgted()
 
589
c************************************************************************
 
590
c     Determines the number of unweighted events which have been written
 
591
c************************************************************************
 
592
      implicit none
 
593
c
 
594
c     Parameter
 
595
c
 
596
      include 'genps.inc'
 
597
      include 'nexternal.inc'
 
598
c
 
599
c     Local
 
600
c
 
601
      integer i
 
602
      double precision xtot, sum
 
603
C     
 
604
C     GLOBAL
 
605
C
 
606
      double precision twgt, maxwgt,swgt(maxevents)
 
607
      integer                             lun, nw, itmin
 
608
      common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
 
609
c-----
 
610
c  Begin Code
 
611
c-----
 
612
 
 
613
c      write(*,*) 'Sorting ',nw
 
614
      if (nw .gt. 1) call dsort(nw,swgt)
 
615
      sum = 0d0
 
616
      do i=1,nw
 
617
         sum=sum+swgt(i)
 
618
      enddo
 
619
      xtot = 0d0
 
620
      i = nw
 
621
      do while (xtot .lt. sum/100d0 .and. i .gt. 2)    !Allow for 1% accuracy
 
622
         xtot = xtot + swgt(i)
 
623
         i=i-1
 
624
      enddo
 
625
      if (i .lt. nw) i = i+1
 
626
c      write(*,*) 'Found ',nw,' events'
 
627
c      write(*,*) 'Integrated weight',sum
 
628
c      write(*,*) 'Maximum wgt',swgt(nw), swgt(i)
 
629
c      write(*,*) 'Average wgt', sum/nw
 
630
c      write(*,*) 'Unweight Efficiency', (sum/nw)/swgt(i)
 
631
      n_unwgted = sum/swgt(i)
 
632
c      write(*,*) 'Number unweighted ',sum/swgt(i), nw
 
633
      if (nw .ge. maxevents) n_unwgted = -sum/swgt(i)
 
634
      end
 
635
 
 
636
 
 
637
      subroutine dsort(n,ra)
 
638
      integer n
 
639
      double precision ra(n),rra
 
640
 
 
641
      l=n/2+1
 
642
      ir=n
 
643
10    continue
 
644
        if(l.gt.1)then
 
645
          l=l-1
 
646
          rra=ra(l)
 
647
        else
 
648
          rra=ra(ir)
 
649
          ra(ir)=ra(1)
 
650
          ir=ir-1
 
651
          if(ir.eq.1)then
 
652
            ra(1)=rra
 
653
            return
 
654
          endif
 
655
        endif
 
656
        i=l
 
657
        j=l+l
 
658
20      if(j.le.ir)then
 
659
          if(j.lt.ir)then
 
660
            if(dabs(ra(j)).lt.dabs(ra(j+1))) j=j+1
 
661
          endif
 
662
          if(dabs(rra).lt.dabs(ra(j)))then
 
663
            ra(i)=ra(j)
 
664
            i=j
 
665
            j=j+j
 
666
          else
 
667
            j=ir+1
 
668
          endif
 
669
        go to 20
 
670
        endif
 
671
        ra(i)=rra
 
672
      go to 10
 
673
      end