← Back to branch summary

~maddevelopers/mg5amcnlo/2.6.6_bug_1813292

« back to all changes in this revision

Viewing changes to Template/NLO/Source/combine_events.f

  • Committer: olivier-mattelaer
  • Date: 2017-05-26 07:48:55 UTC
  • mfrom: (271.1.33 2.5.5)
  • Revision ID: olivier-mattelaer-20170526074855-r463wfxlom110fiu
pass the 2.5.5

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/Source/combine_events.f
1
 
      program test
2
 
c*****************************************************************
3
 
c     tests traversing directories to find all events
4
 
c****************************************************************
5
 
      implicit none
6
 
c
7
 
c     Constants
8
 
c
9
 
      include 'run_config.inc'
10
 
      integer    maxsubprocesses
11
 
      parameter (maxsubprocesses=9999)
12
 
      integer    cmax_events
13
 
      parameter (cmax_events=500000)
14
 
      integer    sfnum
15
 
      parameter (sfnum=17)   !Unit number for scratch file
16
 
      integer    maxexternal
17
 
      parameter (maxexternal=17)
18
 
      integer maxpara
19
 
      parameter (maxpara=1000)
20
 
c     
21
 
c     Local
22
 
c
23
 
      character*30 subname(maxsubprocesses)
24
 
      integer i,j,m,ns,nreq,ievent
25
 
      integer kevent,revent,iarray(cmax_events)
26
 
      double precision sum, xsec, xerr, goal_wgt,xarray(cmax_events)
27
 
      integer i4,r8,record_length
28
 
      integer iseed
29
 
      real xran1
30
 
      double precision wgt,maxwgt
31
 
      double precision p(0:4,maxexternal)
32
 
      integer ic(7,maxexternal),n
33
 
      double precision scale,aqcd,aqed
34
 
      character*20 param(maxpara),value(maxpara)
35
 
      integer npara,nunwgt
36
 
      double precision xtrunc, min_goal,max_goal
37
 
      logical keep(cmax_events),done
38
 
      integer ntry
39
 
      logical gridrun,gridpack
40
 
 
41
 
      character*140 buff
42
 
 
43
 
      data iseed/-1/
44
 
c-----
45
 
c  Begin Code
46
 
c-----
47
 
c
48
 
c     Get requested number of events
49
 
c
50
 
      call load_para(npara,param,value)
51
 
      call get_logical(npara,param,value," gridrun ",gridrun,.false.)
52
 
      call get_logical(npara,param,value," gridpack ",gridpack,.false.)
53
 
      if (gridrun.and.gridpack) then
54
 
         call get_integer(npara,param,value," gevents "  ,nreq  ,2000   )
55
 
      else
56
 
         call get_integer(npara,param,value," nevents "  ,nreq  ,10000   )
57
 
      endif
58
 
 
59
 
c   Get information for the <init> block
60
 
      call setrun
61
 
 
62
 
c      nreq = 10000
63
 
c
64
 
c     Get total cross section
65
 
c
66
 
      xsec = 0d0
67
 
      xerr = 0d0
68
 
      open(unit=15,file='results.dat',status='old',err=21)
69
 
      read(15,*,err=20) xsec,xerr
70
 
      write(*,*) "Results.dat xsec = ",xsec
71
 
 20   close(15)
72
 
 21   if (nreq .gt. 0 .and. xsec .gt. 0) then
73
 
         goal_wgt = xsec/nreq/4d0   !Extra factor of 4 for weighted events
74
 
      else
75
 
         goal_wgt = 0d0    !Write out everything
76
 
      endif
77
 
c
78
 
c     Get list of subprocesses
79
 
c
80
 
      call get_subprocess(subname,ns)
81
 
 
82
 
c
83
 
c     Create scratch file to hold events
84
 
c
85
 
      I4 = 4
86
 
      R8 = 8
87
 
      record_length = 4*I4+maxexternal*I4*7+maxexternal*5*R8+3*R8+
88
 
     &   140
89
 
      open(unit=sfnum,access='direct',file='scratch',err=999,
90
 
     &     recl=record_length)
