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*****************************************************************************
11
parameter (rfile='results.dat')
13
parameter (symfile='symfact.dat')
15
parameter (max_amps=9999)
17
parameter (maxpara=1000)
19
include 'run_config.inc'
29
double precision xsec(max_amps), xerr(max_amps)
30
double precision xerru(max_amps),xerrc(max_amps)
31
double precision xmax(max_amps), eff(max_amps)
32
double precision xlum(max_amps)
33
double precision ysec, yerr, yeff, ymax
34
double precision tsec, terr, teff, tmax, xi
35
integer nw(max_amps), nevents(max_amps), maxit
36
integer nunwgt(max_amps)
37
character*80 fname, gname(max_amps)
39
double precision xtot,errtot,err_goal
40
double precision errtotu,errtotc
41
integer mfact(max_amps)
42
logical parallel, gen_events
43
character*20 param(maxpara),value(maxpara)
44
integer npara, nreq, ngran
46
logical Gridpack,gridrun
51
call load_para(npara,param,value)
52
call get_logical(npara,param,value," gridpack ",gridpack,.false.)
53
if (.not. Gridpack) then
54
write(*,'(a,a)')'Enter fractional accuracy (<1)',
55
& ', or number events (>1) and max processes per job'
56
read(*,*) err_goal, max_np
58
if (err_goal .lt. 1) then
59
write(*,'(a,f8.2,a)') 'Running for accuracy of ',
62
elseif (err_goal .gt. 1) then
63
write(*,'(a,f9.0,a)') 'Generating ',err_goal,
64
& ' unweighted events.'
66
err_goal = err_goal * 1.2 !Extra factor to ensure works
68
write(*,*) 'Error, need non_zero goal'
73
call get_integer(npara,param,value," gevents " ,nreq ,2000 )
75
call get_integer(npara,param,value," gseed " ,iseed ,4321 )
76
call get_integer(npara,param,value," ngran " ,ngran , -1)
77
if (ngran.eq.-1) ngran = int(sqrt(real(nreq)))
78
write(*,*) "Running on Grid to generate ",nreq," additional events"
79
write(*,*) " with granularity equal to ",ngran
82
c Modified to allow for more sequences
83
c iseed can be between 0 and 31328*30081
84
c before patern repeats
86
ij=1802 + mod(iseed,30081)
87
kl=9373 + (iseed/31328)
88
c write(*,'($a,i6,a3,i6)') 'Using random seed offsets',jconfig," : ",ioffset
89
c write(*,*) ' with seed', iseed
90
do while (ij .gt. 31328)
93
do while (kl .gt. 30081)
96
write(*,*) "Using random seeds",ij,kl
99
open(unit=15,file=symfile,status='old',err=999)
106
read(15,*,err=99,end=99) xi,j
109
if ( (xi-int(xi+.01)) .lt. 1d-5) then
112
write(fname,'(a,i1,a,a)') 'G',n,'/',rfile
113
else if (n .lt. 100) then
114
write(fname,'(a,i2,a,a)') 'G',n,'/',rfile
115
else if (n .lt. 1000) then
116
write(fname,'(a,i3,a,a)') 'G',n,'/',rfile
117
else if (n .lt. 10000) then
118
write(fname,'(a,i4,a,a)') 'G',n,'/',rfile
122
write(fname,'(a,f5.3,a,a)') 'G',xi,'/',rfile
123
else if (xi .lt. 100) then
124
write(fname,'(a,f6.3,a,a)') 'G',xi,'/',rfile
125
else if (xi .lt. 1000) then
126
write(fname,'(a,f7.3,a,a)') 'G',xi,'/',rfile
127
else if (xi .lt. 10000) then
128
write(fname,'(a,f8.3,a,a)') 'G',xi,'/',rfile
130
c write(*,*) 'log name ',fname
142
c Read in integration data from run
144
open(unit=25,file=fname,status='old',err=95)
145
read(25,*) xsec(i),xerru(i),xerrc(i),nevents(i),nw(i),maxit,
147
if (xsec(i) .eq. 0d0) xlum(i)=1d99 !zero cross section
148
xlum(i) = xlum(i)/1000 !convert to fb^-1
149
xerr(i)=sqrt(xerru(i)**2+xerrc(i)**2)
157
read(25,*,err=92) l,ysec,yerr,yeff,ymax
158
if (k .gt. 1) tmax = max(tmax,ymax)
163
92 maxit = k-1 !In case of error reading file
165
xerr(i)=sqrt(terr)/maxit
170
xsec(i) = xsec(i)*mfact(i)
171
xerr(i) = xerr(i)*mfact(i)
172
xerru(i) = xerru(i)*mfact(i)
173
xerrc(i) = xerrc(i)*mfact(i)
174
xlum(i) = xlum(i)/mfact(i)
176
eff(i)= xerr(i)*sqrt(real(nevents(i)))/xsec(i)
177
errtotu = errtotu+(xerru(i))**2
178
errtotc = errtotc+(xerrc(i))
179
c xtot = xtot+ xsec(i)*mfact(i)
180
c eff(i)= xerr(i)*sqrt(real(nevents(i)))/xsec(i)
181
c errtot = errtot+(mfact(i)*xerr(i))**2
183
c write(*,*) i,maxit,xsec(i), eff(i)
185
c i=i-1 !This is for case w/ B.W. and optimization
189
errtot=sqrt(errtotc**2+errtotu)
190
if ( .not. gen_events) then
191
call write_bash(xsec,xerru,xerrc,xtot,mfact,err_goal,
194
open(unit=25,file='../results.dat',status='old',err=199)
195
write(*,'(a,e12.3)') 'Reading total xsection ',xtot
197
write(*,'(e12.3)') xtot
200
call write_gen_grid(err_goal,dble(ngran),i,nevents,gname,xlum,xtot,mfact,xsec)
202
call write_gen(err_goal,i,nevents,gname,xlum,xtot,mfact,xsec)
206
999 write(*,*) 'error'
210
subroutine write_bash(xsec,xerru,xerrc,xtot,
211
$ mfact,err_goal,ng,jpoints,gn)
212
c*****************************************************************************
213
c Writes out bash commands for running each channel as needed.
214
c*****************************************************************************
219
include 'run_config.inc'
221
parameter (max_amps=9999)
223
c parameter (max_np = 30)
232
double precision xsec(max_amps), xerru(max_amps),xerrc(max_amps)
233
double precision err_goal,xtot
234
integer mfact(max_amps),jpoints(max_amps)
236
character*(80) gn(max_amps)
240
integer i,j,k, io(max_amps), npoints, ip, np
241
double precision xt(max_amps),elimit
242
double precision yerr,ysec,rerr
250
if (mfact(j) .gt. 0) k=k+1
252
xt(j)= sqrt((xerru(j)+xerrc(j)**2)*mfact(j)) !sort by error
255
c Let's redetermine err_goal based on luminosity
257
write(*,*) 'Cross section pb',xtot
258
write(*,*) 'Desired Goal',err_goal
259
write(*,*) 'Total Error',err_goal
260
c elimit = err_goal*xtot/sqrt(real(k)) !Equal contributions from each
261
elimit = err_goal*xtot/real(k) !Equal contributions from each
266
do while( xt(k) .gt. abs(elimit)) !abs is just in case elimit<0 by mistake
272
c rerr = rerr+xt(j)**2
276
c write(*,*) 'Number of diagrams to fix',k
278
c Now readjust because most don't contribute
280
elimit = sqrt((err_goal*xtot)**2 - rerr)/sqrt(real(k))
286
c yerr = xerr(io(i))*mfact(io(i))
288
c write(*,*) i,xt(i),elimit
289
if (yerr .gt. elimit) then
291
ysec = xsec(io(i)) + yerr
292
npoints=(0.2d0)*jpoints(io(i))*(yerr/elimit)**2
293
npoints = max(npoints,min_events)
294
npoints = min(npoints,max_events)
295
c np = np + 3*npoints
297
if (np .gt. max_np) then
299
write(26,15) 'rm -f run.$script >&/dev/null'
300
write(26,15) 'touch done.$script >&/dev/null'
304
call open_bash_file(26)
309
ip = index(gn(io(i)),'/')
310
write(*,*) 'Channel ',gn(io(i))(2:ip-1),
311
$ yerr, jpoints(io(i)),npoints
313
ip = index(gn(io(i)),'/')
314
write(26,'(2a)') 'j=',gn(io(i))(1:ip-1)
316
if (io(i) .lt. 10) then
317
write(26,'(a,i1)') 'j=G',io(i)
318
elseif (io(i) .lt. 100) then
319
write(26,'(a,i2)') 'j=G',io(i)
320
elseif (io(i) .lt. 1000) then
321
write(26,'(a,i3)') 'j=G',io(i)
322
elseif (io(i) .lt. 10000) then
323
write(26,'(a,i4)') 'j=G',io(i)
327
c Determine estimates for getting the desired accuracy
331
c Now write the commands
333
write(26,20) 'echo $j'
334
write(26,20) 'if [[ ! -e $j ]]; then'
335
write(26,25) 'mkdir $j'
338
write(26,20) 'rm -f $k'
340
write(26,'(5x,a,2i8,a)') 'echo "',npoints,max_iter,
341
$ '" >& input_sg.txt'
342
write(26,'(5x,a,f8.3,a)') 'echo "',max(elimit/ysec,0.001d0),
343
$ '" >> input_sg.txt'
344
write(26,'(5x,a)') 'echo "2" >> input_sg.txt' !Grid
345
write(26,'(5x,a)') 'echo "1" >> input_sg.txt' !Suppress
346
write(26,'(5x,a,i4,a)') 'echo "',nhel_refine,
347
& '" >> input_sg.txt' !Helicity
348
write(26,'(5x,3a)')'echo "',gn(io(i))(2:ip-1),
350
write(26,20) 'time ../madevent >> $k <input_sg.txt'
351
write(26,20) 'mv ftn26 ftn25'
352
c write(26,20) 'rm ftn26'
353
write(26,20) 'cat $k >> log.txt'
354
write(26,20) 'cd ../'
356
enddo !Loop over diagrams
358
write(26,15) 'rm -f run.$script >&/dev/null'
359
write(26,15) 'touch done.$script >&/dev/null'
370
subroutine open_bash_file(lun)
371
c***********************************************************************
372
c Opens bash file for looping including standard header info
373
c which can be used with pbs, or stand alone
374
c***********************************************************************
379
include 'run_config.inc'
397
write(fname(5:5),'(i1)') ic
398
elseif (ic .lt. 100) then
399
write(fname(5:6),'(i2)') ic
400
elseif (ic .lt. 1000) then
401
write(fname(5:7),'(i3)') ic
402
elseif (ic .lt. 10000) then
403
write(fname(5:8),'(i4)') ic
405
write(*,*) 'Error too many jobs in gen_ximprove',ic
408
write(*,*) 'Opening file ',fname
410
open (unit=26, file = fname, status='unknown')
411
write(26,15) '#!/bin/bash'
412
write(26,15) '#PBS -q ' // PBS_QUE
413
write(26,15) '#PBS -o /dev/null'
414
write(26,15) '#PBS -e /dev/null'
415
write(26,15) 'if [[ "$PBS_O_WORKDIR" != "" ]]; then'
416
write(26,15) ' cd $PBS_O_WORKDIR'
418
write(26,15) 'k=run1_app.log'
419
write(lun,15) 'script=' // fname
420
write(lun,15) 'rm -f wait.$script >& /dev/null'
421
write(lun,15) 'touch run.$script'
425
subroutine close_bash_file(lun)
426
c***********************************************************************
428
c***********************************************************************
448
c write(lun,'(a)') ')'
450
c Now write the commands
452
c write(lun,20) 'echo $i'
453
c write(lun,20) 'j=G$i'
454
c write(lun,20) 'if (! -e $j) then'
455
c write(lun,25) 'mkdir $j'
456
c write(lun,20) 'endif'
457
c write(lun,20) 'cd $j'
458
c write(lun,20) 'rm -f ftn25 ftn99'
459
c write(lun,20) 'rm -f $k'
460
c write(lun,20) 'cat ../input_app.txt >& input_app.txt'
461
c write(lun,20) 'echo $i >> input_app.txt'
463
c write(lun,20) 'cp ../../public.sh .'
464
c write(lun,20) 'qsub -N $1$i public.sh >> ../../running_jobs'
466
c write(lun,20) 'time ../madevent > $k <input_app.txt'
467
c write(lun,20) 'rm -f ftn25 ftn99'
468
c write(lun,20) 'cp $k log.txt'
470
c write(lun,20) 'cd ../'
471
c write(lun,15) 'end'
480
subroutine write_gen(goal_lum,ng,jpoints,gn,xlum,xtot,mfact,xsec)
481
c*****************************************************************************
482
c Writes out scripts for achieving unweighted event goals
483
c*****************************************************************************
488
include 'run_config.inc'
490
parameter (max_amps=9999)
496
c integer max_np !now set in run_config.inc
497
c parameter (max_np = 5) !number of channels/job
502
double precision goal_lum, xlum(max_amps), xsec(max_amps),xtot
503
integer jpoints(max_amps), mfact(max_amps)
505
character*(80) gn(max_amps)
509
integer i,j,k, io(max_amps), npoints, ip, nfiles,ifile,npfile
510
double precision xt(max_amps),elimit
511
double precision yerr,ysec,rerr
517
write(*,*) 'Working on creating ', goal_lum, ' events.'
518
goal_lum = goal_lum/(xtot*1000) !Goal luminosity in fb^-1
519
write(*,*) 'Effective Luminosity', goal_lum, ' fb^-1.'
523
xt(j)= goal_lum/xlum(j) !sort by events_needed/have.
524
write(*,*) j,xlum(j),xt(j)
526
c write(*,*) 'Number of channels',ng,k
528
c Let's redetermine err_goal based on luminosity
534
do while( xt(k) .gt. abs(elimit)) !elimit should be >0
535
write(*,*) 'Improving ',k,gn(io(k)),xt(k)
539
c write(*,*) 'Number of diagrams to fix',k
541
c Now readjust because most don't contribute
547
c Want to write channels so that heaviest one (with largest error)
548
c gets grouped with least heavy channels. Complicated ordering for this
549
c follows. np is the present process number.
552
if(mod(k,max_np).gt.0) nfiles=nfiles+1
560
npoints=0.2*jpoints(io(np))*(yerr/elimit)
561
npoints = max(npoints,min_events)
562
npoints = min(npoints,max_events)
564
c np = np + 5*npoints
566
np = nfiles*npfile+1-ifile
567
if (np .gt. k .or. ifile.eq.0) then
569
write(26,15) 'rm -f run.$script >& /dev/null'
570
write(26,15) 'touch done.$script >& /dev/null'
574
call open_bash_file(26)
580
ip = index(gn(io(np)),'/')
581
write(*,*) 'Channel ',gn(io(np))(2:ip-1),
582
$ yerr, jpoints(io(np)),npoints
584
ip = index(gn(io(np)),'/')
585
write(26,'(2a)') 'j=',gn(io(np))(1:ip-1)
587
c Now write the commands
589
write(26,20) 'echo $j'
590
write(26,20) 'if [[ ! -e $j ]]; then'
591
write(26,25) 'mkdir $j'
594
write(26,20) 'rm -f $k'
596
c Now I'll add a check to make sure the grid has been
597
c adjusted (ftn99 or ftn25 exist)
599
write(26,20) 'if [[ -e ftn26 ]]; then'
600
write(26,25) 'cp ftn26 ftn25'
603
write(26,20) 'if [[ ! -e ftn25 ]]; then'
606
write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
607
$ '" >& input_sg.txt'
609
c tjs 8/7/2007 Allow stop when have enough events
611
write(*,*) "Cross section",i,io(np),xsec(io(np)),mfact(io(np))
612
write(26,'(9x,a,f11.3,a)') 'echo "',-goal_lum*xsec(io(np))*1000,
613
$ '" >> input_sg.txt' !Accuracy
614
write(26,'(9x,a)') 'echo "2" >> input_sg.txt' !Grid Adjustment
615
write(26,'(9x,a)') 'echo "1" >> input_sg.txt' !Suppression
616
write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
617
& ' " >> input_sg.txt' !Helicity 0=exact
618
write(26,'(9x,3a)')'echo "',gn(io(np))(2:ip-1),
620
write(26,25) 'time ../madevent >> $k <input_sg.txt'
621
write(26,25) 'cat $k >> log.txt'
622
write(26,25) 'if [[ -e ftn26 ]]; then'
623
write(26,25) ' cp ftn26 ftn25'
627
write(26,25) 'rm -f $k'
629
write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
630
$ '" >& input_sg.txt'
632
c tjs 8/7/2007 Change to request events not accuracy
634
write(26,'(9x,a,f11.3,a)') 'echo "',-goal_lum*xsec(io(np))*1000,
635
$ '" >> input_sg.txt' !Accuracy
636
c write(26,'(9x,a,e12.3,a)') 'echo "',-goal_lum*mfact(io(np)),
637
c $ '" >> input_sg.txt'
638
write(26,'(9x,a)') 'echo "0" >> input_sg.txt'
639
write(26,'(9x,a)') 'echo "1" >> input_sg.txt'
641
write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
642
& ' " >> input_sg.txt' !Helicity 0=exact
644
write(26,'(9x,3a)')'echo "',gn(io(np))(2:ip-1),
648
c write(26,'(9x,a)') 'echo "1" >> input_sg.txt' !Helicity 0=exact
650
c write(26,'(5x,3a)')'echo "',gn(io(np))(2:ip-1),
651
c $ '" >>input_sg.txt'
652
c write(26,20) 'cp ../../public_sg.sh .'
653
c write(26,20) 'qsub -N $1$j public_sg.sh >> ../../running_jobs'
654
write(26,25) 'if [[ -e ftn26 ]]; then'
655
write(26,25) ' cp ftn26 ftn25'
657
write(26,25) 'time ../madevent >> $k <input_sg.txt'
658
write(26,25) 'cat $k >> log.txt'
660
write(26,20) 'cd ../'
663
write(26,15) 'rm -f run.$script'
664
write(26,15) 'touch done.$script'
675
subroutine write_gen_grid(goal_lum,ngran,ng,jpoints,gn,xlum,xtot,mfact,xsec)
676
c*****************************************************************************
677
c Writes out scripts for achieving unweighted event goals
678
c*****************************************************************************
683
include 'run_config.inc'
685
parameter (max_amps=9999)
694
double precision goal_lum, xlum(max_amps), xsec(max_amps),xtot
695
double precision ngran !Granularity.... min # points from channel
696
integer jpoints(max_amps), mfact(max_amps)
698
character*(80) gn(max_amps)
702
integer i,j,k, npoints, ip
703
double precision xt(max_amps),elimit
704
double precision yerr,ysec,rerr
707
double precision rvec
715
c These random #'s should be changed w/ run
720
write(*,*) 'Working on creating ', goal_lum, ' events.'
722
np = max_np !Flag to open csh file
725
ip = index(gn(i),'/')
726
fname = gn(i)(1:ip) // 'gscalefact.dat'
727
open(unit=27,file=fname,status='unknown',err=91)
728
if (goal_lum * xsec(i)/xtot .ge. rvec*ngran ) then !need events
729
write(*,*) 'Requesting events from ',gn(i)(1:ip-1),xsec(i),xtot/goal_lum
730
if (xsec(i) .gt. xtot*ngran/goal_lum) then
733
write(27,*) xtot*ngran/xsec(i)/goal_lum
735
npoints = goal_lum * xsec(i) / xtot
736
if (npoints .lt. min_gevents_wu) npoints = min_gevents_wu
738
if (np .gt. max_np) then
740
write(26,15) 'rm -f run.$script >& /dev/null'
741
write(26,15) 'touch done.$script >& /dev/null'
745
call open_bash_file(26)
748
ip = index(gn(i),'/')
749
write(*,*) 'Channel ',gn(i)(2:ip-1),
750
$ yerr, jpoints(i),npoints
752
ip = index(gn(i),'/')
753
write(26,'(2a)') 'j=',gn(i)(1:ip-1)
755
c Now write the commands
757
write(26,20) 'echo $j'
758
write(26,20) 'if [[ ! -e $j ]]; then'
759
write(26,25) 'mkdir $j'
762
write(26,20) 'rm -f $k'
764
c Now I'll add a check to make sure the grid has been
765
c adjusted (ftn99 or ftn25 exist)
767
write(26,20) 'if [[ -e ftn26 ]]; then'
768
write(26,25) 'cp ftn26 ftn25'
771
write(26,20) 'if [[ ! -e ftn25 ]]; then'
774
write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
775
$ '" >& input_sg.txt'
777
c tjs 8/7/2007 Allow stop when have enough events
779
write(*,*) "Cross section",i,xsec(i),mfact(i)
780
write(26,'(9x,a,f11.3,a)') 'echo "',-npoints*1.2,
781
$ '" >> input_sg.txt' !Accuracy
782
write(26,'(9x,a)') 'echo "2" >> input_sg.txt' !Grid Adjustment
783
write(26,'(9x,a)') 'echo "1" >> input_sg.txt' !Suppression
784
write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
785
& ' " >> input_sg.txt' !Helicity 0=exact
786
write(26,'(9x,3a)')'echo "',gn(i)(2:ip-1),
788
write(26,25) 'time ../madevent >> $k <input_sg.txt'
789
write(26,25) 'cat $k >> log.txt'
790
write(26,25) 'if [[ -e ftn26 ]]; then'
791
write(26,25) ' cp ftn26 ftn25'
795
write(26,25) 'rm -f $k'
797
write(26,'(9x,a,2i8,a)') 'echo "',npoints,max_iter,
798
$ '" >& input_sg.txt'
800
c tjs 8/7/2007 Change to request events not accuracy
802
write(26,'(9x,a,f11.3,a)') 'echo "',-npoints*1.2,
803
$ '" >> input_sg.txt' !Accuracy
804
write(26,'(9x,a)') 'echo "0" >> input_sg.txt'
805
write(26,'(9x,a)') 'echo "1" >> input_sg.txt'
807
write(26,'(9x,a,i4,a)') 'echo "',nhel_refine,
808
& ' " >> input_sg.txt' !Helicity 0=exact
810
write(26,'(9x,3a)')'echo "',gn(i)(2:ip-1),
813
write(26,25) 'if [[ -e ftn26 ]]; then'
814
write(26,25) ' cp ftn26 ftn25'
816
write(26,25) 'time ../madevent >> $k <input_sg.txt'
817
write(26,25) 'cat $k >> log.txt'
819
write(26,20) 'cd ../'
820
else !No events from this channel
821
write(*,*) 'Skipping channel:',gn(i)(1:ip-1),xsec(i)*goal_lum/xtot,rvec
827
write(26,15) 'rm -f run.$script'
828
write(26,15) 'touch done.$script'
838
subroutine sort2(array,aux1,n)
843
double precision array(n)
846
double precision temp
855
if (array(k) .lt. array(k+1)) then
857
array(k) = array(k+1)
870
subroutine get_xsec_log(xsec,xerr,eff,xmax)
871
c*********************************************************************
872
c Reads from output file, gets cross section and maxwgt from
873
c first two iterations
874
c*********************************************************************
879
double precision xsec(2),xerr(2),eff(2),xmax(2)
892
read(25,'(a80)',err=99) buff
893
if (buff(1:4) .eq. 'Iter') then
894
read(buff(11:16),'(i5)') i
895
if (i .eq. 1 .or. i .eq. 2) then
896
read(buff(61:70),*) xmax(i)
897
read(buff(21:33),*) xsec(i)
898
xmax(i)=xmax(i)/xsec(i)
899
c read(buff(48:59),*) xerr(i)
900
c read(buff(48:59),*) xmax(i)
902
read(25,'(a80)',err=99) buff
903
read(buff(1:6),'(i5)') i
904
if (i .eq. 1 .or. i .eq. 2) then
905
read(buff(6:17),*) xsec(i)
906
read(buff(20:31),*) xerr(i)
907
read(buff(34:40),*) eff(i)
909
write(*,'(i4,4f12.3)') i,xsec(i),xerr(i),eff(i),xmax(i)
added
removed
Template/Source/gen_ximprove.f
