← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_129738/PROC_129738/Source/gen_ximprove.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_129738/PROC_129738/Source/gen_ximprove.f
 
1
      program gen_ximprove
 
2
c*****************************************************************************
 
3
c     Program to combine results from all of the different sub amplitudes 
 
4
c     and given total  cross section and error.
 
5
c*****************************************************************************
 
6
      implicit none
 
7
c
 
8
c     Constants
 
9
c
 
10
      character*(*) rfile
 
11
      parameter (rfile='results.dat')
 
12
      character*(*) symfile
 
13
      parameter (symfile='symfact.dat')
 
14
 
 
15
      include 'run_config.inc'
 
16
      include 'maxparticles.inc'
 
17
      include 'maxconfigs.inc'
 
18
c
 
19
c     global
 
20
c
 
21
      integer max_np,min_iter
 
22
      common/max_np/max_np,min_iter
 
23
 
 
24
c
 
25
c     local
 
26
c
 
27
      double precision xsec(lmaxconfigs), xerr(lmaxconfigs)
 
28
      double precision xerru(lmaxconfigs),xerrc(lmaxconfigs)
 
29
      double precision xmax(lmaxconfigs), eff(lmaxconfigs)
 
30
      double precision xlum(lmaxconfigs)
 
31
      double precision ysec, yerr, yeff, ymax
 
32
      double precision tsec, terr, teff, tmax, xi
 
33
      integer nw(lmaxconfigs), nevents(lmaxconfigs), maxit
 
34
      integer nunwgt(lmaxconfigs)
 
35
      character*80 fname, gname(lmaxconfigs)
 
36
      integer i,j,k,l,n,ipp
 
37
      double precision xtot,errtot,err_goal
 
38
      double precision errtotu,errtotc
 
39
      integer mfact(lmaxconfigs)
 
40
      logical parallel, gen_events
 
41
      character*20 param(maxpara),value(maxpara)
 
42
      integer npara, nreq, ngran, nhel_refine
 
43
      integer ij, kl, ioffset
 
44
      integer*8 iseed     !tjs 20/6/2012 to avoid integer overflow
 
45
      logical Gridpack,gridrun
 
46
      logical split_channels
 
47
      common /to_split/split_channels
 
48
      integer ncode,npos
 
49
      character*20 formstr
 
50
      logical file_exists
 
51
      character*30 filename
 
52
 
 
53
c-----
 
54
c  Begin Code
 
55
c-----
 
56
      call load_para(npara,param,value)
 
57
      call get_logical(npara,param,value," gridpack ",gridpack,.false.)
 
58
      call get_integer(npara,param,value," nhel ",nhel_refine,0)
 
59
c     If different card options set for nhel_refine and nhel_survey:
 
60
      call get_integer(npara,param,value," nhel_refine ",nhel_refine,
 
61
     $     1*nhel_refine)
 
62
      if (.not. Gridpack) then
 
63
         write(*,'(a,a)')'Enter fractional accuracy (<1)',
 
64
     &        ', or number events (>1), max processes per job',
 
65
     &        ', and whether to split channels (T/F)'
 
66
         read(5,*) err_goal, max_np, split_channels
 
67
         min_iter=3
 
68
         parallel = .false.
 
69
         if (err_goal .lt. 1) then
 
70
            write(*,'(a,f8.2,a)') 'Running for accuracy of ',
 
71
     $           err_goal*100,'%'
 
72
            gen_events=.false.         
 
73
         elseif (err_goal .gt. 1) then
 
74
            write(*,'(a,f9.0,a)') 'Generating ',err_goal,
 
75
     &           ' unweighted events.'
 
76
            gen_events=.true.         
 
77
            err_goal = err_goal * 1.2 !Extra factor to ensure works
 
78
         else
 
79
            write(*,*) 'Error, need non_zero goal'
 
80
            stop
 
81
         endif
 
82
      else
 
83
         gen_events=.true.
 
84
         split_channels=.false.
 
85
c        Allow all the way down to a single iteration for gridruns
 
86
         min_iter=1
 
87
         call get_integer(npara,param,value," gevents "  ,nreq  ,2000   )
 
88
         err_goal = 1.2*nreq ! extra factor to ensure works
 
89
         call get_int8(npara,param,value," gseed "  ,iseed  ,4321   )
 
90
         call get_integer(npara,param,value," ngran "  ,ngran  , -1)
 
91
         if (ngran.eq.-1) ngran = 1
 
92
         write(*,*) "Running on Grid to generate ",nreq," events"
 
93
         write(*,*) "   with granularity equal to ",ngran
 
94
c
 
95
c     TJS 3/13/2008
 
96
c     Modified to allow for more sequences
 
97
c     iseed can be between 0 and 30081*30081
 
98
c     before patern repeats
 
99
c     JA 11/2/2011: Check for ioffset, as in ntuple (ranmar.f)
 
100
c     TJS  20/6/2012 changed mod value to 30081 to avoid duplicate sequences
 
101
c
 
102
         call get_offset(ioffset)
 
103
         iseed = iseed * 31300       
 
104
         ij=1802 + mod(iseed,30081)      
 
105
         kl=9373 + (iseed/30081) + ioffset
 
106
         write(*,'($a,i6,a3,i6)') 'Using random seed offset: ',ioffset
 
107
         write(*,*) ' with seed', iseed
 
108
         do while (ij .gt. 31328)
 
109
            ij = ij - 31328
 
110
         enddo
 
111
         do while (kl .gt. 30081)
 
112
            kl = kl - 30081
 
113
         enddo
 
114
         write(*,*) "Using random seeds",ij,kl
 