91
 
c
92
 
c     Loop through subprocesses filling up scratch file with events
93
 
c
94
 
      sum=0d0
95
 
      kevent=0
96
 
      revent=0
97
 
      maxwgt=0d0
98
 
      write(*,*) 'SubProcess/Channel     kept   read   xsec '
99
 
      do i=1,ns
100
 
c         write(*,*) 'Subprocess: ',subname(ns)
101
 
         call read_channels(subname(i),sum,kevent,revent,goal_wgt,maxwgt)
102
 
      enddo 
103
 
c
104
 
c     Get Random order for events
105
 
c
106
 
      do i=1,kevent
107
 
         iarray(i)=i
108
 
         xarray(i)=xran1(iseed)
109
 
      enddo
110
 
      call sortO3(xarray,iarray,kevent)
111
 
c
112
 
c     Write out the events in iarray order
113
 
c
114
 
      open(unit=15,file='events.lhe',status='unknown',err=98)
115
 
      call writebanner(15,kevent,sum,maxwgt,xerr)
116
 
      do i=1,kevent
117
 
         read(sfnum,rec=iarray(i)) wgt,n,
118
 
     &        ((ic(m,j),j=1,maxexternal),m=1,7),ievent,
119
 
     &        ((p(m,j),m=0,4),j=1,maxexternal),scale,aqcd,aqed,
120
 
     &     buff
121
 
         call write_event(15,P,wgt,n,ic,ievent,scale,aqcd,aqed,buff)
122
 
      enddo
123
 
      close(15)
124
 
c
125
 
c     Now select unweighted events.
126
 
c
127
 
      goal_wgt = sum/(nreq*1.03)
128
 
      min_goal = goal_wgt/5d0
129
 
      max_goal = goal_wgt*5d0
130
 
      ntry = 1
131
 
c
132
 
c     Loop to refine guess for goal_wgt while keeping xtrunc<0.01
133
 
c
134
 
      done=.false.
135
 
      do while(.not. done)
136
 
         done=.true.
137
 
         nunwgt=0
138
 
         xtrunc=0d0
139
 
         do i=1,kevent
140
 
            read(sfnum,rec=iarray(i)) wgt,n,
141
 
     &           ((ic(m,j),j=1,maxexternal),m=1,7),ievent,
142
 
     &           ((p(m,j),m=0,4),j=1,maxexternal),scale,aqcd,aqed,
143
 
     &        buff
144
 
            if (wgt .gt. goal_wgt*xran1(iseed)) then
145
 
               keep(i) = .true.
146
 
               nunwgt=nunwgt+1
147
 
               if (wgt .gt. goal_wgt) then
148
 
                  xtrunc=xtrunc+wgt-goal_wgt
149
 
               endif
150
 
            else
151
 
               keep(i)=.false.
152
 
            endif
153
 
         enddo
154
 
         if (xtrunc .gt. 0.01d0*sum) then
155
 
            done=.false.
156
 
            min_goal = max(goal_wgt,min_goal)
157
 
            goal_wgt = goal_wgt*1.3d0            
158
 
            write(*,*) 'Iteration ',ntry, ' too large truncation ',xtrunc/sum,nunwgt
159
 
c            write(*,*) min_goal,goal_wgt,max_goal
160
 
         elseif (nunwgt .lt. nreq) then
161
 
            done=.false.
162
 
            max_goal = min(goal_wgt,max_goal)
163
 
            goal_wgt = goal_wgt*0.95d0
164
 
            write(*,*) 'Iteration ',ntry, ' too few events ',xtrunc/sum,nunwgt
165
 
c            write(*,*) min_goal,goal_wgt,max_goal
166
 
            if (goal_wgt .lt. min_goal) then
167
 
               done=.true.
168
 
               write(*,*) 'Failed to find requested number of unweighted events',nreq,nunwgt
169
 
            endif
170
 
         endif
171
 
         ntry=ntry+1
172
 
         if (ntry .gt. 20) done=.true.
173
 
      enddo
174
 
      if (nunwgt .lt. nreq) then