115
        call rmarin(ij,kl)
 
116
      endif
 
117
      open(unit=15,file=symfile,status='old',err=999)
 
118
      errtot=0d0
 
119
      errtotu=0d0
 
120
      errtotc=0d0
 
121
      xtot = 0d0
 
122
      i = 0
 
123
c     ncode is number of digits needed for the bw coding
 
124
      ncode=int(dlog10(3d0)*(max_particles-3))+1
 
125
      do while (.true.)
 
126
         read(15,*,err=99,end=99) xi,j
 
127
         if (j .gt. 0) then
 
128
            i = i+1
 
129
            k = int(xi*(1+10**(-ncode)))
 
130
            npos=int(dlog10(dble(k)))+1
 
131
            if ( (xi-k) .eq. 0) then
 
132
c              Write with correct number of digits
 
133
               write(formstr,'(a,i1,a)') '(a,i',npos,',a,a)'
 
134
               write(fname, formstr) 'G',k,'/',rfile
 
135
            else               !Handle B.W.
 
136
c              Write with correct number of digits
 
137
               write(formstr,'(a,i1,a,i1,a)') '(a,f',npos+ncode+1,
 
138
     $                 '.',ncode,',a,a)'
 
139
               write(fname, formstr) 'G',xi,'/',rfile
 
140
            endif
 
141
c              write(*,*) 'log name ',fname
 
142
         endif
 
143
         if (j .gt. 0) then
 
144
            gname(i)=fname
 
145
            nevents(i)=0d0
 
146
            xsec(i)=0d0
 
147
            xerr(i)=0d0
 
148
            nw(i)  =0d0
 
149
            mfact(i)=j
 
150
 
 
151
c
 
152
c     Read in integration data from run
 
153
c
 
154
            open(unit=25,file=fname,status='old',err=95)
 
155
            read(25,*,err=94,end=94) xsec(i),xerru(i),xerrc(i),nevents(i),nw(i),maxit,
 
156
     &           nunwgt(i),xlum(i)
 
157
            if (xsec(i) .eq. 0d0) xlum(i)=1d99     !zero cross section
 
158
            xlum(i) = xlum(i)/1000   !convert to fb^-1 
 
159
            xerr(i)=sqrt(xerru(i)**2+xerrc(i)**2)
 
160
            if (.false.) then
 
161
c            maxit = 2
 
162
               tmax = -1d0
 
163
               terr = 0d0
 
164
               teff = 0d0
 
165
               tsec = 0d0
 
166
               do k=1,maxit
 
167
                  read(25,*,err=92) l,ysec,yerr,yeff,ymax
 
168
                  if (k .gt. 1) tmax = max(tmax,ymax)
 
169
                  tsec = tsec + ysec
 
170
                  terr = terr +yerr**2
 
171
                  teff = teff + yeff
 
172
               enddo
 
173
 92            maxit = k-1      !In case of error reading file
 
174
               xsec(i)=tsec/maxit
 
175
               xerr(i)=sqrt(terr)/maxit
 
176
               xmax(i)=tmax/xsec(i)
 
177
            endif
 
178
c            tmax
 
179
            xmax(i) = -1d0
 
180
            xsec(i) = xsec(i)*mfact(i)
 
181
            xerr(i) = xerr(i)*mfact(i)
 
182
            xerru(i) = xerru(i)*mfact(i)
 
183
            xerrc(i) = xerrc(i)*mfact(i)
 
184
            xlum(i) = xlum(i)/mfact(i)
 
185
            xtot = xtot+ xsec(i)
 
186
            eff(i)= xerr(i)*sqrt(real(nevents(i)))/(xsec(i)+1d-99)
 
187
            errtotu = errtotu+(xerru(i))**2
 
188
            errtotc = errtotc+(xerrc(i))
 
189
c            xtot = xtot+ xsec(i)*mfact(i)
 
190
c            eff(i)= xerr(i)*sqrt(real(nevents(i)))/xsec(i)
 
191
c            errtot = errtot+(mfact(i)*xerr(i))**2
 
192
            goto 95
 
193
 94         continue
 
194
c        There was an error reading an existing results.dat file
 
195
c        Stop generation with error message
 
196
            filename='../../error'
 
197
            INQUIRE(FILE="../../RunWeb", EXIST=file_exists)
 
198
            if(.not.file_exists) filename = '../' // filename
 
199
            open(unit=26,file=filename,status='unknown')
 
200
            write(26,*) 'Bad results.dat file for channel ',xi
 
201
 95         close(25)
 
202
c            write(*,*) i,maxit,xsec(i), eff(i)
 
203
         else
 
204
c            i=i-1   !This is for case w/ B.W. and optimization
 
205
         endif
 
206
      enddo
 
207
 99   close(15)
 
208
      errtot=sqrt(errtotc**2+errtotu)
 
209
      if ( .not. gen_events) then
 
210
         call write_bash(xsec,xerru,xerrc,xtot,mfact,err_goal,
 
211
     $        i,nevents,gname,nhel_refine)
 
212
      else
 
213
         open(unit=25,file='../results.dat',status='old',err=199)
 
214
         read(25,*) xtot
 
215
         write(*,'(a,e12.3)') 'Reading total xsection ',xtot
 
216
 199     close(25)
 
217
         if (gridpack) then
 
218
            call write_gen_grid(err_goal,dble(ngran),i,nevents,gname,
 
219
     $           xlum,xtot,mfact,xsec,nhel_refine)
 
220
         else
 
221
            call write_gen(err_goal,i,nevents,gname,xlum,xtot,mfact,
 
222
     $           xsec,xerr,nhel_refine)
 
223
         endif
 
224
      endif
 
225
      stop
 
226
 999  write(*,*) 'error'
 
227
      end
 
228
 
 
229
 
 
230
      subroutine write_bash(xsec,xerru,xerrc,xtot,
 
231
     $     mfact,err_goal,ng,jpoints,gn,nhel_refine)
 
232
c*****************************************************************************
 
233
c     Writes out bash commands for running each channel as needed.
 
234
c*****************************************************************************
 
235
      implicit none
 
236
c
 
237
c     Constants
 
238
c
 
239
      include 'run_config.inc'
 
240
      include 'maxconfigs.inc'
 
241
 
 
242
c      integer    max_np
 
243
c      parameter (max_np = 30)
 
244
c
 
245
c     global
 
246
c
 
247
      integer max_np,min_iter
 
248
      common/max_np/max_np,min_iter
 
249
c
 
250
c     Arguments
 
251
c
 
252
      double precision xsec(lmaxconfigs), xerru(lmaxconfigs),xerrc(lmaxconfigs)
 
253
      double precision err_goal,xtot
 
254
      integer mfact(lmaxconfigs),jpoints(lmaxconfigs),nhel_refine
 
255
      integer ng
 
256
      character*(80) gn(lmaxconfigs)
 
257
c
 
258
c     Local
 
259
c
 
260
      integer i,j,k, io(lmaxconfigs), npoints, ip, np
 
261
      double precision xt(lmaxconfigs),elimit
 
262
      double precision yerr,ysec,rerr
 
263
      logical fopened
 
264
c-----
 
265
c  Begin Code
 
266
c-----
 
267
      fopened = .false.
 
268
      k=0
 
269
      do j=1,ng
 
270
         if (mfact(j) .gt. 0) k=k+1
 
271
         io(j) = j
 
272
         xt(j)= sqrt((xerru(j)+xerrc(j)**2)*mfact(j))     !sort by error
 
273
      enddo
 
274
c
 
275
c     Let's redetermine err_goal based on luminosity
 
276
c
 
277
      write(*,*) 'Cross section pb',xtot
 
278
      write(*,*) 'Desired Goal',err_goal      
 
279
      write(*,*) 'Total Error',err_goal
 
280
c      elimit = err_goal*xtot/sqrt(real(k)) !Equal contributions from each
 
281
      elimit = err_goal*xtot/real(k) !Equal contributions from each
 
282
 
 
283
      call sort2(xt,io,ng)
 
284
      k=1
 
285
      xt(ng+1) = 0
 
286
      do while( xt(k) .gt. abs(elimit))  !abs is just in case elimit<0 by mistake
 
287
         k=k+1
 
288
      enddo
 
289
      k=k-1
 
290
      rerr=0d0
 
291
      do j=k+1,ng
 
292
c         rerr = rerr+xt(j)**2
 
293
         rerr = rerr+xt(j)
 
294
      enddo
 
295
      rerr=rerr**2
 
296
c      write(*,*) 'Number of diagrams to fix',k
 
297
c
 
298
c     Now readjust because most don't contribute
 
299
c
 
300
      elimit = sqrt((err_goal*xtot)**2 - rerr)/sqrt(real(k))
 
301
 
 
302
      
 
303
      np = max_np
 
304
      do i=1,k
 
305
 
 
306
c         yerr = xerr(io(i))*mfact(io(i))
 
307
         yerr = xt(i)
 
308
c         write(*,*) i,xt(i),elimit
 
309
         if (yerr .gt. elimit) then
 
310
 
 
311
         ysec = xsec(io(i)) + yerr
 
312
         npoints=(0.2d0)*jpoints(io(i))*(yerr/elimit)**2
 
313
         npoints = max(npoints,min_events)
 
314
         npoints = min(npoints,max_events)
 
315
c         np = np + 3*npoints
 
316
         np = np +1
 
317
         if (np .gt. max_np) then
 
318
            if (fopened) then
 
319
               call close_bash_file(26)
 
320
            endif
 
321
            fopened=.true.
 
322
            call open_bash_file(26)
 
323
c            np = 3*npoints
 
324
            np = 1
 
325
         endif
 
326
 
 
327
         ip = index(gn(io(i)),'/')
 
328
         write(*,*) 'Channel ',gn(io(i))(2:ip-1),
 
329
     $        yerr, jpoints(io(i)),npoints
 
330
 
 
331
         ip = index(gn(io(i)),'/')
 
332
         write(26,'(2a)') 'j=',gn(io(i))(1:ip-1)
 
333
c
 
334
c     Determine estimates for getting the desired accuracy
 
335
c
 
336
 
 
337
c
 
338
c     Now write the commands
 
339
c      
 
340
         write(26,20) 'if [[ ! -e $j ]]; then'
 
341
         write(26,25) 'mkdir $j'
 
342
         write(26,20) 'fi'
 
343
         write(26,20) 'cd $j'
 
344
         write(26,20) 'rm -f $k'
 
345
c         write(26,20) 'rm -f moffset.dat'
 
346
 
 
347
         write(26,'(5x,a,3i8,a)') 'echo "',npoints,max_iter,min_iter,
 
348
     $        '" >& input_sg.txt' 
 
349
         write(26,'(5x,a,f8.3,a)') 'echo "',max(elimit/ysec,0.001d0),
 
350
     $        '" >> input_sg.txt'
 
351
         write(26,'(5x,a)') 'echo "2" >> input_sg.txt'  !Grid
 
352
         write(26,'(5x,a)') 'echo "1" >> input_sg.txt'  !Suppress
 