175
 
         write(*,*) 'Unable to get ',nreq,' events. Writing ',nunwgt
176
 
         nreq = nunwgt
177
 
      else
178
 
         write(*,*) 'Found ',nunwgt,' events writing first ',nreq
179
 
      endif
180
 
      write(*,*) 'Unweighting selected ',nreq, ' events.'
181
 
      write(*,'(a,f5.2,a)') 'Truncated ',xtrunc*100./sum, '% of cross section'
182
 
      open(unit=15,file='unweighted_events.lhe',status='unknown',err=99)
183
 
      call writebanner_u(15,nreq,xsec,xtrunc,xerr)
184
 
      ntry = 0
185
 
      do i=1,kevent
186
 
         if (keep(i) .and. ntry .lt. nreq) then
187
 
            read(sfnum,rec=iarray(i)) wgt,n,
188
 
     &           ((ic(m,j),j=1,maxexternal),m=1,7),ievent,
189
 
     &           ((p(m,j),m=0,4),j=1,maxexternal),scale,aqcd,aqed,
190
 
     $        buff
191
 
            call write_event(15,P,xsec/nreq,n,ic,ievent,scale,aqcd,aqed,buff)
192
 
            ntry=ntry+1
193
 
         endif
194
 
      enddo
195
 
      close(15)
196
 
      close(sfnum)
197
 
      return
198
 
 98   write(*,*) 'Error writing events.dat' 
199
 
      return
200
 
 99   write(*,*) 'Error writing unweighted_events.dat' 
201
 
      return
202
 
 999  write(*,*) 'Error opening scratch file'
203
 
      end
204
 
 
205
 
 
206
 
      subroutine writebanner(lunw,nevent,sum,maxwgt,xerr)
207
 
c**************************************************************************************
208
 
c     Writes out banner information at top of event file
209
 
c**************************************************************************************
210
 
      implicit none
211
 
c
212
 
c     Arguments
213
 
c     
214
 
      integer lunw,nevent
215
 
      double precision sum,maxwgt,xerr
216
 
c
217
 
c     Local
218
 
c
219
 
      integer i,j
220
 
 
221
 
c
222
 
c     Les Houches init block (for the <init> info)
223
 
c
224
 
      integer maxpup
225
 
      parameter(maxpup=100)
226
 
      integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
227
 
      double precision ebmup,xsecup,xerrup,xmaxup
228
 
      common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
229
 
     &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
230
 
     &     xmaxup(maxpup),lprup(maxpup)
231
 
 
232
 
c
233
 
c     Global
234
 
c
235
 
c      double precision etmin(3:nexternal),etamax(3:nexternal)
236
 
c      double precision                    r2min(3:nexternal,3:nexternal)
237
 
c      double precision s_min(nexternal,nexternal)
238
 
c      common/to_cuts/  etmin     ,etamax     , r2min, s_min
239
 
 
240
 
c-----
241
 
c  Begin Code
242
 
c-----
243
 
c
244
 
c     gather the info
245
 
c
246
 
c      call setpara('param_card.dat')
247
 
c      call setcuts
248
 
c
249
 
c     write it out
250
 
c
251
 
c      call write_para(lunw)
252
 
c      write(lunw,'(a70)') '##                                                                    '
253
 
c      write(lunw,'(a70)') '##-------------------                                                 '
254
 
c      write(lunw,'(a70)') '## Run-time options                                                   '
255
 
c      write(lunw,'(a70)') '##-------------------                                                 '
256
 
c      write(lunw,'(a70)') '##                                                                    '
257
 
c     write(lunw,'(a70)') '##********************************************************************'     
258
 
c     write(lunw,'(a70)') '## Standard Cuts                                                     *'
259
 
c     write(lunw,'(a70)') '##********************************************************************'    
260
 
c      write(lunw,'(a13,8i8)')   '## Particle  ',(i,i=3,nexternal)
261
 
c      write(lunw,'(a13,8f8.1)') '## Et       >',(etmin(i),i=3,nexternal)
262
 
c      write(lunw,'(a13,8f8.1)') '## Eta      <',(etamax(i),i=3,nexternal)
263
 