353
         write(26,'(5x,a,i4,a)') 'echo "',nhel_refine,
 
354
     &        '"  >> input_sg.txt' !Helicity 
 
355
         write(26,'(5x,3a)')'echo "',gn(io(i))(2:ip-1),
 
356
     $        '" >>input_sg.txt'
 
357
         write(26,20) '../madevent >> $k <input_sg.txt'
 
358
         write(26,20) 'rm ftn25 ftn26'
 
359
         write(26,20) 'cat $k >> log.txt'
 
360
         write(26,20) 'cd ../'
 
361
      endif
 
362
      enddo  !Loop over diagrams
 
363
      if (fopened) then
 
364
         call close_bash_file(26)
 
365
      endif
 
366
      fopened=.false.
 
367
 15   format(a)
 
368
 20   format(5x,a)
 
369
 25   format(10x,a)
 
370
 999  close(26)
 
371
      end
 
372
 
 
373
 
 
374
      subroutine open_bash_file(lun)
 
375
c***********************************************************************
 
376
c     Opens bash file for looping including standard header info
 
377
c     which can be used with pbs, or stand alone
 
378
c***********************************************************************
 
379
      implicit none
 
380
c
 
381
c     Constants
 
382
c
 
383
      include 'run_config.inc'
 
384
c
 
385
c     Arguments
 
386
c
 
387
      integer lun
 
388
c
 
389
c     local
 
390
c
 
391
      character*30 fname
 
392
      integer ic, npos
 
393
      character*10 formstr
 
394
 
 
395
      data ic/0/
 
396
c-----
 
397
c  Begin Code
 
398
c-----
 
399
      ic=ic+1
 
400
      fname='ajob'
 
401
c     Write ic with correct number of digits
 
402
      npos=int(dlog10(dble(ic)))+1
 
403
      write(formstr,'(a,i1,a)') '(I',npos,')'
 
404
      write(fname(5:(5+npos-1)),formstr) ic
 
405
 
 
406
      write(*,*) 'Opening file ',fname
 
407
      open (unit=26, file = fname, status='unknown')
 
408
      write(26,15) '#!/bin/bash'
 
409
c      write(26,15) '#PBS -q ' // PBS_QUE
 
410
c      write(26,15) '#PBS -o /dev/null'
 
411
c      write(26,15) '#PBS -e /dev/null'
 
412
c      write(26,15) 'if [[ "$PBS_O_WORKDIR" != "" ]]; then' 
 
413
c      write(26,15) '    cd $PBS_O_WORKDIR'
 
414
c      write(26,15) 'fi'
 
415
      write(26,15) 'k=run1_app.log'
 
416
      write(lun,15) 'script=' // fname
 
417
c      write(lun,15) 'rm -f wait.$script >& /dev/null'
 
418
c      write(lun,15) 'touch run.$script'
 
419
 15   format(a)
 
420
      end
 
421
 
 
422
      subroutine close_bash_file(lun)
 
423
c***********************************************************************
 
424
c     Closes bash file
 
425
c***********************************************************************
 
426
      implicit none
 
427
c
 
428
c     Constants
 
429
c
 
430
c
 
431
c     Arguments
 
432
c
 
433
      integer lun
 
434
c
 
435
c     local
 
436
c
 
437
      character*30 fname
 
438
      integer ic
 
439
 
 
440
      data ic/0/
 
441
c-----
 
442
c  Begin Code
 
443
c-----
 
444
 
 
445
c      write(lun,'(a)') ')'
 
446
c
 
447
c     Now write the commands
 
448
c      
 
449
c      write(lun,20) 'j=G$i'
 
450
c      write(lun,20) 'if (! -e $j) then'
 
451
c      write(lun,25) 'mkdir $j'
 
452
c      write(lun,20) 'endif'
 
453
c      write(lun,20) 'cd $j'
 
454
c      write(lun,20) 'rm -f ftn25 ftn99'
 
455
c      write(lun,20) 'rm -f $k'
 
456
c      write(lun,20) 'cat ../input_app.txt >& input_app.txt'
 
457
c      write(lun,20) 'echo $i >> input_app.txt'
 
458
c      if (.false.) then
 
459
c         write(lun,20) 'cp ../../public.sh .'
 
460
c         write(lun,20) 'qsub -N $1$i public.sh >> ../../running_jobs'
 
461
c      else
 
462
c         write(lun,20) '../madevent > $k <input_app.txt'
 
463
c         write(lun,20) 'rm -f ftn25 ftn99'
 
464
c         write(lun,20) 'cp $k log.txt'
 
465
c      endif
 
466
c      write(lun,20) 'cd ../'
 
467
c      write(lun,15) 'end'
 
468
c      write(lun,15) 'rm -f run.$script >&/dev/null'
 
469
c      write(lun,15) 'touch done.$script >&/dev/null'
 
470
 15   format(a)
 
471
 20   format(5x,a)
 
472
 25   format(10x,a)
 
473
      close(lun)
 
474
      end
 
475
 
 
476
 
 
477
 
 
478
      subroutine write_gen(goal_lum,ng,jpoints,gn,xlum,xtot,mfact,xsec,
 
479
     $     xerr,nhel_refine)
 
480
c*****************************************************************************
 
481
c     Writes out scripts for achieving unweighted event goals
 
482
c*****************************************************************************
 
483
      implicit none
 
484
c
 
485
c     Constants
 
486
c
 
487
      include 'run_config.inc'
 
488
      include 'maxconfigs.inc'
 
489
c
 
490
c     global
 
491
c
 
492
      integer max_np,min_iter
 