c      do j=3,nexternal-1
264
 
c         write(lunw,'(a,i2,a,8f8.1)') '## d R #',j,'  >',(-0.0,i=3,j),
265
 
c     &        (r2min(i,j),i=j+1,nexternal)
266
 
c         do i=j+1,nexternal
267
 
c            r2min(i,j)=r2min(i,j)**2 !Since r2 returns distance squared
268
 
c         enddo
269
 
c      enddo
270
 
c      do j=3,nexternal-1
271
 
c         write(lunw,'(a,i2,a,8f8.1)') '## s min #',j,'>',
272
 
c     &        (s_min(i,j),i=3,nexternal)
273
 
c      enddo
274
 
c      write(lunw,'(a70)') '#********************************************************************'    
275
 
c
276
 
c     Now write out specific information on the event set
277
 
c
278
 
c
279
 
      write(lunw,'(a)') '<MGGenerationInfo>'
280
 
      write(lunw,'(a30,i10)')   '#  Number of Events        :  ',nevent
281
 
      write(lunw,'(a30,e10.5)') '#  Integrated weight (pb)  :  ',sum
282
 
      write(lunw,'(a30,e10.5)') '#  Max wgt                 :  ',maxwgt
283
 
      write(lunw,'(a30,e10.5)') '#  Average wgt             :  ',sum/nevent
284
 
      write(lunw,'(a)') '</MGGenerationInfo>'
285
 
   
286
 
    
287
 
 
288
 
C   Write out compulsory init info
289
 
      write(lunw,'(a)') '</header>'
290
 
      write(lunw,'(a)') '<init>'
291
 
      write(lunw,90) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
292
 
     $   (pdfsup(i),i=1,2),2,nprup
293
 
      do i=1,nprup
294
 
         write(lunw,91) xsecup(i),xerr*xsecup(i)/sum,maxwgt,lprup(i) ! FACTOR OF nevts for maxwgt and wgt? error?
295
 
      enddo
296
 
      write(lunw,'(a)') '</init>'
297
 
 90   FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)
298
 
 91   FORMAT(3e19.11,i4)
299
 
      end
300
 
 
301
 
 
302
 
      subroutine writebanner_u(lunw,nevent,sum,maxwgt,xerr)
303
 
c**************************************************************************************
304
 
c     Writes out banner information at top of event file
305
 
c**************************************************************************************
306
 
      implicit none
307
 
c
308
 
c     Arguments
309
 
c     
310
 
      integer lunw,nevent
311
 
      double precision sum,maxwgt,xerr
312
 
c
313
 
c     Local
314
 
c
315
 
      integer i,j
316
 
c
317
 
c     Les Houches init block (for the <init> info)
318
 
c
319
 
      integer maxpup
320
 
      parameter(maxpup=100)
321
 
      integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
322
 
      double precision ebmup,xsecup,xerrup,xmaxup
323
 
      common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
324
 
     &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
325
 
     &     xmaxup(maxpup),lprup(maxpup)
326
 
c
327
 
c     Global
328
 
c
329
 
c      double precision etmin(3:nexternal),etamax(3:nexternal)
330
 
c      double precision                    r2min(3:nexternal,3:nexternal)
331
 
c      double precision s_min(nexternal,nexternal)
332
 
c      common/to_cuts/  etmin     ,etamax     , r2min, s_min
333
 
 
334
 
c-----
335
 
c  Begin Code
336
 
c-----
337
 
c
338
 
c     gather the info
339
 
c
340
 
c      call setpara('param_card.dat')
341
 
c      call setcuts
342
 
c
343
 
c     write it out
344
 
c
345
 
c      call write_para(lunw)
346
 
c      write(lunw,'(a70)') '##                                                                    '
347
 
c      write(lunw,'(a70)') '##-------------------                                                 '
348
 
c      write(lunw,'(a70)') '## Run-time options                                                   '
349
 
c      write(lunw,'(a70)') '##-------------------                                                 '
350
 
c      write(lunw,'(a70)') '##                                                                    '
351
 
c     write(lunw,'(a70)') '##********************************************************************'     
352
 
c     write(lunw,'(a70)') '## Standard Cuts                                                     *'
353
 
c     write(lunw,'(a70)') '##********************************************************************'    
354
 
c      write(lunw,'(a13,8i8)')   '## Particle  ',(i,i=3,nexternal)
355
 
c      write(lunw,'(a13,8f8.1)') '## Et       >',(etmin(i),i=3,nexternal)
356
 
c      write(lunw,'(a13,8f8.1)') '## Eta      <',(etamax(i),i=3,nexternal)
357
 
c      do j=3,nexternal-1
358
 
c         write(lunw,'(a,i2,a,8f8.1)') '## d R #',j,'  >',(-0.0,i=3,j),
359
 
c     &        (r2min(i,j),i=j+1,nexternal)
360
 
c         do i=j+1,nexternal
361
 
c            r2min(i,j)=r2min(i,j)**2 !Since r2 returns distance squared
362
 
c         enddo
363
 
c      enddo
364
 
c      do j=3,nexternal-1
365
 
c         write(lunw,'(a,i2,a,8f8.1)') '## s min #',j,'>',
366
 
c     &        (s_min(i,j),i=3,nexternal)
367
 
c      enddo
368
 
c      write(lunw,'(a70)') '##********************************************************************'    
369
 
c
370
 
c     Now write out specific information on the event set
371
 
c
372
 
 
373
 
      write(lunw,'(a)') '<MGGenerationInfo>'
374
 
      write(lunw,'(a30,i10)')   '#  Number of Events        :  ',nevent
375
 
      write(lunw,'(a30,e10.5)') '#  Integrated weight (pb)  :  ',sum
376
 
      write(lunw,'(a30,e10.5)') '#  Truncated wgt (pb)      :  ',maxwgt
377
 
      write(lunw,'(a30,e10.5)') '#  Unit wgt                :  ',sum/nevent
378
 
      write(lunw,'(a)') '</MGGenerationInfo>'
379
 
 
380
 
C   Write out compulsory init info
381
 
      write(lunw,'(a)') '</header>'
382
 
      write(lunw,'(a)') '<init>'
383
 
      write(lunw,90) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
384
 
     $   (pdfsup(i),i=1,2),3,nprup
385
 
      do i=1,nprup
386
 
         write(lunw,91) xsecup(i),xerr*xsecup(i)/sum,sum/nevent,lprup(i) ! FACTOR OF nevts for maxwgt and wgt? error?
387
 
      enddo
388
 
      write(lunw,'(a)') '</init>'
389
 
 90   FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)
390
 
 91   FORMAT(3e19.11,i4)
391
 
 
392
 
      end
393
 
 
394
 
 
395
 
      subroutine read_channels(dir,sum,kevent,revent,goal_wgt,maxwgt)
396
 
c*****************************************************************
397
 
c     tests traversing directories to find all events
398
 
c****************************************************************
399
 
      implicit none
400
 
c
401
 
c     Constants
402
 
c
403
 
      character*(*) symfile
404
 
      parameter (symfile='symfact.dat')
405
 
c
406
 
c     Arguments
407
 
c
408
 
      character*(30) dir
409
 
      integer kevent,revent
410
 
      double precision sum,goal_wgt,maxwgt
411
 
c
412
 
c     Local
413
 
c
414
 
      integer i,j, k, ip
415
 
      double precision xi
416
 
      character*50 dirname,dname,channame
417
 
c-----
418
 
c  Begin Code
419
 
c-----
420
 
       i = index(dir," ")
421
 
      dname = dir(1:i-1)// "/" // symfile
422
 
      open(unit=35, file=dname ,status='old',err=59)
423
 
      do while (.true.)
424
 
         read(35,*,err=99,end=99) xi,j
425
 
         if (j .gt. 0) then
426
 
            if ( (xi-int(xi+.01)) .lt. 1d-5) then
427
 
               k = int(xi+.01)
428
 
               if (k .lt. 10) then
429
 
                  write(dirname,'(a,i1,a)') 'G',k,'/'