493
      common/max_np/max_np,min_iter
 
494
c      integer    max_np     !now set in run_config.inc
 
495
c      parameter (max_np = 5)  !number of channels/job
 
496
 
 
497
c
 
498
c     Arguments
 
499
c
 
500
      double precision goal_lum, xlum(lmaxconfigs), xsec(lmaxconfigs),xtot
 
501
      double precision xerr(lmaxconfigs)
 
502
      integer jpoints(lmaxconfigs), mfact(lmaxconfigs)
 
503
      integer ng, np, nhel_refine
 
504
      character*(80) gn(lmaxconfigs)
 
505
c
 
506
c     Local
 
507
c
 
508
      integer i,j,k,kk, io(lmaxconfigs), npoints, ip, nfiles,ifile,npfile
 
509
      double precision xt(lmaxconfigs+1),elimit
 
510
      double precision yerr,ysec,rerr
 
511
      logical fopened
 
512
      character*26 cjobs
 
513
      integer mjobs,ijob,jc
 
514
      character*150 fname
 
515
 
 
516
      logical split_channels
 
517
      common /to_split/split_channels
 
518
 
 
519
      data cjobs/"abcdefghijklmnopqrstuvwxyz"/
 
520
 
 
521
c-----
 
522
c  Begin Code
 
523
c-----
 
524
      fopened=.false.
 
525
      write(*,*) 'Working on creating ', goal_lum, ' events.'
 
526
      goal_lum = goal_lum/(xtot*1000) !Goal luminosity in fb^-1
 
527
      write(*,*) 'Effective Luminosity', goal_lum, ' fb^-1.'
 
528
      k=0
 
529
      do j=1,ng
 
530
         io(j) = j
 
531
         xt(j)= goal_lum/(xlum(j)+1d-99)       !sort by events_needed/have.
 
532
         write(*,*) j,xlum(j),xt(j)
 
533
      enddo
 
534
c      write(*,*) 'Number of channels',ng,k
 
535
 
 
536
c    Reset multijob.dat for all channels
 
537
      do j=1,ng
 
538
        jc = index(gn(j),"/")
 
539
        fname = gn(j)(1:jc)// "multijob.dat"
 
540
        write(*,*) 'Resetting ' // fname
 
541
        open(unit=15,file=fname,status="unknown",err=10)
 
542
        write(15,*) 0
 
543
 10     close(15)
 
544
      enddo
 
545
c
 
546
c     Let's redetermine err_goal based on luminosity
 
547
c
 
548
      elimit = 1d0
 
549
      call sort2(xt,io,ng)
 
550
      k=1
 
551
      xt(ng+1) = 0
 
552
      do while( xt(k) .gt. abs(elimit)) !elimit should be >0
 
553
         write(*,*) 'Improving ',k,gn(io(k)),xt(k)
 
554
         k=k+1
 
555
      enddo
 
556
      kk=k
 
557
c     Check error for the rest of the channels - rerun if 
 
558
c     bigger than channel xsec and bigger than 1% of largest channel
 
559
      do while( kk .le. ng)
 
560
         if (xerr(io(kk)).gt.max(xsec(io(kk)),0.01*xsec(io(1)))) then
 
561
            write(*,*) 'Improving for error ',kk,gn(io(kk)),xt(kk),xsec(io(kk)),xerr(io(kk))
 
562
            io(k)=io(kk)
 
563
            xt(k)=xt(kk)
 
564
            k=k+1
 
565
         endif
 
566
         kk=kk+1
 
567
      enddo
 
568
      k=k-1
 
569
      write(*,*) 'Number of diagrams to fix',k
 
570
c
 
571
c     Now readjust because most don't contribute
 
572
c
 
573
 
 
574
c      np = max_np
 
575
 
 
576
c
 
577
c     Want to write channels so that heaviest one (with largest error)
 
578
c     gets grouped with least heavy channels. Complicated ordering for this
 
579
c     follows. np is the present process number.
 
580
c
 
581
      nfiles = k/max_np
 
582
      if(mod(k,max_np).gt.0) nfiles=nfiles+1
 
583
      ifile  = 0
 
584
      npfile = 0
 
585
      np = 1
 
586
      
 
587
 
 
588
      do i=1,k
 
589
         yerr = xt(np)
 
590
         npoints=0.2*jpoints(io(np))*(yerr/elimit)
 
591
         npoints = max(npoints,min_events)
 
592
         npoints = min(npoints,max_events)
 
593
 
 
594
         npfile=npfile+1
 
595
c         np = nfiles*npfile+1-ifile   !Fancy order for combining channels removed 12/6/2010 by tjs
 
596
         np = i
 
597
c
 
598
c     tjs 12/5/2010
 
599
c     Add loop to allow for multiple jobs on a single channel
 
600
c
 
601
         mjobs = (goal_lum*xsec(io(np))*1000 / MaxEventsPerJob + 0.9)
 
602
c         write(*,*) "Working on Channel ",i,io(np),xt(np), goal_lum*xsec(io(np))*1000 /maxeventsperjob
 
603
         if (mjobs .gt. 26)  then
 
604
            write(*,*) 'Error in gen_ximprove.f, too many events requested ',mjobs*maxeventsperjob
 
605
            mjobs=26
 
606
         endif
 
607
         if (mjobs .lt. 1 .or. .not. split_channels)  mjobs=1
 
608
c
 
609
c        write multijob.dat file for combine_runs.f 
 
610
c
 
611
         jc = index(gn(io(np)),"/")
 
612
         fname = gn(io(np))(1:jc)// "multijob.dat"
 
613
c            write(*,*) "Writing file ", fname
 
614
         open(unit=15,file=fname,status="unknown",err=11)
 
615
         if (mjobs .gt. 1) then
 
616
           write(15,*) mjobs
 
617
         else
 
618
           write(15,*) 0
 
619
         endif
 
620
 11      close(15)
 
621
         do ijob = 1, mjobs
 
622
c---
 
623
c tjs
 
624
c---
 
625
         if (npfile .gt. max_np .or. ifile.eq.0 .or. mjobs .gt. 1) then
 
626
            if (fopened) then
 
627
               call close_bash_file(26)
 
628
            endif
 
629
            fopened=.true.
 
630
            call open_bash_file(26)
 
631
            ifile=ifile+1
 
632
            npfile=1
 
633
c            if (ijob .eq. 1)  np = ifile !Only increment once / source channel
 
634
         endif
 
635
         ip = index(gn(io(np)),'/')
 
636
         write(*,*) 'Channel ',gn(io(np))(2:ip-1),
 
637
     $        yerr, jpoints(io(np)),npoints
 
638
 
 
639
         ip = index(gn(io(np)),'/')
 
640
         if (mjobs .gt. 1) then
 
641
            write(26,'(3a)') 'j=',gn(io(np))(1:ip-1),cjobs(ijob:ijob)
 
642
         else
 
643
            write(26,'(3a)') 'j=',gn(io(np))(1:ip-1)
 
644
         endif
 
645
c
 
646
c     Now write the commands
 
647
c      
 
648
         write(26,20) 'if [[ ! -e $j ]]; then'
 
649
         write(26,25) 'mkdir $j'
 
650
         write(26,20) 'fi'
 
651
         write(26,20) 'cd $j'
 
652
         write(26,20) 'rm -f $k'
 
653
         write(26,20) 'rm -f moffset.dat >& /dev/null'
 
654
         write(26,*) '     echo ',ijob, ' > moffset.dat'
 
655
 
 
656
c
 
657
c     
 
658
c
 
659
c
 
660
c     Now I'll add a check to make sure the grid has been
 
661
c     adjusted  (ftn99 or ftn25 exist)
 
662
c
 
663
         write(26,20) 'if  [[ -e ftn26 ]]; then'
 
664
         write(26,25) 'cp ftn26 ftn25'
 
665
         write(26,20) 'fi'
 
666
 
 
667
         write(26,20) 'if [[ ! -e ftn25 ]]; then'
 
668
 
 
669
 
 
670
         write(26,'(9x,a,3i8,a)') 'echo "',npoints,max_iter,min_iter,
 
671
     $        '" >& input_sg.txt' 
 
672
c
 
673
c     tjs 8/7/2007-JA 8/17/11 Allow stop when have enough luminocity
 
674
c
 
675
         write(*,*) "Cross section",i,io(np),xsec(io(np)),mfact(io(np))
 
676
         write(26,'(9x,a,e13.5,a)') 'echo "',-goal_lum*1000/mjobs,
 
677
     $        '" >> input_sg.txt'                       !Luminocity
 
678
         write(26,'(9x,a)') 'echo "2" >> input_sg.txt'  !Grid Adjustment
 
679
         write(26,'(9x,a)') 'echo "1" >> input_sg.txt'  !Suppression
 
680
         write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
 
681
     &        ' " >> input_sg.txt' !Helicity 0=exact
 
682
         write(26,'(9x,3a)')'echo "',gn(io(np))(2:ip-1),
 
683
     $        '" >>input_sg.txt'
 
684
         write(26,25) '../madevent >> $k <input_sg.txt'
 
685
         write(26,25) 'cat $k >> log.txt'
 
686
         write(26,25) 'if [[ -e ftn26 ]]; then'
 
687
         write(26,25) '     cp ftn26 ftn25'
 
688
         write(26,25) 'fi'
 
689
         write(26,20) 'else'
 
690
 
 
691
         write(26,25) 'rm -f $k'
 
692
 
 
693
         write(26,'(9x,a,3i8,a)') 'echo "',npoints,max_iter,min_iter,
 
694
     $        '" >& input_sg.txt' 
 
695
c
 
696
c tjs 8/7/2007-JA 8/17/11    Change to request luminocity not accuracy
 
697
c
 
698
         write(26,'(9x,a,e13.5,a)') 'echo "',-goal_lum*1000/mjobs,
 
699
     $        '" >> input_sg.txt'                       !Luminocity
 
700
c         write(26,'(9x,a,e12.3,a)') 'echo "',-goal_lum*mfact(io(np)),
 
701
c     $        '" >> input_sg.txt'
 
702
         write(26,'(9x,a)') 'echo "0" >> input_sg.txt'
 
703
         write(26,'(9x,a)') 'echo "1" >> input_sg.txt'
 
704
 
 
705
         write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
 
706
     &        ' " >> input_sg.txt' !Helicity 0=exact
 
707
 
 
708
         write(26,'(9x,3a)')'echo "',gn(io(np))(2:ip-1),
 
709
     $        '" >>input_sg.txt'
 
710
 
 
711
 
 
712
c         write(26,'(9x,a)') 'echo "1" >> input_sg.txt' !Helicity 0=exact
 
713
 
 
714
c         write(26,'(5x,3a)')'echo "',gn(io(np))(2:ip-1),
 
715
c     $        '" >>input_sg.txt'
 
716
c         write(26,20) 'cp ../../public_sg.sh .'
 
717
c         write(26,20) 'qsub -N $1$j public_sg.sh >> ../../running_jobs'
 
718
         write(26,25) 'if [[ -e ftn26 ]]; then'
 