430
 
               else if (k .lt. 100) then
431
 
                  write(dirname,'(a,i2,a)') 'G',k,'/'
432
 
               else if (k .lt. 1000) then
433
 
                  write(dirname,'(a,i3,a)') 'G',k,'/'
434
 
               else if (k .lt. 10000) then
435
 
                  write(dirname,'(a,i4,a)') 'G',k,'/'
436
 
               endif
437
 
            else               !Handle B.W.
438
 
               if (xi .lt. 10) then
439
 
                  write(dirname,'(a,f5.3,a,a)') 'G',xi,'/'
440
 
               else if (xi .lt. 100) then
441
 
                  write(dirname,'(a,f6.3,a,a)') 'G',xi,'/'
442
 
               else if (xi .lt. 1000) then
443
 
                  write(dirname,'(a,f7.3,a,a)') 'G',xi,'/'
444
 
               else if (xi .lt. 10000) then
445
 
                  write(dirname,'(a,f8.3,a,a)') 'G',xi,'/'
446
 
               endif
447
 
            endif     
448
 
            ip = index(dirname,'/')
449
 
            channame = dname(1:i-1)// "/" //dirname(1:ip)
450
 
            call read_dir_events(channame(1:i+ip),j,kevent,revent,sum,goal_wgt,maxwgt)
451
 
            write(*,'(a,2i8,e10.3)') channame(1:i+ip),kevent,revent,sum
452
 
         endif
453
 
 98   enddo
454
 
 99   close(35)
455
 
      return
456
 
c
457
 
c     Come here if there isn't a symfact file. Means we will work on
458
 
c     this file alone
459
 
c
460
 
 59   dirname="./"
461
 
      j = 1
462
 
      ip = 2
463
 
      channame = dirname(1:ip)
464
 
      call read_dir_events(channame,j,kevent,revent,sum,goal_wgt,maxwgt)
465
 
      write(*,'(a30,i8,e10.3)') channame(1:i+ip),kevent,sum
466
 
      return
467
 
      end
468
 
 
469
 
      subroutine read_dir_events(channame,nj,kevent,revent,sum,goal_wgt,maxwgt)
470
 
c********************************************************************
471
 
c********************************************************************
472
 
      implicit none
473
 
c
474
 
c     parameters
475
 
c     
476
 
      integer    sfnum
477
 
      parameter (sfnum=17)   !Unit number for scratch file
478
 
      character*(*) scaled_file
479
 
      parameter (scaled_file='events.lhe')
480
 
      integer maxexternal
481
 
      parameter (maxexternal=17)
482
 
      include 'run_config.inc'
483
 
      integer    max_read
484
 
      parameter (max_read = 2000000)
485
 
c
486
 
c     Arguments
487
 
c
488
 
      character*(*) channame
489
 
      integer nj,kevent,revent
490
 
      double precision sum,goal_wgt,maxwgt
491
 
c
492
 
c     Local
493
 
c
494
 
      double precision wgt
495
 
      double precision p(0:4,maxexternal)
496
 
      double precision gsfact
497
 
      real xwgt(max_read),xtot
498
 
      integer i,j,k,m, ic(7,maxexternal),n
499
 
      double precision scale,aqcd,aqed
500
 
      integer ievent,iseed
501
 
      logical done,found
502
 
      character*140 buff
503
 
      character*60 fullname
504
 
c
505
 
c     Les Houches init block (for the <init> info)
506
 
c
507
 
      integer maxpup
508
 
      parameter(maxpup=100)
509
 
      integer idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
510
 
      double precision ebmup,xsecup,xerrup,xmaxup
511
 
      common /heprup/ idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
512
 
     &     idwtup,nprup,xsecup(maxpup),xerrup(maxpup),
513
 
     &     xmaxup(maxpup),lprup(maxpup)
514
 
      data nprup/0/
515
 
c
516
 
c     external
517
 
c
518
 
      real xran1
519
 
c
520
 
c     data
521
 
c
522
 
      data iseed/-1/
523
 
c-----
524
 
c  Begin Code
525
 