719
         write(26,25) '     cp ftn26 ftn25'
 
720
         write(26,25) 'fi'
 
721
         write(26,25) '../madevent >> $k <input_sg.txt'
 
722
         write(26,25) 'cat $k >> log.txt'
 
723
         write(26,20) 'fi'
 
724
         write(26,20) 'cd ../'
 
725
c------
 
726
c  tjs  end loop over split process   
 
727
c------
 
728
      enddo  !(ijob, split channel)         
 
729
 
 
730
      enddo !(k  each channel)
 
731
      if (fopened) then
 
732
         call close_bash_file(26)
 
733
      endif
 
734
c      write(26,15) 'end'
 
735
 15   format(a)
 
736
 20   format(5x,a)
 
737
 25   format(10x,a)
 
738
 999  close(26)
 
739
      end
 
740
 
 
741
 
 
742
      subroutine write_gen_grid(goal_lum,ngran,ng,jpoints,gn,xlum,xtot,mfact,xsec,nhel_refine)
 
743
c*****************************************************************************
 
744
c     Writes out scripts for achieving unweighted event goals
 
745
c*****************************************************************************
 
746
      implicit none
 
747
c
 
748
c     Constants
 
749
c
 
750
      include 'run_config.inc'
 
751
      include 'maxconfigs.inc'
 
752
c
 
753
c   global
 
754
c
 
755
      integer max_np,min_iter
 
756
      common/max_np/max_np,min_iter
 
757
c
 
758
c     Arguments
 
759
c
 
760
      double precision goal_lum, xlum(lmaxconfigs), xsec(lmaxconfigs),xtot
 
761
      double precision ngran   !Granularity.... min # points from channel
 
762
      integer jpoints(lmaxconfigs), mfact(lmaxconfigs)
 
763
      integer ng, np, nhel_refine
 
764
      character*(80) gn(lmaxconfigs)
 
765
c
 
766
c     Local
 
767
c
 
768
      integer i,j,k, npoints, ip
 
769
      double precision xt(lmaxconfigs),elimit
 
770
      double precision yerr,ysec,rerr
 
771
      character*72 fname
 
772
      logical fopened
 
773
      double precision rvec
 
774
c-----
 
775
c  Begin Code
 
776
c-----
 
777
 
 
778
c      data ngran /10/
 
779
      fopened=.false.
 
780
c
 
781
c     These random #'s should be changed w/ run
 
782
c
 
783
c      ij=2134
 
784
c      kl = 4321
 
785
      rvec=0d0
 
786
      write(*,*) 'Working on creating ', goal_lum, ' events.'
 
787
      max_np = 1
 
788
      np = max_np   !Flag to open csh file
 
789
      do i=1,ng
 
790
         call ranmar(rvec)
 
791
         ip = index(gn(i),'/')
 
792
         fname = gn(i)(1:ip) // 'gscalefact.dat'
 
793
         open(unit=27,file=fname,status='unknown',err=91)
 
794
         if (goal_lum * xsec(i)/xtot  .ge. rvec*ngran ) then !need events
 
795
            write(*,*) 'Requesting events from ',gn(i)(1:ip-1),xsec(i),xtot/goal_lum
 
796
            if (xsec(i) .gt. xtot*ngran/goal_lum) then
 
797
               write(27,*) 1d0
 
798
            else
 
799
               write(27,*) xtot*ngran/xsec(i)/goal_lum
 
800
            endif
 
801
            npoints = goal_lum * xsec(i) / xtot
 
802
            if (npoints .lt. ngran) npoints = ngran
 
803
            np = np+1
 
804
            if (np .gt. max_np) then
 
805
               if (fopened) then
 
806
                  call close_bash_file(26)
 
807
               endif
 
808
               fopened=.true.
 
809
               call open_bash_file(26)
 
810
               np = 1
 
811
            endif
 
812
            ip = index(gn(i),'/')
 
813
            write(*,*) 'Channel ',gn(i)(2:ip-1), goal_lum * xsec(i) / xtot,
 
814
     $           npoints
 
815
 
 
816
            ip = index(gn(i),'/')
 
817
            write(26,'(2a)') 'j=',gn(i)(1:ip-1)
 
818
c
 
819
c           Now write the commands
 
820
c      
 
821
            write(26,20) 'if [[ ! -e $j ]]; then'
 
822
            write(26,25) 'mkdir $j'
 
823
            write(26,20) 'fi'
 
824
            write(26,20) 'cd $j'
 
825
            write(26,20) 'rm -f $k'
 
826
c
 
827
c     Now I'll add a check to make sure the grid has been
 
828
c     adjusted  (ftn99 or ftn25 exist)
 
829
c
 
830
            write(26,20) 'if  [[ -e ftn26 ]]; then'
 
831
            write(26,25) 'cp ftn26 ftn25'
 
832
            write(26,20) 'fi'
 
833
 
 
834
            write(26,20) 'if [[ ! -e ftn25 ]]; then'
 
835
 
 
836
 
 
837
            write(26,'(9x,a,3i8,a)') 'echo "',max(npoints,min_events),
 
838
     $           max_iter,min_iter,'" >& input_sg.txt' 
 
839
c
 
840
c     tjs 8/7/2007  Allow stop when have enough events
 
841
c
 
842
            write(*,*) "Cross section",i,xsec(i),mfact(i)
 
843
            write(26,'(9x,a,e13.5,a)') 'echo "',-npoints/xsec(i),
 
844
     $        '" >> input_sg.txt'                       !Luminocity
 
845
            write(26,'(9x,a)') 'echo "2" >> input_sg.txt' !Grid Adjustment
 