c-----     
526
 
      fullname = channame // "gscalefact.dat"
527
 
      gsfact = 1d0
528
 
      open (unit=15,file=fullname,status='old',err=12)
529
 
      read(15,*) gsfact    !Scale factor for grid runs that only use some channels
530
 
 12   close(15)
531
 
      if (gsfact .eq. 0d0) return
532
 
      fullname = channame // scaled_file      
533
 
      open(unit=15,file=fullname, status='old',err=999)
534
 
      done=.false.
535
 
c
536
 
c     Start by initializing all event variables to zero (not really necessary)
537
 
c
538
 
      do j=1,maxexternal
539
 
         do i=1,7
540
 
            ic(i,j)=0
541
 
         enddo
542
 
         do i=0,4
543
 
            p(i,j) = 0d0
544
 
         enddo
545
 
      enddo
546
 
c
547
 
c     Now loop through events
548
 
c
549
 
      do while (.not. done)
550
 
         call read_event(15,P,wgt,n,ic,ievent,scale,aqcd,aqed,buff,done)
551
 
         if (.not. done) then
552
 
            revent = revent+1
553
 
            wgt = wgt*nj*gsfact                 !symmetry factor * grid factor
554
 
            if (wgt .gt. maxwgt) maxwgt=wgt
555
 
            if (wgt .ge. goal_wgt*xran1(iseed)) then
556
 
               kevent=kevent+1
557
 
               if (wgt .lt. goal_wgt) wgt = goal_wgt
558
 
               write(sfnum,rec=kevent) wgt,n,
559
 
     &              ((ic(i,j),j=1,maxexternal),i=1,7),ievent,
560
 
     &              ((p(i,j),i=0,4),j=1,maxexternal),scale,aqcd,aqed,buff
561
 
               sum=sum+wgt
562
 
               found=.false.
563
 
               do i=1,nprup
564
 
                  if(ievent.eq.lprup(i))then
565
 
                     xsecup(i)=xsecup(i)+wgt
566
 
                     found=.true.
567
 
                  endif
568
 
               enddo
569
 
               if(.not.found)then
570
 
                  nprup=nprup+1
571
 
                  lprup(nprup)=ievent
572
 
                  xsecup(nprup)=wgt
573
 
               endif
574
 
            endif
575
 
         endif
576
 
         if (kevent .ge. max_read) then
577
 
            write(*,*) 'Error too many events to read in select_events'
578
 
     $           , kevent
579
 
            write(*,*) 'Reset max_read in Source/select_events.f'
580
 
            stop
581
 
         endif
582
 
      enddo
583
 
 99   close(15)
584
 
 55   format(i3,4e19.11)         
585
 
c      write(*,*) 'Found ',kevent,' events'
586
 
c      write(*,*) 'Integrated weight',sum
587
 
      return
588
 
 999  write(*,*) 'Error opening file ',channame,scaled_file
589
 
 
590
 
      end
591
 
 
592
 
 
593
 
 
594
 
      subroutine get_subprocess(subname,ns)
595
 
c*****************************************************************
596
 
c     tests traversing directories to find all events
597
 
c****************************************************************
598
 
      implicit none
599
 
c
600
 
c     Constants
601
 
c
602
 
      character*(*) plist
603
 
      parameter (plist='subproc.mg')
604
 
      integer    maxsubprocesses
605
 
      parameter (maxsubprocesses=999)
606
 
c
607
 
c     Arguments
608
 
c
609
 
      character*30 subname(maxsubprocesses)
610
 
      integer ns
611
 
c-----
612
 
c  Begin Code
613
 
c-----
614
 
      ns = 1
615
 
      open(unit=15, file=plist,status='old',err=99)
616
 
      do while (.true.)
617
 
         read(15,*,err=999,end=999) subname(ns)
618
 
         ns=ns+1
619
 
      enddo
620
 
 99   subname(ns) = './'
621
 
      write(*,*) "Did not find ", plist
622
 
      return
623
 
 999  ns = ns-1
624
 
      write(*,*) "Found ", ns," subprocesses"
625
 
      close(15)