846
            write(26,'(9x,a)') 'echo "1" >> input_sg.txt' !Suppression
 
847
            write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
 
848
     &           ' " >> input_sg.txt' !Helicity 0=exact
 
849
            write(26,'(9x,3a)')'echo "',gn(i)(2:ip-1),
 
850
     $           '" >>input_sg.txt'
 
851
            write(26,25) '../madevent >> $k <input_sg.txt'
 
852
            write(26,25) 'cat $k >> log.txt'
 
853
            write(26,25) 'if [[ -e ftn26 ]]; then'
 
854
            write(26,25) '     cp ftn26 ftn25'
 
855
            write(26,25) 'fi'
 
856
            write(26,20) 'else'
 
857
 
 
858
            write(26,25) 'rm -f $k'
 
859
            
 
860
            write(26,'(9x,a,3i8,a)') 'echo "',max(npoints,min_events),
 
861
     $           max_iter,min_iter,'" >& input_sg.txt' 
 
862
c
 
863
c tjs 8/7/2007    Change to request events not accuracy
 
864
c
 
865
            write(26,'(9x,a,e13.5,a)') 'echo "',-npoints / xsec(i),
 
866
     $           '" >> input_sg.txt' ! Luminocity
 
867
            write(26,'(9x,a)') 'echo "0" >> input_sg.txt'
 
868
            write(26,'(9x,a)') 'echo "1" >> input_sg.txt'
 
869
 
 
870
            write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
 
871
     &           ' " >> input_sg.txt' !Helicity 0=exact
 
872
 
 
873
            write(26,'(9x,3a)')'echo "',gn(i)(2:ip-1),
 
874
     $           '" >>input_sg.txt'
 
875
 
 
876
            write(26,25) 'if [[ -e ftn26 ]]; then'
 
877
            write(26,25) '     cp ftn26 ftn25'
 
878
            write(26,25) 'fi'
 
879
            write(26,25) '../madevent >> $k <input_sg.txt'
 
880
            write(26,25) 'cat $k >> log.txt'
 
881
            write(26,20) 'fi'
 
882
            write(26,20) 'cd ../'
 
883
         else    !No events from this channel
 
884
            write(*,*) 'Skipping channel:',gn(i)(1:ip-1),xsec(i)*goal_lum/xtot,rvec
 
885
            write(27,*) 0d0
 
886
         endif
 
887
         close(27)
 
888
 91      cycle
 
889
      enddo
 
890
      call close_bash_file(26)
 
891
 15   format(a)
 
892
 20   format(5x,a)
 
893
 25   format(10x,a)
 
894
 999  close(26)
 
895
      close(27)
 
896
      end
 
897
 
 
898
 
 
899
      subroutine sort2(array,aux1,n)
 
900
      implicit none
 
901
! Arguments
 
902
      integer n
 
903
      integer aux1(n)
 
904
      double precision array(n)
 
905
!  Local Variables
 
906
      integer i,k
 
907
      double precision temp
 
908
      logical done
 
909
 
 
910
!-----------
 
911
! Begin Code
 
912
!-----------
 
913
      do i=n-1,1,-1
 
914
         done = .true.
 
915
         do k=1,i
 
916
            if (array(k) .lt. array(k+1)) then
 
917
               temp = array(k)
 
918
               array(k) = array(k+1)
 
919
               array(k+1) = temp
 
920
               temp = aux1(k)
 
921
               aux1(k) = aux1(k+1)
 
922
               aux1(k+1) = temp
 
923
               done = .false.
 
924
            end if
 
925
         end do
 
926
         if (done) return
 
927
      end do
 
928
      end 
 
929
 
 
930
 
 
931
      subroutine get_xsec_log(xsec,xerr,eff,xmax)
 
932
c*********************************************************************
 
933
c     Reads from output file, gets cross section and maxwgt from
 
934
c     first two iterations
 
935
c*********************************************************************
 
936
      implicit none
 
937
c
 
938
c     Arguments
 
939
c
 
940
      double precision xsec(2),xerr(2),eff(2),xmax(2)
 
941
c
 
942
c     Local
 
943
c     
 
944
      character*78 buff
 
945
      integer i
 
946
c-----
 
947
c  Begin Code
 
948
c-----
 
949
      xsec(1) = 0d0
 
950
      xerr(1) = 0d0
 
951
      xmax(1) = 0d0
 
952
      do while (.true.)
 
953
         read(25,'(a80)',err=99) buff
 
954
         if (buff(1:4) .eq. 'Iter') then
 
955
            read(buff(11:16),'(i5)') i
 
956
            if (i .eq. 1 .or. i .eq. 2) then
 
957
               read(buff(61:70),*) xmax(i)
 
958
               read(buff(21:33),*) xsec(i)
 
959
               xmax(i)=xmax(i)/xsec(i)
 
960
c               read(buff(48:59),*) xerr(i)
 
961
c               read(buff(48:59),*) xmax(i)
 
962
            endif
 
963
            read(25,'(a80)',err=99) buff
 
964
            read(buff(1:6),'(i5)') i
 
965
            if (i .eq. 1 .or. i .eq. 2) then
 
966
               read(buff(6:17),*) xsec(i)
 
967
               read(buff(20:31),*) xerr(i)
 
968
               read(buff(34:40),*) eff(i)               
 
969
            endif
 
970
            write(*,'(i4,4f12.3)') i,xsec(i),xerr(i),eff(i),xmax(i)
 
971
         endif
 
972
      enddo
 
973
 99   end
 
974
 
 
975
 
 
976
 
 
977