626
 
      end
627
 
 
628
 
 
629
 
      function xran1(idum)
630
 
      dimension r(97)
631
 
      parameter (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6)
632
 
      parameter (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6)
633
 
      parameter (m3=243000,ia3=4561,ic3=51349)
634
 
      data iff /0/
635
 
      save r, ix1,ix2,ix3
636
 
      if (idum.lt.0.or.iff.eq.0) then
637
 
        iff=1
638
 
        ix1=mod(ic1-idum,m1)
639
 
        ix1=mod(ia1*ix1+ic1,m1)
640
 
        ix2=mod(ix1,m2)
641
 
        ix1=mod(ia1*ix1+ic1,m1)
642
 
        ix3=mod(ix1,m3)
643
 
        do 11 j=1,97
644
 
          ix1=mod(ia1*ix1+ic1,m1)
645
 
          ix2=mod(ia2*ix2+ic2,m2)
646
 
          r(j)=(float(ix1)+float(ix2)*rm2)*rm1
647
 
11      continue
648
 
        idum=1
649
 
      endif
650
 
      ix1=mod(ia1*ix1+ic1,m1)
651
 
      ix2=mod(ia2*ix2+ic2,m2)
652
 
      ix3=mod(ia3*ix3+ic3,m3)
653
 
      j=1+(97*ix3)/m3
654
 
      if(j.gt.97.or.j.lt.1)then
655
 
         write(*,*) 'j is bad in ran1.f',j, 97d0*ix3/m3
656
 
         STOP
657
 
      endif
658
 
      xran1=r(j)
659
 
      r(j)=(float(ix1)+float(ix2)*rm2)*rm1
660
 
      return
661
 
      end
662
 
 
663
 
 
664
 
      subroutine sort2(array,aux1,n)
665
 
      implicit none
666
 
! Arguments
667
 
      integer n
668
 
      integer aux1(n)
669
 
      double precision array(n)
670
 
!  Local Variables
671
 
      integer i,k
672
 
      double precision temp
673
 
      logical done
674
 
 
675
 
!-----------
676
 
! Begin Code
677
 
!-----------
678
 
      do i=n-1,1,-1
679
 
         done = .true.
680
 
         do k=1,i
681
 
            if (array(k) .lt. array(k+1)) then
682
 
               temp = array(k)
683
 
               array(k) = array(k+1)
684
 
               array(k+1) = temp
685
 
               temp = aux1(k)
686
 
               aux1(k) = aux1(k+1)
687
 
               aux1(k+1) = temp
688
 
               done = .false.
689
 
            end if
690
 
         end do
691
 
         if (done) return
692
 
      end do
693
 
      end 
694
 
 
695
 
      subroutine sortO3(array,aux1,n)
696
 
 
697
 
c   O-Sort Version 3, Sorting routine by Erik Oosterwal
698
 
c   http://www.geocities.com/oosterwal/computer/sortroutines.html
699
 
 
700
 
      implicit none
701
 
 
702
 
! Arguments
703
 
      integer n
704
 
      integer aux1(n)
705
 
      double precision array(n)
706
 
!  Local Variables
707
 
      integer step,i,itemp
708
 
      double precision SngPhi,SngFib
709
 
 
710
 
      SngPhi = 0.78             ! Define phi value
711
 
      SngFib = n * SngPhi       ! Set initial real step size
712
 
      step = int(SngFib)        ! set initial integer step size
713
 
 
714
 
      do while (step > 0)
715
 
        do i = 1,n-step         ! Set the range of the lower search cells
716
 
          if (array(aux1(i))<array(aux1(i+step))) then ! Compare cells
717
 
            itemp = aux1(i)     !                       \
718
 
            aux1(i) = aux1(i+step) !                     | Swap cells
719
 
            aux1(i+step) = itemp !                      /
720
 
          end if
721
 
        enddo
722
 
        
723
 
        SngFib = SngFib * SngPhi ! Decrease the Real step size
724
 
        Step = Int(SngFib)      ! Set the integer step value
725
 
 
726
 
      enddo
727
 
 
728
 
      end