« back to all changes in this revision
Viewing changes to Template/Source/dsample.f
- browse files at revision 50.9.3
- compare with another revision
- download diff
- download tarball
- view history from revision 50.9.3
- Committer: Johan Alwall
- Date: 2010-10-30 03:26:37 UTC
- mto: (50.16.2 new_color_ordering) (143.1.2 FKS5) (74.1.37 aloha_paper)
- mto: This revision was merged to the branch mainline in revision 80.
- Revision ID: alwall@alwall-laptop-20101030032637-fchuklnv5s6fww65
Added Template and HELAS into bzr
- files added:
- HELAS
- HELAS/lib
- Template
- Template/Cards
- Template/Events
- Template/HTML
- Template/Source
- Template/Source/CERNLIB
- Template/Source/DHELAS
- Template/Source/MODEL
- Template/Source/PDF
- Template/SubProcesses
- Template/bin
- Template/lib
- Template/lib/Pdfdata
added
removed
Template/Source/dsample.f
1
subroutine sample_full(ndim,ncall,itmax,dsig,ninvar,nconfigs)
2
c**************************************************************************
3
c Driver for sample which does complete integration
4
c This is done in double precision, and should be told the
5
c number of possible phasespace choices.
6
c Arguments:
7
c ndim Number of dimensions for integral(number or random #'s/point)
8
c ncall Number of times to evaluate the function/iteration
9
c itmax Number of iterations
10
c ninvar Number of invarients to keep grids on (s,t,u, s',t' etc)
11
c nconfigs Number of different pole configurations
12
c dsig Function to be integrated
13
c**************************************************************************
14
implicit none
15
include 'genps.inc'
16
c
17
c Arguments
18
c
19
integer ndim,ncall,itmax,ninvar,nconfigs
20
external dsig
21
double precision dsig
22
c
23
c Local
24
c
25
double precision x(maxinvar),wgt,p(4*maxdim/3+14)
26
double precision tdem, chi2
27
integer ievent,kevent,nwrite,iter,nun
28
integer jmax,i,j,ipole
29
integer itmax_adjust
30
c
31
c External
32
c
33
integer n_unwgted
34
external n_unwgted
35
c
36
c Global
37
c
38
integer nsteps
39
character*40 result_file,where_file
40
common /sample_status/result_file,where_file,nsteps
41
double precision fx
42
common /to_fx/ fx
43
44
integer mincfig, maxcfig
45
common/to_configs/mincfig, maxcfig
46
47
double precision xmean(99),xsigma(99),xwmax(99),xeff(99)
48
common/to_iterations/xmean, xsigma, xwmax, xeff
49
50
double precision accur
51
common /to_accuracy/accur
52
53
double precision twgt, maxwgt,swgt(maxevents)
54
integer lun, nw
55
common/to_unwgt/twgt, maxwgt, swgt, lun, nw
56
57
integer nzoom
58
double precision tx(1:3,maxinvar)
59
common/to_xpoints/tx, nzoom
60
61
double precision xzoomfact
62
common/to_zoom/ xzoomfact
63
64
double precision tmean, tsigma
65
integer dim, events, itm, kn, cur_it, invar, configs
66
common /sample_common/
67
. tmean, tsigma, dim, events, itm, kn, cur_it, invar, configs
68
69
integer use_cut
70
common /to_weight/use_cut
71
72
integer icor
73
common/to_correlated/icor
74
75
integer neventswritten
76
common /to_eventswritten/ neventswritten
77
78
c
79
c External
80
c
81
logical pass_point
82
c
83
c Data
84
c
85
data result_file,where_file,nsteps/'SAMPLE','WHERE.AMI',100/
86
data accur/-1d0/
87
data mincfig /1/
88
data maxcfig /1/
89
data twgt/-1d0/ !Dont write out events
90
data lun/27/ !Unit number for events
91
data maxwgt/0d0/
92
data nw/0/ !Number of events written
93
94
95
c-----
96
c Begin Code
97
c-----
98
ievent = 0
99
kevent = 0
100
nzoom = 0
101
xzoomfact = 1d0
102
if (nsteps .lt. 1) nsteps=1
103
nwrite = itmax*ncall/nsteps
104
c open(unit=66,file='.sample_warn',status='unknown')
105
c write(66,*) 'Warnings from sample run.',itmax,ncall
106
c close(66)
107
call sample_init(ndim,ncall,itmax,ninvar,nconfigs)
108
call graph_init
109
do i=1,itmax
110
xmean(i)=0d0
111
xsigma(i)=0d0
112
enddo
113
c mincfig=1
114
c maxcfig=nconfigs
115
wgt = 0d0
116
c
117
c Main Integration Loop
118
c
119
iter = 1
120
do while(iter .le. itmax)
121
c
122
c Get integration point
123
c
124
call sample_get_config(wgt,iter,ipole)
125
if (iter .le. itmax) then
126
ievent=ievent+1
127
call x_to_f_arg(ndim,ipole,mincfig,maxcfig,ninvar,wgt,x,p)
128
if (pass_point(p)) then
129
xzoomfact = 1d0
130
fx = dsig(p,wgt) !Evaluate function
131
if (xzoomfact .gt. 0d0) then
132
wgt = wgt*fx*xzoomfact
133
else
134
wgt = -xzoomfact
135
endif
136
if (wgt .gt. 0d0) call graph_point(p,wgt) !Update graphs
137
else
138
fx =0d0
139
wgt=0d0
140
endif
141
if (nzoom .le. 0) then
142
call sample_put_point(wgt,x(1),iter,ipole) !Store result
143
else
144
nzoom = nzoom -1
145
ievent=ievent-1
146
endif
147
endif
148
if (wgt .gt. 0d0) kevent=kevent+1
149
c
150
c Write out progress/histograms
151
c
152
if (kevent .ge. nwrite) then
153
nwrite = nwrite+ncall*itmax/nsteps
154
nwrite = min(nwrite,ncall*itmax)
155
call graph_store
156
endif
157
99 enddo
158
c
159
c All done
160
c
161
open(unit=66,file='results.dat',status='unknown')
162
i=1
163
do while(xmean(i) .ne. 0 .and. i .lt. cur_it)
164
i=i+1
165
enddo
166
cur_it = i
167
i = cur_it - 3
168
if (i .gt. 0) then
169
tmean = 0d0
170
tsigma = 0d0
171
tdem = 0d0
172
do while (xmean(i) .ne. 0 .and. i .lt. cur_it)
173
tmean = tmean+xmean(i)*xmean(i)**2/xsigma(i)**2
174
tdem = tdem+xmean(i)**2/xsigma(i)**2
175
tsigma = tsigma + xmean(i)**2/ xsigma(i)**2
176
i=i+1
177
enddo
178
tmean = tmean/tsigma
179
tsigma= tmean/sqrt(tsigma)
180
c nun = n_unwgted()
181
182
nun = neventswritten
183
184
chi2 = 0d0
185
do i = cur_it-3,cur_it-1
186
chi2 = chi2+(xmean(i)-tmean)**2/xsigma(i)**2
187
enddo
188
chi2 = chi2/2d0 !Since using only last 3, n-1=2
189
write(*,'(a)') '-----------------------------------------------------'
190
write(*,'(a)') '---------------------------'
191
write(*,'(a,e12.4)') ' Results Last 3 iters: Integral = ',tmean
192
write(*,'(25x,a,e12.4)') 'Std dev = ',tsigma
193
write(*,'(17x,a,f12.4)') 'Chi**2 per DoF. =',chi2
194
write(*,'(a)') '-----------------------------------------------------'
195
write(*,'(a)') '---------------------------'
196
197
if (nun .lt. 0) nun=-nun !Case when wrote maximun number allowed
198
if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
199
if (icor .eq. 0) then
200
write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,tsigma,0.0,kevent,nw,
201
& cur_it-1,nun, nun/max(tmean,1d-99)
202
else
203
write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,0.0,tsigma,kevent,nw,
204
& cur_it-1,nun, nun/max(tmean,1d-99)
205
endif
206
c do i=1,cur_it-1
207
do i=cur_it-3,cur_it-1
208
write(66,'(i4,4e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i)
209
enddo
210
close(66)
211
else
212
open(unit=66,file='results.dat',status='unknown')
213
write(66,'(3e12.5,2i9,i5,i9,e10.3)')0,0,0.0,kevent,nw,
214
& 1,0, 0
215
write(66,'(i4,4e15.5)') 1,0,0,0,0
216
close(66)
217
218
endif
219
c
220
c Now let's check to see if we got all of the events we needed
221
c if not, will give it another try with 5 iterations to set
222
c the grid, and 4 more to try and get the appropriate number of
223
c unweighted events.
224
c
225
write(*,*) "Status",accur, cur_it, itmax
226
if (accur .ge. 0d0 .or. cur_it .gt. itmax+3) return
227
228
if (neventswritten .gt. -accur) then
229
write(*,*) "We found enough events",neventswritten, -accur*1000*tmean
230
return
231
endif
232
233
c
234
c Need to start from scratch. This is clunky but I'll just
235
c remove the grid, so we are clean
236
c
237
write(*,*) "Trying w/ fresh grid"
238
open(unit=25,file='ftn25',status='unknown',err=102)
239
write(25,*) ' '
240
102 close(25)
241
242
c
243
c First few iterations will allow the grid to adjust
244
c
245
c
246
c Reset counters
247
c
248
ievent = 0
249
kevent = 0
250
nzoom = 0
251
xzoomfact = 1d0
252
253
ncall = ncall*4 ! / 2**(itmax-2)
254
write(*,*) "Starting w/ ncall = ", ncall
255
itmax = 8
256
call sample_init(ndim,ncall,itmax,ninvar,nconfigs)
257
do i=1,itmax
258
xmean(i)=0d0
259
xsigma(i)=0d0
260
enddo
261
wgt = 0d0
262
call clear_events
263
call set_peaks
264
c
265
c Main Integration Loop
266
c
267
iter = 1
268
c itmax = 8
269
itmax_adjust = 5
270
use_cut = 2 !Start adjusting grid
271
do while(iter .le. itmax)
272
if (iter .gt. itmax_adjust .and. use_cut .ne. 0) then
273
use_cut=0 !Fix grid
274
write(*,*) 'Fixing grid'
275
endif
276
c
277
c Get integration point
278
c
279
call sample_get_config(wgt,iter,ipole)
280
if (iter .le. itmax) then
281
ievent=ievent+1
282
call x_to_f_arg(ndim,ipole,mincfig,maxcfig,ninvar,wgt,x,p)
283
if (pass_point(p)) then
284
xzoomfact = 1d0
285
fx = dsig(p,wgt) !Evaluate function
286
if (xzoomfact .gt. 0d0) then
287
wgt = wgt*fx*xzoomfact
288
else
289
wgt = -xzoomfact
290
endif
291
if (wgt .gt. 0d0) call graph_point(p,wgt) !Update graphs
292
else
293
fx =0d0
294
wgt=0d0
295
endif
296
if (nzoom .le. 0) then
297
call sample_put_point(wgt,x(1),iter,ipole) !Store result
298
else
299
nzoom = nzoom -1
300
ievent=ievent-1
301
endif
302
endif
303
if (wgt .gt. 0d0) kevent=kevent+1
304
199 enddo
305
c
306
c All done
307
c
308
open(unit=66,file='results.dat',status='unknown')
309
i=1
310
do while(xmean(i) .ne. 0 .and. i .lt. cur_it)
311
i=i+1
312
enddo
313
cur_it = i
314
i = cur_it - 3
315
if (i .gt. 0) then
316
tmean = 0d0
317
tsigma = 0d0
318
tdem = 0d0
319
do while (xmean(i) .ne. 0 .and. i .lt. cur_it)
320
tmean = tmean+xmean(i)*xmean(i)**2/xsigma(i)**2
321
tdem = tdem+xmean(i)**2/xsigma(i)**2
322
tsigma = tsigma + xmean(i)**2/ xsigma(i)**2
323
i=i+1
324
enddo
325
tmean = tmean/tsigma
326
tsigma= tmean/sqrt(tsigma)
327
c nun = n_unwgted()
328
c
329
c tjs 8/7/2007
330
c
331
nun = neventswritten
332
333
chi2 = 0d0
334
do i = cur_it-3,cur_it-1
335
chi2 = chi2+(xmean(i)-tmean)**2/xsigma(i)**2
336
enddo
337
chi2 = chi2/2d0 !Since using only last 3, n-1=2
338
write(*,'(a)') '-----------------------------------------------------'
339
write(*,'(a)') '---------------------------'
340
write(*,'(a,e12.4)') ' Results Last 3 iters: Integral = ',tmean
341
write(*,'(25x,a,e12.4)') 'Std dev = ',tsigma
342
write(*,'(17x,a,f12.4)') 'Chi**2 per DoF. =',chi2
343
write(*,'(a)') '-----------------------------------------------------'
344
write(*,'(a)') '---------------------------'
345
346
if (nun .lt. 0) nun=-nun !Case when wrote maximun number allowed
347
if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
348
if (icor .eq. 0) then
349
write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,tsigma,0.0,kevent,nw,
350
& cur_it-1,nun, nun/max(tmean,1d-99)
351
else
352
write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,0.0,tsigma,kevent,nw,
353
& cur_it-1,nun, nun/max(tmean,1d-99)
354
endif
355
c do i=1,cur_it-1
356
do i=cur_it-3,cur_it-1
357
write(66,'(i4,4e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i)
358
enddo
359
close(66)
360
else
361
open(unit=66,file='results.dat',status='unknown')
362
write(66,'(3e12.5,2i9,i5,i9,e10.3)')0,0,0.0,kevent,nw,
363
& 1,0, 0
364
write(66,'(i4,4e15.5)') 1,0,0,0,0
365
close(66)
366
367
endif
368
end
369
370
subroutine sample_fullx(ndim,ncall,itmax,dsig,ninvar,nconfigs)
371
c**************************************************************************
372
c Driver for sample which does complete integration
373
c This is done in double precision, and should be told the
374
c number of possible phasespace choices.
375
c Arguments:
376
c ndim Number of dimensions for integral(number or random #'s/point)
377
c ncall Number of times to evaluate the function/iteration
378
c itmax Number of iterations
379
c ninvar Number of invarients to keep grids on (s,t,u, s',t' etc)
380
c nconfigs Number of different pole configurations
381
c dsig Function to be integrated
382
c**************************************************************************
383
implicit none
384
include 'genps.inc'
385
c
386
c Arguments
387
c
388
integer ndim,ncall,itmax,ninvar,nconfigs
389
external dsig
390
double precision dsig
391
c
392
c Local
393
c
394
double precision x(maxinvar),wgt,p(4*maxdim/3+14)
395
double precision tdem, chi2
396
integer ievent,kevent,nwrite,iter,nun
397
integer jmax,i,j,ipole
398
c
399
c External
400
c
401
integer n_unwgted
402
external n_unwgted
403
c
404
c Global
405
c
406
integer nsteps
407
character*40 result_file,where_file
408
common /sample_status/result_file,where_file,nsteps
409
double precision fx
410
common /to_fx/ fx
411
412
integer mincfig, maxcfig
413
common/to_configs/mincfig, maxcfig
414
415
double precision xmean(99),xsigma(99),xwmax(99),xeff(99)
416
common/to_iterations/xmean, xsigma, xwmax, xeff
417
418
double precision accur
419
common /to_accuracy/accur
420
421
double precision twgt, maxwgt,swgt(maxevents)
422
integer lun, nw
423
common/to_unwgt/twgt, maxwgt, swgt, lun, nw
424
425
integer nzoom
426
double precision tx(1:3,maxinvar)
427
common/to_xpoints/tx, nzoom
428
429
double precision xzoomfact
430
common/to_zoom/ xzoomfact
431
432
double precision tmean, tsigma
433
integer dim, events, itm, kn, cur_it, invar, configs
434
common /sample_common/
435
. tmean, tsigma, dim, events, itm, kn, cur_it, invar, configs
436
437
integer icor
438
common/to_correlated/icor
439
c
440
c External
441
c
442
logical pass_point
443
c
444
c Data
445
c
446
c data result_file,where_file,nsteps/'SAMPLE','WHERE.AMI',100/
447
c data accur/-1d0/
448
c data mincfig /1/
449
c data maxcfig /1/
450
c data twgt/-1d0/ !Dont write out events
451
c data lun/27/ !Unit number for events
452
c data maxwgt/0d0/
453
c data nw/0/ !Number of events written
454
455
456
c-----
457
c Begin Code
458
c-----
459
ievent = 0
460
kevent = 0
461
nzoom = 0
462
xzoomfact = 1d0
463
if (nsteps .lt. 1) nsteps=1
464
nwrite = itmax*ncall/nsteps
465
c open(unit=66,file='.sample_warn',status='unknown')
466
c write(66,*) 'Warnings from sample run.',itmax,ncall
467
c close(66)
468
call sample_init(ndim,ncall,itmax,ninvar,nconfigs)
469
call graph_init
470
do i=1,itmax
471
xmean(i)=0d0
472
xsigma(i)=0d0
473
enddo
474
c mincfig=1
475
c maxcfig=nconfigs
476
wgt = 0d0
477
c
478
c Main Integration Loop
479
c
480
iter = 1
481
do while(iter .le. itmax)
482
c
483
c Get integration point
484
c
485
call sample_get_config(wgt,iter,ipole)
486
if (iter .le. itmax) then
487
ievent=ievent+1
488
call x_to_f_arg(ndim,ipole,mincfig,maxcfig,ninvar,wgt,x,p)
489
if (pass_point(p)) then
490
xzoomfact = 1d0
491
fx = dsig(p,wgt) !Evaluate function
492
if (xzoomfact .gt. 0d0) then
493
wgt = wgt*fx*xzoomfact
494
else
495
wgt = -xzoomfact
496
endif
497
if (wgt .gt. 0d0) call graph_point(p,wgt) !Update graphs
498
else
499
fx =0d0
500
wgt=0d0
501
endif
502
if (nzoom .le. 0) then
503
call sample_put_point(wgt,x(1),iter,ipole) !Store result
504
else
505
nzoom = nzoom -1
506
ievent=ievent-1
507
endif
508
endif
509
if (wgt .gt. 0d0) kevent=kevent+1
510
c
511
c Write out progress/histograms
512
c
513
if (kevent .ge. nwrite) then
514
nwrite = nwrite+ncall*itmax/nsteps
515
nwrite = min(nwrite,ncall*itmax)
516
c open(unit=22,file=where_file,status='old',
517
c & access='append',err=99)
518
c write(22,'(2i15)') ievent,kevent
519
c close(22)
520
call graph_store
521
endif
522
99 enddo
523
c
524
c All done
525
c
526
c open(unit=66,file='.sample_warn',status='old',access='append')
527
c write(66,*) 'Finished sample',ievent,kevent,wgt
528
c if (wgt .ne. -1) then
529
c do j=1,5
530
c jmax = min(ndim,4*j)
531
c write(66,'(4e19.12)') (x(i),i=(j-1)*4+1,jmax)
532
c enddo
533
c endif
534
c close(66)
535
536
open(unit=66,file='results.dat',status='unknown')
537
i=1
538
do while(xmean(i) .ne. 0 .and. i .lt. cur_it)
539
i=i+1
540
enddo
541
cur_it = i
542
i = cur_it - 3
543
if (i .gt. 0) then
544
tmean = 0d0
545
tsigma = 0d0
546
tdem = 0d0
547
do while (xmean(i) .ne. 0 .and. i .lt. cur_it)
548
tmean = tmean+xmean(i)*xmean(i)**2/xsigma(i)**2
549
tdem = tdem+xmean(i)**2/xsigma(i)**2
550
tsigma = tsigma + xmean(i)**2/ xsigma(i)**2
551
i=i+1
552
enddo
553
c tmean = tmean/dble(i-1)
554
c tsigma= sqrt(tsigma)/dble(i-1)
555
tmean = tmean/tsigma
556
c tsigma= sqrt(tsigma)/dble(3)
557
tsigma= tmean/sqrt(tsigma)
558
nun = n_unwgted()
559
560
chi2 = 0d0
561
do i = cur_it-3,cur_it-1
562
chi2 = chi2+(xmean(i)-tmean)**2/xsigma(i)**2
563
enddo
564
chi2 = chi2/2d0 !Since using only last 3, n-1=2
565
c tsigma = tsigma*sqrt(chi2)
566
c write(*,*) "chi2 / dof=", chi2
567
write(*,'(a)') '-----------------------------------------------------'
568
write(*,'(a)') '---------------------------'
569
write(*,'(a,e12.4)') ' Results Last 3 iters: Integral = ',tmean
570
write(*,'(25x,a,e12.4)') 'Std dev = ',tsigma
571
write(*,'(17x,a,f12.4)') 'Chi**2 per DoF. =',chi2
572
write(*,'(a)') '-----------------------------------------------------'
573
write(*,'(a)') '---------------------------'
574
575
if (nun .lt. 0) nun=-nun !Case when wrote maximun number allowed
576
if (chi2 .gt. 1) tsigma=tsigma*sqrt(chi2)
577
if (icor .eq. 0) then
578
write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,tsigma,0.0,kevent,nw,
579
& cur_it-1,nun, nun/max(tmean,1d-99)
580
else
581
write(66,'(3e12.5,2i9,i5,i9,e10.3)')tmean,0.0,tsigma,kevent,nw,
582
& cur_it-1,nun, nun/max(tmean,1d-99)
583
endif
584
c do i=1,cur_it-1
585
do i=cur_it-3,cur_it-1
586
write(66,'(i4,4e15.5)') i,xmean(i),xsigma(i),xeff(i),xwmax(i)
587
enddo
588
close(66)
589
else
590
open(unit=66,file='results.dat',status='unknown')
591
write(66,'(3e12.5,2i9,i5,i9,e10.3)')0,0,0.0,kevent,nw,
592
& 1,0, 0
593
write(66,'(i4,4e15.5)') 1,0,0,0,0
594
close(66)
595
596
endif
597
598
end
599
600
601
subroutine sample_writehtm()
602
c***********************************************************************
603
c Writes out results of run in html format
604
c***********************************************************************
605
implicit none
606
c
607
c Constants
608
c
609
character*(*) htmfile
610
parameter (htmfile='results.html')
611
integer lun
612
parameter (lun=26)
613
c
614
c Local
615
c
616
character*4 cpref
617
double precision scale
618
integer i
619
c
620
c Global
621
c
622
double precision xmean(99),xsigma(99),xwmax(99),xeff(99)
623
common/to_iterations/xmean, xsigma, xwmax, xeff
624
625
c-----
626
c Begin Code
627
c-----
628
return
629
c
630
c Here we determine the appropriate units. Assuming the results
631
c were written in picobarns
632
c
633
if (xmean(1) .ge. 1e4) then !Use nano barns
634
scale=1d-3
635
cpref='(nb)'
636
elseif (xmean(1) .ge. 1e1) then !Use pico barns
637
scale=1d0
638
cpref='(pb)'
639
else !Use fempto
640
scale=1d+3
641
cpref='(fb)'
642
endif
643
open(unit=lun,file=htmfile,status='unknown',err=999)
644
write(lun,50) '<head><title>Results_head</title></head>'
645
write(lun,50) '<body><h2>Results for Process</h2>'
646
write(lun,50) '<table border>'
647
write(lun,50) '<Caption> Caption Results'
648
write(lun,49) '<tr><th>Iteration</th>'
649
write(lun,48)'<th>Cross Sect',cpref,'</th><th>Error',cpref,'</th>'
650
write(lun,49) '<th>Events (K)</th><th>Eff</th>'
651
write(lun,50) '<th>Wrote</th><th>Unwgt</th></tr>'
652
653
c write(lun,60) '<tr><th>AVG</th><th>',xtot*scale
654
c $ ,'</th><th>',errtot*scale,'</th><th align=right>',
655
c $ ntot/1000,'</th><th align=right>',teff,'</th></tr>'
656
i=1
657
do while(xmean(i) .gt. 0d0)
658
write(lun,'(a)') '<tr>'
659
write(lun,45) '<td align=right>',i,'</tr>'
660
write(lun,46) '<td align=right>',xmean(i)*scale,'</td>'
661
write(lun,46) '<td align=right>',xsigma(i)*scale,'</td>'
662
write(lun,46) '<td align=right>',xeff(i)*scale,'</td>'
663
write(lun,'(a)') '</tr>'
664
i=i+1
665
enddo
666
write(lun,50) '</table></body>'
667
999 close(lun)
668
45 format(a,i4,a)
669
46 format(a,f12.3,a)
670
48 format(a,a,a,a)
671
49 format(a)
672
50 format(a)
673
end
674
675
676
677
subroutine sample_init(p1, p2, p3, p4, p5)
678
c************************************************************************
679
c Initialize grid and random number generators
680
c************************************************************************
681
implicit none
682
c
683
c Constants
684
c
685
include 'genps.inc'
686
c
687
c Arguments
688
c
689
integer p1, p2, p3, p4, p5
690
691
c
692
c Local
693
c
694
integer i, j
695
c
696
c Global
697
c
698
integer nsteps
699
character*40 result_file,where_file
700
common /sample_status/result_file,where_file,nsteps
701
702
double precision tmean, tsigma
703
integer dim, events, iter, kn, cur_it, invar, configs
704
common /sample_common/
705
. tmean, tsigma, dim, events, iter, kn, cur_it, invar, configs
706
707
double precision grid(2, ng, 0:maxinvar)
708
common /data_grid/ grid
709
integer Minvar(maxdim,lmaxconfigs)
710
common /to_invar/ Minvar
711
double precision psect(maxconfigs),alpha(maxconfigs)
712
common/to_mconfig2/psect ,alpha
713
logical first_time
714
common/to_first/first_time
715
integer use_cut
716
common /to_weight/use_cut
717
integer ituple
718
common /to_random/ituple
719
720
logical flat_grid
721
common/to_readgrid/flat_grid !Tells if grid read from file
722
723
integer icor
724
common/to_correlated/icor
725
726
logical zooming
727
common /to_zoomchoice/zooming
728
729
data use_cut/2/ !Grid: 0=fixed , 1=standard, 2=non-zero
730
data ituple/1/ !1=htuple, 2=sobel
731
data Minvar(1,1)/-1/ !No special variable mapping
732
733
c-----
734
c Begin Code
735
c-----
736
icor = 0
737
If (use_cut .eq. 0) then
738
icor = 1 !Assume correlated unless grid read
739
print*,'Keeping grid fixed.'
740
elseif(use_cut .eq. 1) then
741
print*,'Using standard SAMPLE grid deformation.'
742
elseif(use_cut .eq. 2) then
743
print*,'Using non-zero grid deformation.'
744
elseif(use_cut .eq. 3) then
745
print*,'Using fluctuation for grid deformation.'
746
elseif(use_cut .eq. 4) then
747
print*,'Generating unweighted event shape.'
748
elseif(use_cut .eq. 5) then
749
print*,'Using constant plus linear grid deformation.'
750
elseif(use_cut .eq. 6) then
751
print*,'Using power law grid deformation.'
752
else
753
print*,'Using unknown grid deformation:',use_cut
754
endif
755
c open(unit=22,file=result_file,status='unknown')
756
c write(22,*) 'Sample Status ',p2,p3,nsteps
757
c close(22)
758
c open(unit=22,file=where_file,status='unknown')
759
c write(22,*) 'Sample Progress ',p2,p3,nsteps
760
c close(22)
761
762
dim = p1
763
events = p2
764
iter = p3
765
invar = p4
766
configs = p5
767
first_time = .true.
768
769
if (dim .gt. maxdim) then
770
write(*,*) 'Too many dimensions requested from Sample()'
771
stop
772
endif
773
c if (dim .gt. invar) then
774
c write(*,*) 'Too many dimensions dim > invar',dim,invar
775
c stop
776
c endif
777
if (p4 .gt. maxinvar) then
778
write(*,*) 'Too many invarients requested from Sample()',p4
779
stop
780
endif
781
if (p5 .gt. maxconfigs) then
782
write(*,*) 'Too many configs requested from Sample()',p5
783
stop
784
endif
785
786
write(*,'(i3,a,i7,a,i3,a,i3,a,i3,a)') dim, ' dimensions', events,
787
& ' events',p4,' invarients',iter, ' iterations',
788
& p5,' config(s), (0.99)'
789
790
if (ituple .eq. 1) then
791
print*,'Using h-tuple random number sequence.'
792
elseif (ituple .eq. 2) then
793
print*,'Using Sobel quasi-random number sequence.'
794
write(*,*) 'Sorry cant use sobel'
795
stop
796
c call isobel(dim)
797
else
798
print*,'Unknown random number generator',ituple
799
endif
800
c
801
c See if need mapping between dimensions in different configurations
802
c (ie using s,t,u type invarients)
803
c
804
if (Minvar(1,1) .eq. -1) then
805
print*,'No invarient mapping defined, using 1 to 1.'
806
do i=1,configs
807
do j=1,dim
808
Minvar(j,i) = j+(i-1)*dim
809
enddo
810
enddo
811
endif
812
c
813
c Reset counters
814
c
815
tmean = 0d0
816
tsigma = 0d0
817
kn = 0
818
cur_it = 1
819
do j=1,ng
820
grid(2,j,0) = xgmin+(xgmax-xgmin)*j/dble(ng)
821
enddo
822
c
823
c Try to read grid from file
824
c
825
flat_grid=.true.
826
open(unit=25,file='ftn25',status='unknown',err=102)
827
read(25,fmt='(4f20.17)', err=101, end=101)
828
. ((grid(2,i,j),i=1,ng),j=1,maxinvar)
829
c write(*,*) 'Got the grid OK, now getting alpha'
830
c read(25,fmt='(4f20.17)',err=101,end=101)(alpha(i),i=1,maxconfigs)
831
write(*,*) 'Grid read from file'
832
flat_grid=.false.
833
close(25)
834
c
835
c Determine weighting for each configuration
836
c
837
if (.not. flat_grid) icor = 0 !0 = not correlated
838
zooming = (.not. flat_grid .and. use_cut .eq. 0) !only zoom if grid already adjusted and not changing more
839
c
840
c tjs 5/22/07 turn off zooming
841
c
842
zooming = .false.
843
if (configs .eq. 1) then
844
do i=1,maxconfigs
845
alpha(i) = 1
846
enddo
847
else
848
write(*,*) 'Using uniform alpha',alpha(1)
849
c tot=0d0
850
c do i=1,configs
851
c tot=tot+alpha(i)
852
c enddo
853
do i=1,maxconfigs
854
if(i .le. configs) then
855
alpha(i)=1d0/dble(configs)
856
else
857
alpha(i)=0d0
858
endif
859
enddo
860
endif
861
goto 103
862
101 close(25)
863
c write(*,*) 'Tried reading it',i,j
864
102 write(*,*) 'Error opening grid'
865
c
866
c Unable to read grid, using uniform grid and equal points in
867
c each configuration
868
c
869
write(*,*) 'Using Uniform Grid!'
870
do j = 1, maxinvar
871
do i = 1, ng
872
grid(2, i, j) = xgmin+ (xgmax-xgmin)*(i / dble(ng))**1
873
end do
874
end do
875
do j=1,maxconfigs
876
if (j .le. configs) then
877
alpha(j)=1d0/dble(configs)
878
else
879
alpha(j)=0d0
880
endif
881
enddo
882
write(*,*) 'Using uniform alpha',alpha(1)
883
c write(*,*) 'Forwarding random number generator'
884
885
103 write(*,*) 'Grid defined OK'
886
end
887
888
subroutine setgrid(j,xo,a,itype)
889
c*************************************************************************
890
c Presets the grid for a 1/(x-a)^itype distribution down to xo
891
c*************************************************************************
892
implicit none
893
c
894
c Constants
895
c
896
include 'genps.inc'
897
c
898
c Arguments
899
c
900
integer j, itype !grid number
901
double precision xo !minimum value
902
double precision a !offset for peak
903
c
904
c Local
905
c
906
integer i,k
907
integer ngu, ngd
908
c
909
c Global
910
c
911
double precision grid(2, ng, 0:maxinvar)
912
common /data_grid/ grid
913
914
logical flat_grid
915
common/to_readgrid/flat_grid !Tells if grid read from file
916
917
c-----
918
c Begin Code
919
c-----
920
if (flat_grid) then
921
write(*,'(a,i4,2e15.5,i4)') 'Setting grid',j,xo,a,itype
922
if (a .ge. xo) then
923
write(*,*) 'Can not integrate over singularity'
924
write(*,*) 'Set grid',j,xo,a
925
return
926
endif
927
c grid(2,1,j) = xo
928
grid(2,ng,j)=xgmax
929
if (itype .eq. 1) then
930
c
931
c We'll use most for the peak, but save some for going down
932
c
933
ngu = ng *0.9
934
ngd = ng-ngu
935
do i=1,ngu-1
936
c-------------------
937
c tjs 6/30/2009
938
c New form for setgrid
939
c-------------------
940
c grid(2,i+ngd,j)=((1d0-a)/(xo-a))**(1d0-dble(i)/dble(ngu))
941
c grid(2,i+ngd,j)=1d0/grid(2,i+ngd,j)+a
942
grid(2,i+ngd,j) = xo + ((dble(i)+xo-a)/(dble(ngu)+xo-a))**2
943
944
enddo
945
c
946
c Now lets go down the other side
947
c
948
grid(2,ngd,j) = xo
949
do i=1,ngd-1
950
c grid(2,i,j) = ((1d0-a)/(xo-a))**(1d0-dble(i)/dble(ngd))
951
grid(2,ngd-i,j) = xo-(grid(2,ngd+i,j)-xo)
952
if (grid(2,ngd-i,j) .lt. -1d0) then
953
write(*,*) 'Error grid set too low',grid(2,ngd-i,j)
954
do k=1,ng
955
write(*,*) k,grid(2,k,j)
956
enddo
957
stop
958
endif
959
enddo
960
c
961
c tjs 5/11/2009
962
c Make sure sample all the way down to zero
963
c
964
if (xo .gt. 0) grid(2,1,j) = 0d0
965
c write(*,*) "Adjusted bin 1 to zero"
966
967
elseif (itype .eq. 2) then
968
do i=2,ng-1
969
grid(2,i,j)=(1d0/(xo-a))*(1d0-dble(i)/dble(ng))+
970
$ (dble(i)/dble(ng))*(1d0/(1d0-a))
971
grid(2,i,j)=1d0/grid(2,i,j)+a
972
enddo
973
else
974
write(*,*) 'No modification in setgrid',itype
975
endif
976
do i=1,ng
977
c write(*,*) j,i,grid(2,i,j)
978
enddo
979
call sample_write_g(j,'_0')
980
else
981
write(*,*) 'No modification is setgrid, grid read from file'
982
endif
983
end
984
985
subroutine sample_get_config(wgt, iteration, iconfig)
986
c************************************************************************
987
c
988
c INPUTS:
989
c
990
c OUTPUTS: wgt == 1/nevents*niterations
991
c iteration == Current iteration
992
c iconfig == configuration to use
993
c
994
c************************************************************************
995
implicit none
996
c
997
c Constants
998
c
999
include 'genps.inc'
1000
c
1001
c Arguments
1002
c
1003
double precision wgt
1004
integer iteration, iconfig
1005
c
1006
c Local
1007
c
1008
integer idum
1009
real xrnd
1010
double precision tot
1011
c
1012
c External
1013
c
1014
real ran1
1015
c
1016
c Global
1017
c
1018
double precision tmean, tsigma
1019
integer dim, events, iter, kn, cur_it, invar, configs
1020
common /sample_common/
1021
. tmean, tsigma, dim, events, iter, kn, cur_it, invar, configs
1022
double precision psect(maxconfigs),alpha(maxconfigs)
1023
common/to_mconfig2/psect ,alpha
1024
data idum/0/
1025
1026
integer mincfig, maxcfig
1027
common/to_configs/mincfig, maxcfig
1028
1029
c-----
1030
c Begin Code
1031
c-----
1032
iteration = cur_it
1033
if (cur_it .gt. iter) then
1034
wgt = -1d0
1035
else
1036
wgt = 1d0 / (dble(events) * dble(iter))
1037
c
1038
c Choose configuration
1039
c
1040
if (configs .gt. 1) then
1041
xrnd = ran1(idum)
1042
iconfig=1
1043
tot = alpha(iconfig)
1044
do while (tot .lt. xrnd .and. iconfig .lt. configs)
1045
iconfig=iconfig+1
1046
tot = tot+alpha(iconfig)
1047
enddo
1048
else
1049
iconfig=mincfig
1050
endif
1051
endif
1052
end
1053
1054
subroutine sample_get_x(wgt, x, j, ipole, xmin, xmax)
1055
c************************************************************************
1056
c Returns maxdim random numbers between 0 and 1, and the wgt
1057
c associated with this set of points, and the iteration number
1058
c This routine chooses the point within the range specified by
1059
c xmin and xmax for dimension j in configuration ipole
1060
c************************************************************************
1061
implicit none
1062
c
1063
c Constants
1064
c
1065
include 'genps.inc'
1066
c
1067
c Arguments
1068
c
1069
double precision wgt, x, xmin, xmax
1070
integer j, ipole
1071
c
1072
c Local
1073
c
1074
integer im, ip,ij,icount,it_warned
1075
double precision xbin_min,xbin_max,ddum(maxdim),xo,y
1076
c
1077
c External
1078
c
1079
double precision xbin
1080
external xbin
1081
c
1082
c Global
1083
c
1084
double precision tmean, tsigma
1085
integer dim, events, iter, kn, cur_it, invar, configs
1086
common /sample_common/
1087
. tmean, tsigma, dim, events, iter, kn, cur_it, invar, configs
1088
1089
double precision grid(2, ng, 0:maxinvar)
1090
common /data_grid/ grid
1091
integer Minvar(maxdim,lmaxconfigs)
1092
common /to_invar/ Minvar
1093
1094
integer ituple
1095
common /to_random/ituple
1096
1097
double precision spole(maxinvar),swidth(maxinvar),bwjac
1098
common/to_brietwigner/spole ,swidth ,bwjac
1099
1100
integer nzoom
1101
double precision tx(1:3,maxinvar)
1102
common/to_xpoints/tx, nzoom
1103
1104
data ddum/maxdim*0d0/
1105
data icount/0/
1106
data it_warned/0/
1107
c-----
1108
c Begin Code
1109
c-----
1110
if (it_warned .ne. cur_it) then
1111
icount=0
1112
it_warned = cur_it
1113
endif
1114
if (ituple .eq. 2) then !Sobel generator
1115
print*,'Sorry Sobel generator disabled'
1116
stop
1117
c call sobel(ddum)
1118
c write(*,'(7f11.5)')(ddum(j)*real(ng),j=1,dim)
1119
endif
1120
if (ituple .eq. 1) then
1121
c write(*,*) 'Getting variable',ipole,j,minvar(j,ipole)
1122
xbin_min = xbin(xmin,minvar(j,ipole))
1123
xbin_max = xbin(xmax,minvar(j,ipole))
1124
if (xbin_min .gt. xbin_max) then
1125
c write(*,'(a,4e15.4)') 'Bad limits',xbin_min,xbin_max,
1126
c & xmin,xmax
1127
c xbin_max=xbin_min+1d-10
1128
xbin_min = xbin(xmin,minvar(j,ipole))
1129
xbin_max = xbin(xmax,minvar(j,ipole))
1130
endif
1131
c
1132
c Line which allows us to keep choosing same x
1133
c
1134
c if (swidth(j) .ge. 0) then
1135
if (nzoom .le. 0) then
1136
call ntuple(ddum(j), xbin_min,xbin_max, j, ipole)
1137
else
1138
c write(*,*) 'Reusing num',j,nzoom,tx(2,j)
1139
1140
call ntuple(ddum(j),max(xbin_min,dble(int(tx(2,j)))),
1141
$ min(xbin_max,dble(int(tx(2,j))+1)),j,ipole)
1142
1143
if(max(xbin_min,dble(int(tx(2,j)))).gt.
1144
$ min(xbin_max,dble(int(tx(2,j))+1))) then
1145
c write(*,*) 'not good'
1146
endif
1147
1148
c write(*,'(2i6,4e15.5)') nzoom,j,ddum(j),tx(2,j),
1149
c $ max(xbin_min,dble(int(tx(2,j)))),
1150
c $ min(xbin_max,dble(int(tx(2,j))+1))
1151
1152
c ddum(j) = tx(2,j) !Use last value
1153
1154
1155
endif
1156
tx(1,j) = xbin_min
1157
tx(2,j) = ddum(j)
1158
tx(3,j) = xbin_max
1159
elseif (ituple .eq. 2) then
1160
if (ipole .gt. 1) then
1161
print*,'Sorry Sobel not configured for multi-pole.'
1162
stop
1163
endif
1164
ddum(j)=ddum(j)*dble(ng)
1165
else
1166
print*,'Error unknown random number generator.',ituple
1167
stop
1168
endif
1169
1170
im = ddum(j)
1171
ip = im + 1
1172
ij = Minvar(j,ipole)
1173
c
1174
c New method of choosing x from bins
1175
c
1176
if (ip .eq. 1) then !This is in the first bin
1177
xo = grid(2, ip, ij)-xgmin
1178
x = grid(2, ip, ij) - xo * (dble(ip) - ddum(j))
1179
else
1180
xo = grid(2, ip, ij)-grid(2,im,ij)
1181
x = grid(2, ip, ij) - xo * (dble(ip) - ddum(j))
1182
endif
1183
c
1184
c Now we transform x if there is a B.W., S, or T pole
1185
c
1186
if (ij .gt. 0) then
1187
if (swidth(ij) .gt. 0d0) then
1188
y = x !Takes uniform y and returns
1189
c write(*,*) 'Tranpole called',ij,swidth(ij)
1190
call transpole(spole(ij),swidth(ij),y,x,wgt) !x on BW pole
1191
endif
1192
endif
1193
c
1194
c Simple checks to see if we got the right point note 1e-3 corresponds
1195
c to the fact that the grids are required to be separated by 1e-14. Since
1196
c double precision is about 18 digits, we expect things to agree to
1197
c 3 digit accuracy.
1198
c
1199
if (abs(ddum(j)-xbin(x,ij))/(ddum(j)+1d-22) .gt. 1e-3) then
1200
if (icount .lt. 5) then
1201
write(*,'(a,i4,2e14.6,1e12.4)')
1202
& 'Warning xbin not returning correct x', ij,
1203
& ddum(j),xbin(x,ij),xo
1204
elseif (icount .eq. 5) then
1205
write(*,'(a,a)')'Warning xbin still not working well. ',
1206
& 'Last message this iteration.'
1207
endif
1208
icount=icount+1
1209
endif
1210
if (x .lt. xmin .or. x .gt. xmax) then
1211
c write(*,'(a,4i4,2f24.16,1e10.2)') 'Bad x',ij,int(xbin_min),ip,
1212
c & int(xbin_max),xmin,x,xmax-xmin
1213
endif
1214
1215
wgt = wgt * xo * dble(xbin_max-xbin_min)
1216
c print*,'Returning x',ij,ipole,j,x
1217
end
1218
1219
subroutine sample_get_wgt(wgt, x, j, ipole, xmin, xmax)
1220
c************************************************************************
1221
c Returns the wgt for a point x in grid j of configuration
1222
c ipole between xmin and xmax
1223
c************************************************************************
1224
implicit none
1225
c
1226
c Constants
1227
c
1228
include 'genps.inc'
1229
c
1230
c Arguments
1231
c
1232
double precision wgt, x, xmin, xmax
1233
integer j, ipole
1234
c
1235
c Local
1236
c
1237
integer im, ip,ij
1238
double precision xbin_min,xbin_max,xbin2
1239
double precision xo
1240
c
1241
c External
1242
c
1243
double precision xbin
1244
external xbin
1245
c
1246
c Global
1247
c
1248
double precision tmean, tsigma
1249
integer dim, events, iter, kn, cur_it, invar, configs
1250
common /sample_common/
1251
. tmean, tsigma, dim, events, iter, kn, cur_it, invar, configs
1252
1253
double precision grid(2, ng, 0:maxinvar)
1254
common /data_grid/ grid
1255
integer Minvar(maxdim,lmaxconfigs)
1256
common /to_invar/ Minvar
1257
integer ituple
1258
common /to_random/ituple
1259
double precision spole(maxinvar),swidth(maxinvar),bwjac
1260
common/to_brietwigner/spole ,swidth ,bwjac
1261
1262
c-----
1263
c Begin Code
1264
c-----
1265
if (xmin .gt. x) then
1266
if (xmin-x .lt. 1d-13) then
1267
x=xmin
1268
else
1269
write(*,'(a,2i4,4e10.4)') 'Error x out of range in get_wgt',
1270
$ j,minvar(j,ipole),xmin,x,xmax,x-xmin
1271
return
1272
endif
1273
endif
1274
if (xmax .lt. x) then
1275
if (x-xmax .lt. 1d-13) then
1276
x=xmax
1277
else
1278
write(*,'(a,2i4,4f8.4)') 'Error x out of range in get_wgt',
1279
$ j,minvar(j,ipole),xmin,x,xmax,x-xmin
1280
return
1281
endif
1282
endif
1283
if (ituple .eq. 1) then
1284
xbin_min = xbin(xmin,minvar(j,ipole))
1285
xbin_max = xbin(xmax,minvar(j,ipole))
1286
xbin2 = xbin(x,minvar(j,ipole)) !This must be last one for bwjac
1287
if (xbin_min .gt. xbin_max) then
1288
write(*,'(a,2e15.3,i6,2e15.3)') 'Error xbinmin>xbinmax'
1289
& ,xbin_min,
1290
& xbin_max,minvar(j,ipole),xmin,xmax
1291
endif
1292
else
1293
print*,'Error unknown random number generator.',ituple
1294
stop
1295
endif
1296
im = xbin2
1297
ip = im + 1
1298
ij = Minvar(j,ipole)
1299
c
1300
c New method for finding bin
1301
c
1302
if (ip .eq. 1) then
1303
xo=grid(2,ip,ij)-xgmin
1304
else
1305
xo=grid(2,ip,ij)-grid(2,im,ij)
1306
endif
1307
wgt = wgt * xo * dble(xbin_max-xbin_min)*bwjac
1308
if (wgt .le. 0d0) then
1309
c write(*,'(a,3i4,2f6.1,3e15.3)') 'Error wgt<0',j,ij,ip,
1310
c & xbin_min,xbin_max,xo,xmin,xmax
1311
c write(*,'(2e25.15)') grid(2, ip, ij),grid(2, im, ij)
1312
c write(*,'(a,5e15.5)') 'Wgt',wgt,xo,
1313
c & dble(xbin_max-xbin_min),bwjac
1314
endif
1315
end
1316
1317
subroutine sample_result(mean, sigma)
1318
implicit none
1319
double precision mean, sigma
1320
integer i,cur_it
1321
double precision tsigma,tmean,tsig,tdem
1322
1323
double precision xmean(99),xsigma(99),xwmax(99),xeff(99)
1324
common/to_iterations/xmean, xsigma, xwmax, xeff
1325
1326
1327
i=1
1328
do while(xmean(i) .ne. 0 .and. i .lt. 99)
1329
i=i+1
1330
enddo
1331
cur_it = i
1332
i = cur_it - 3
1333
tmean = 0d0
1334
tsigma = 0d0
1335
if (i .gt. 0) then
1336
tdem = 0d0
1337
do while (xmean(i) .ne. 0 .and. i .lt. cur_it)
1338
tmean = tmean+xmean(i)*xmean(i)**2/xsigma(i)**2
1339
tdem = tdem+xmean(i)**2/xsigma(i)**2
1340
tsigma = tsigma + xmean(i)**2/ xsigma(i)**2
1341
i=i+1
1342
enddo
1343
tmean = tmean/tsigma
1344
tsigma= tmean/sqrt(tsigma)
1345
endif
1346
1347
mean = tmean
1348
sigma = tsigma
1349
1350
end
1351
1352
1353
subroutine sample_put_point(wgt, point, iteration,ipole)
1354
c**************************************************************************
1355
c Given point(maxinvar),wgt and iteration, updates the grid.
1356
c If at the end of an iteration, reforms the grid as necessary
1357
c and outputs current results
1358
c**************************************************************************
1359
implicit none
1360
c
1361
c Constants
1362
c
1363
include 'genps.inc'
1364
integer max_events
1365
parameter (max_events=500000) !Maximum # events before get non_zero
1366
c
1367
c Arguments
1368
c
1369
integer iteration,ipole
1370
double precision wgt, point(maxinvar)
1371
c
1372
c Local
1373
c
1374
integer i, j, k, knt, non_zero, nun
1375
double precision vol,xnmin,xnmax,tot
1376
double precision rc, dr, xo, xn, x(maxinvar), dum(ng)
1377
save vol,knt
1378
double precision chi2
1379
save chi2, non_zero
1380
double precision wmax1,ddumb
1381
save wmax1
1382
double precision twgt1,xchi2,xmean
1383
integer iavg,navg
1384
save twgt1,iavg,navg
1385
c
1386
c External
1387
c
1388
double precision binwidth,xbin
1389
logical rebin
1390
integer n_unwgted
1391
external binwidth,xbin,rebin,n_unwgted
1392
c
1393
c Global
1394
c
1395
double precision accur
1396
common /to_accuracy/accur
1397
1398
double precision ymean(99),ysigma(99),ywmax(99),yeff(99)
1399
common/to_iterations/ymean, ysigma, ywmax, yeff
1400
1401
double precision mean,sigma
1402
common/to_result/mean,sigma
1403
1404
double precision grid2(0:ng,maxinvar)
1405
integer inon_zero(ng,maxinvar)
1406
common/to_grid2/grid2,inon_zero
1407
1408
double precision tmean, tsigma
1409
integer dim, events, iter, kn, cur_it, invar, configs
1410
common /sample_common/
1411
. tmean, tsigma, dim, events, iter, kn, cur_it, invar, configs
1412
1413
double precision grid(2, ng, 0:maxinvar)
1414
common /data_grid/ grid
1415
integer nsteps
1416
character*40 result_file,where_file
1417
common /sample_status/result_file,where_file,nsteps
1418
logical first_time
1419
common/to_first/first_time
1420
integer use_cut
1421
common /to_weight/use_cut
1422
double precision xmin(maxinvar),xmax(maxinvar)
1423
common /to_extreme/xmin ,xmax
1424
double precision reliable(ng,maxdim)
1425
common /to_error/reliable
1426
1427
double precision twgt, maxwgt,swgt(maxevents)
1428
integer lun, nw
1429
common/to_unwgt/twgt, maxwgt, swgt, lun, nw
1430
1431
1432
real*8 wmax !This is redundant
1433
common/to_unweight/wmax
1434
1435
double precision fx
1436
common /to_fx/ fx
1437
double precision prb(maxconfigs,maxpoints,maxplace)
1438
double precision fprb(maxinvar,maxpoints,maxplace)
1439
integer jpnt,jplace
1440
common/to_mconfig1/prb ,fprb,jpnt,jplace
1441
double precision psect(maxconfigs),alpha(maxconfigs)
1442
common/to_mconfig2/psect ,alpha
1443
double precision spole(maxinvar),swidth(maxinvar),bwjac
1444
common/to_brietwigner/spole ,swidth ,bwjac
1445
1446
integer neventswritten
1447
common /to_eventswritten/ neventswritten
1448
1449
data prb/maxprb*1d0/
1450
data fprb/maxfprb*1d0/
1451
data jpnt,jplace /1,1/
1452
1453
c-----
1454
c Begin Code
1455
c-----
1456
if (first_time) then
1457
first_time = .false.
1458
twgt1 = 0d0 !
1459
iavg = 0 !Vars for averging to increase err estimate
1460
navg = 1 !
1461
wmax1= 99d99
1462
wmax = -1d0
1463
mean = 0d0
1464
sigma = 0d0
1465
chi2 = 0d0
1466
non_zero = 0
1467
vol = 1d0 / dble(events * iter)
1468
knt = events
1469
1470
do i=1,maxconfigs
1471
psect(i)=0d0
1472
enddo
1473
do i=1,invar
1474
xmin(i) = xgmax
1475
xmax(i) = xgmin
1476
do j=1,ng
1477
inon_zero(j,i)=0
1478
grid(1,j,i) =0d0
1479
grid2(j,i) =0d0
1480
enddo
1481
enddo
1482
endif
1483
1484
if (iteration .eq. cur_it) then
1485
kn = kn + 1
1486
if (.true.) then !Average points to increase error estimate
1487
twgt1=twgt1+wgt !This doesn't change anything should remove
1488
iavg = iavg+1
1489
if (iavg .ge. navg) then
1490
sigma=sigma+twgt1**2
1491
iavg = 0
1492
twgt1=0d0
1493
endif
1494
else
1495
sigma = sigma + wgt**2
1496
endif
1497
if (wgt .gt. 0.) then
1498
if (wgt*iter*events .gt. wmax) then
1499
wmax=wgt*iter*events
1500
endif
1501
non_zero = non_zero + 1
1502
mean = mean + wgt
1503
if (.true. ) then
1504
c psect(ipole)=psect(ipole)+wgt*wgt/alpha(ipole) !Ohl
1505
c psect(ipole)=1d0 !Not doing multi_config
1506
else
1507
tot = 0d0
1508
do i=1,configs
1509
tot=tot+prb(i,jpnt,jplace)*alpha(i)
1510
enddo
1511
do i=1,configs
1512
if (tot .gt. 0d0) then !Pittau hep-ph/9405257
1513
psect(i)=psect(i)+wgt*wgt*prb(i,jpnt,jplace)/tot
1514
else
1515
psect(i)=psect(i)+wgt*wgt*alpha(i) !prb not set....
1516
endif
1517
enddo
1518
endif
1519
c write(123,'(2i6,1e15.5)') 1,1,wgt
1520
c write(123,'(5e15.9)') (fprb(i,jpnt,jplace),i=1,invar)
1521
c write(123,'(5e15.9)') (prb(i,jpnt,jplace),i=1,configs)
1522
do j = 1, invar
1523
i = int(xbin(point(j),j))+1
1524
if (i .gt. ng) then
1525
print*,'error i>ng',i,j,ng,point(j)
1526
i=ng
1527
endif
1528
grid(1, i, j) = grid(1, i, j) + abs(wgt)
1529
grid2(i, j) = grid2(i, j) + wgt**2
1530
c
1531
c Lines below are for multiconfiguration
1532
c
1533
c grid(1, i, j) = grid(1, i, j) +
1534
c & (abs(wgt)**2)*fprb(j,jpnt,jplace)
1535
c grid2(i, j) = grid2(i, j) + wgt**4*fprb(j,jpnt,jplace)
1536
if (abs(wgt) .gt. 0) inon_zero(i,j) = inon_zero(i,j)+1
1537
c
1538
c Here we need to look out for point(j) which has been transformed
1539
c for Briet-Wigner pole
1540
c
1541
if (j .gt. 0) then
1542
if (swidth(j) .gt. 0d0) then
1543
ddumb=0d0
1544
call untranspole(spole(j),swidth(j),
1545
& point(j),point(j),ddumb)
1546
if (point(j) .lt. 0d0) then
1547
print*,'Warning point<0',j,point(j)
1548
endif
1549
endif
1550
endif
1551
if (abs(wgt) .gt. 0) xmin(j)=min(xmin(j),point(j))
1552
if (abs(wgt) .gt. 0) xmax(j)=max(xmax(j),point(j))
1553
if (xmin(j) .lt. xgmin) then
1554
print*,'Warning xmin<0',j,xmin(j),point(j)
1555
endif
1556
xmin(j)=max(xmin(j),xgmin)
1557
end do
1558
endif
1559
c
1560
c Now if done with an iteration, print out stats, rebin, reset
1561
c
1562
c if (kn .eq. events) then
1563
if (kn .ge. max_events .and. non_zero .lt. 5) then
1564
call none_pass(max_events)
1565
endif
1566
if (non_zero .eq. events .or. (kn .gt. 200*events .and.
1567
$ non_zero .gt. 5)) then
1568
mean=mean*dble(events)/dble(non_zero)
1569
twgt1=twgt1*dble(events)/dble(non_zero)
1570
sigma=sigma+twgt1**2 !This line for averaging over points
1571
if (non_zero .eq. 0) then
1572
write(*,*) 'Error no points passed the cuts.'
1573
write(*,*) 'Try running with more points or looser cuts.'
1574
stop
1575
endif
1576
c mean = mean * iter !Used if don't have non_zero
1577
if (.true.) then
1578
mean = mean * iter *dble(non_zero)/dble(kn)
1579
knt = kn
1580
endif
1581
c
1582
c Need to fix this if averaging over navg events
1583
c
1584
c write(*,*) (sigma/vol/vol-knt*mean*mean)/dble(knt-1)/dble(knt),
1585
c & (sigma/vol/vol-knt*mean*mean*navg)/dble(knt-1)/ dble(knt)
1586
1587
if (.true.) then
1588
c vol = 1d0/(knt*iter)
1589
sigma = (sigma/vol/vol-non_zero*mean*mean*navg) !knt replaced by non_zero
1590
. / dble(knt-1) / dble(knt)
1591
else
1592
sigma = (sigma/vol/vol - knt*mean*mean)
1593
. / dble(knt-1) / dble(knt)
1594
endif
1595
1596
tmean = tmean + mean * (mean**2 / sigma)
1597
tsigma = tsigma + mean**2 / sigma
1598
chi2 = chi2 + mean**2 * (mean**2 / sigma)
1599
sigma = sqrt(abs(sigma))
1600
1601
if (cur_it .lt. 100) then
1602
ymean(cur_it) = mean
1603
ysigma(cur_it) = sigma
1604
ywmax(cur_it)= wmax*dble(non_zero)/dble(kn)
1605
yeff(cur_it)= sigma*sqrt(dble(non_zero))/mean
1606
c call sample_writehtm()
1607
endif
1608
write(*,222) 'Iteration',cur_it,'Mean: ',mean,
1609
& ' Fluctuation: ',sigma,
1610
& wmax*(dble(non_zero)/dble(kn)),
1611
& dble(non_zero)/dble(kn)*100.,'%'
1612
222 format(a10,I3,3x,a6,e10.4,a16,e10.3,e12.3,3x,f5.1,a1)
1613
1614
write(*,223) cur_it, mean, ' +- ', sigma,
1615
& sigma*sqrt(dble(non_zero))/mean
1616
223 format( i3,3x,e10.4,a,e10.4,f10.2)
1617
tot=0d0
1618
do i=1,configs
1619
tot=tot+psect(i)
1620
enddo
1621
if (configs .gt. 1)
1622
& write(*,'(8f10.5)') (psect(i)/tot, i=1,configs)
1623
c
1624
c Now set things up for generating unweighted events
1625
c
1626
if (twgt .eq. -2d0) then
1627
twgt = mean *kn/ (dble(iter)*dble(events)*dble(events))
1628
c
1629
c now scale twgt, in case have large fluctuations
1630
c
1631
1632
c twgt = twgt * max(1d0, yeff(cur_it))
1633
1634
c
1635
c For small number of events only write about 1% of events
1636
c
1637
c if (events .le. 2500) then
1638
c twgt = mean *kn*100 /
1639
c $ (dble(iter)*dble(events)*dble(events))
1640
c endif
1641
c twgt = max(twgt, maxwgt/10d0)
1642
write(*,*) 'Writing out events',twgt, yeff(cur_it)
1643
c write(*,*) mean, kn, iter, events
1644
endif
1645
c
1646
c This tells it to write out a file for unweighted events
1647
c
1648
c if(wmax*(dble(non_zero)/dble(kn)) .lt. wmax1) then
1649
if(sigma/(mean+1d-99) .lt. wmax1 .and. use_cut .ne. 0) then
1650
c wmax1 = wmax*(dble(non_zero)/dble(kn))
1651
wmax1 = sigma/(mean+1d-99)
1652
c open(26, file='ftn99',status='unknown')
1653
c write(26,fmt='(4f20.17)')
1654
c $ ((grid(2,i,j),i=1,ng),j=1,maxinvar)
1655
c write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
1656
c close(26)
1657
endif
1658
tot=0d0
1659
if (use_cut .ne. 0) then
1660
c write(*,*) 'Keeping alpha fixed'
1661
if (configs .gt. 1) then
1662
do i=1,configs
1663
alpha(i)=alpha(i)*sqrt(sqrt(psect(i))) !Pittau
1664
tot = tot+alpha(i)
1665
psect(i)=0d0
1666
enddo
1667
do i=1,configs
1668
alpha(i)=alpha(i)/tot
1669
enddo
1670
write(*,'(A)') 'Configs:'
1671
write(*,'(8f10.5)') (alpha(i),i=1,configs)
1672
endif
1673
endif
1674
c open(unit=22,file=result_file,status='old',access='append',
1675
c & err=23)
1676
c write(22,222) 'Iteration',cur_it,'Mean: ',mean,
1677
c & ' Fluctuation: ',sigma,
1678
c & wmax*(dble(non_zero)/dble(kn)),
1679
c & dble(non_zero)/dble(kn)*100.,'%'
1680
c close(22)
1681
c------
1682
c Here we will double the number of events requested for the next run
1683
c-----
1684
23 events = 2 * events
1685
vol = 1d0/dble(events*iter)
1686
knt = events
1687
twgt = mean / (dble(iter)*dble(events))
1688
c write(*,*) 'New number of events',events,twgt
1689
1690
mean = 0d0
1691
sigma = 0d0
1692
cur_it = cur_it + 1
1693
kn = 0
1694
wmax = -1d0
1695
c
1696
c do so adjusting of weights according to number of events in bin
1697
c
1698
do j=1,invar
1699
do i = 1, ng
1700
if (use_cut .ne. 2 .and.
1701
& use_cut .ne. 3 .and. use_cut .ne. 5)
1702
$ inon_zero(i,j) = 0
1703
if (use_cut .eq. 3) grid(1,i,j)=grid2(i,j)
1704
if (inon_zero(i,j) .ne. 0) then
1705
grid(1,i,j) = grid(1,i,j)
1706
& *dble(min((real(non_zero)/real(inon_zero(i,j))),
1707
$ 10000.))
1708
grid2(i,j) = grid2(i,j)
1709
& *dble(min((real(non_zero)/real(inon_zero(i,j))),
1710
$ 10000.))**2
1711
if (real(non_zero)/real(inon_zero(i,j))
1712
& .gt. 100000) then
1713
c if (j .eq. 1) then
1714
print*,'Exceeded boost',j,i,
1715
& real(non_zero)/real(inon_zero(i,j))
1716
c endif
1717
endif
1718
inon_zero(i,j) = 0
1719
endif
1720
if (use_cut .eq. 4)
1721
& reliable(i,j)=dsqrt(grid2(i,j))/grid(1,i,j)
1722
enddo
1723
enddo
1724
if (use_cut .eq. 4) then
1725
use_cut=0
1726
endif
1727
do j = 1, invar
1728
k=1
1729
c
1730
c special routines to deal with xmin cutoff
1731
c
1732
do while(grid(1,k,j) .le. 0d0 .and. k .lt. ng)
1733
k=k+1
1734
enddo
1735
1736
c if (j .eq. 1) then
1737
c open(unit=22,file='x1.dat',status='unknown')
1738
c do i=1,ng
1739
c write(22,'(i6,2e20.8)') i,grid(1,i,j),
1740
c $ dsqrt(grid2(i,j))
1741
c enddo
1742
c close(22)
1743
c endif
1744
1745
x(j)=0d0
1746
do i=1,ng
1747
x(j)=x(j)+grid(1,i,j)
1748
enddo
1749
1750
call average_grid(j,k,grid,grid2,x)
1751
1752
c if (j .eq. 1 .and. .true.) then
1753
c open(unit=22,file='x1avg.dat',status='unknown')
1754
c do i=1,ng
1755
c write(22,'(i6,2e20.8)') i,grid(1,i,1),
1756
c $ dsqrt(grid2(i,1))
1757
c enddo
1758
c close(22)
1759
c endif
1760
1761
c
1762
c Now take logs to help the rebinning converge quicker
1763
c
1764
rc = 0d0
1765
do i= k, ng
1766
xo = (1.0d-14) + grid(1, i, j) / x(j)
1767
grid(1, i, j) = ((xo - 1d0) / log(xo))**1.5 !this is 1.5
1768
rc = rc + grid(1, i, j)
1769
c write(*,*) i,rc
1770
end do
1771
rc = rc / dble(ng)
1772
k = 0
1773
xn = xgmin
1774
dr = 0d0
1775
i = 0
1776
c
1777
c Special lines to deal with xmin .ne. 0 cutoffs
1778
c
1779
c
1780
c These assume one endpoints are xgmin and xgmax
1781
c
1782
c
1783
1784
xnmin = xgmin !Endpoints for grid usually 0d0
1785
xnmax = xgmax !Endpoint for grid usually 1d0
1786
if (xmin(j)-xgmin .gt. (grid(2,2,j)-grid(2,1,j)))then
1787
xnmin = xmin(j)-(grid(2,2,j)-grid(2,1,j))/5d0
1788
i = 1
1789
dum(i)= xnmin
1790
xn = xnmin
1791
rc = rc * dble(ng)/dble(ng-i)
1792
endif
1793
dum(ng-1) = -1d0
1794
if (xgmax-xmax(j).gt.(grid(2,ng-1,j)-grid(2,ng-2,j)))then
1795
xnmax = xmax(j)+(grid(2,ng-1,j)-grid(2,ng-2,j))/5d0
1796
dum(ng-1)= xnmax
1797
rc = rc * dble(ng-i)/dble(ng-i-1)
1798
c print*,'xmax',j,xmax(j),dum(ng-1)
1799
endif
1800
1801
25 k = k + 1
1802
dr = dr + grid(1, k, j)
1803
xo = xn
1804
xn = max(grid(2, k, j),xnmin)
1805
xn = min(xn,xnmax)
1806
26 if (rc .gt. dr) goto 25
1807
1808
i = i + 1
1809
dr = dr - rc
1810
dum(i) = xn - (xn - xo) * dr / grid(1, k, j)
1811
c
1812
c Put in check for 0 width bin NEED TO FIX THIS
1813
c
1814
if (dum(ng-1) .eq. -1) then
1815
if (i .lt. ng - 1 ) goto 26
1816
else
1817
if (i .lt. ng - 2 ) goto 26
1818
endif
1819
c
1820
c Here is another fix for 0 width bins
1821
c
1822
do i=1,ng-1
1823
if (dum(i+1)-dum(i) .le. 1d-14) then
1824
c write(*,'(a,2i4,2f24.17,1e10.3)') 'Bin too small',
1825
c & j,i,dum(i),dum(i+1),dum(i+1)-dum(i)
1826
dum(i+1)=dum(i)+1d-14
1827
if (dum(i+1) .gt. xgmax) then
1828
write(*,*) 'Error in rebin',i,dum(i),dum(i+1)
1829
endif
1830
endif
1831
enddo
1832
c
1833
c Now reset counters and set new grid as necessary
1834
c
1835
do i = 1, ng - 1
1836
grid(1, i, j) = 0d0
1837
grid2(i,j) = 0d0
1838
if (use_cut .ne. 0 .and. j .gt. 0)
1839
$ grid(2, i, j) = dum(i)
1840
end do
1841
grid(1, ng, J) = 0d0
1842
grid(2, ng, J) = xgmax
1843
grid2(ng,j) = 0d0
1844
non_zero = 0
1845
1846
call sample_write_g(j,'_1')
1847
1848
end do
1849
c write(*,*) (irebin(j),j=1,dim)
1850
c open(unit=26,file='grid.dat',status='unknown')
1851
c do j=1,maxinvar
1852
c do i=1,ng
1853
c write(26,*) grid(2,i,j),j,i
1854
c enddo
1855
c enddo
1856
c close(26)
1857
c
1858
c Add test to see if we have achieved desired accuracy
1859
c
1860
if (tsigma .gt. 0d0 .and. cur_it .gt. 5 .and. accur .gt. 0d0) then
1861
1862
xmean = tmean/tsigma
1863
xchi2 = (chi2/xmean/xmean-tsigma)/dble(cur_it-2)
1864
write(*,'(a,4f8.3)') 'We got it',sqrt(xchi2/tsigma),
1865
& accur,1/sqrt(tsigma),xchi2
1866
c write(*,*) 'We got it',1d0/sqrt(tsigma), accur
1867
c if (1d0/sqrt(tsigma) .lt. accur) then
1868
if (sqrt(xchi2/tsigma) .lt. accur) then
1869
write(*,*) 'Finished due to accuracy',sqrt(xchi2/tsigma), accur
1870
tmean = tmean / tsigma
1871
if (cur_it .gt. 2) then
1872
chi2 = (chi2/tmean/tmean-tsigma)/dble(cur_it-2)
1873
else
1874
chi2=0d0
1875
endif
1876
tsigma = tmean / sqrt(tsigma)
1877
write(*, 80) real(tmean), real(tsigma), real(chi2)
1878
if (use_cut .ne. 0) then
1879
open(26, file='ftn26',status='unknown')
1880
write(26,fmt='(4f20.17)')
1881
$ ((grid(2,i,j),i=1,ng),j=1,maxinvar)
1882
write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
1883
close(26)
1884
endif
1885
call sample_writehtm()
1886
c open(unit=22,file=result_file,status='old',
1887
c $ access='append',err=122)
1888
c write(22, 80) real(tmean), real(tsigma), real(chi2)
1889
c 122 close(22)
1890
tsigma = tsigma*sqrt(chi2) !This gives the 68% confidence cross section
1891
call store_events
1892
cur_it = iter+2
1893
return
1894
endif
1895
endif
1896
c
1897
c New check to see if we need to keep integrating this one or not.
1898
c
1899
if (cur_it .gt. 3 .and. accur .lt. 0d0) then !Check # unweighted
1900
c
1901
c Lets get the actual number instead
1902
c tjs 5/22/2007
1903
c
1904
c nun = n_unwgted()
1905
c write(*,*) 'Estimated events',nun, accur
1906
call store_events
1907
nun = neventswritten
1908
write(*,*) "Checking number of events",accur,nun
1909
if (nun .gt. -accur)then
1910
tmean = tmean / tsigma
1911
if (cur_it .gt. 2) then
1912
chi2 = (chi2/tmean/tmean-tsigma)/dble(cur_it-2)
1913
else
1914
chi2=0d0
1915
endif
1916
tsigma = tmean / sqrt(tsigma)
1917
write(*, 80) real(tmean), real(tsigma), real(chi2)
1918
if (use_cut .ne. 0) then
1919
open(26, file='ftn26',status='unknown')
1920
write(26,fmt='(4f20.17)')
1921
$ ((grid(2,i,j),i=1,ng),j=1,maxinvar)
1922
write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
1923
close(26)
1924
endif
1925
call sample_writehtm()
1926
1927
c open(unit=22,file=result_file,status='old',
1928
c $ access='append',err=129)
1929
c write(22, 80) real(tmean), real(tsigma), real(chi2)
1930
c 129 close(22)
1931
tsigma = tsigma*sqrt(chi2) !This gives the 68% confidence cross section
1932
cur_it = iter+20
1933
return
1934
endif
1935
endif
1936
1937
1938
if (cur_it .gt. iter) then
1939
call store_events
1940
tmean = tmean / tsigma
1941
chi2 = (chi2 / tmean / tmean - tsigma) / dble(iter - 1)
1942
tsigma = tmean / sqrt(tsigma)
1943
write(*, 80) real(tmean), real(tsigma), real(chi2)
1944
80 format(/1X,79(1H-)/1X,23HAccumulated results: ,
1945
. 10HIntegral =,e12.4/24X,10HStd dev =,e12.4/
1946
. 13X,21HChi**2 per DoF. =,f12.4/1X,79(1H-))
1947
if (use_cut .ne. 0) then
1948
open(26, file='ftn26',status='unknown')
1949
write(26,fmt='(4f20.17)')
1950
$ ((grid(2,i,j),i=1,ng),j=1,maxinvar)
1951
write(26,fmt='(4f20.17)') (alpha(i),i=1,maxconfigs)
1952
close(26)
1953
endif
1954
call sample_writehtm()
1955
c open(unit=22,file=result_file,status='old',
1956
c $ access='append',err=123)
1957
c write(22, 80) real(tmean), real(tsigma), real(chi2)
1958
c 123 close(22)
1959
tsigma = tsigma*sqrt(chi2) !This gives the 68% confidence cross section
1960
else
1961
c
1962
c Starting new iteration, should clean out stored events
1963
c and start fresh
1964
c
1965
c nun = n_unwgted()
1966
c write(*,*) 'Estimated unweighted events ', nun
1967
call clear_Events
1968
endif
1969
endif
1970
else
1971
endif
1972
end
1973
1974
subroutine none_pass(max_events)
1975
c*************************************************************************
1976
c Special break to handle case where no events are passing cuts
1977
c We'll set the cross section to zero here.
1978
c*************************************************************************
1979
implicit none
1980
c
1981
c Constants
1982
c
1983
include 'genps.inc'
1984
c
1985
c Arguments
1986
c
1987
integer max_events
1988
c
1989
c Global
1990
c
1991
integer nsteps
1992
character*40 result_file,where_file
1993
common /sample_status/result_file,where_file,nsteps
1994
1995
c----
1996
c Begin Code
1997
c----
1998
write(*,*) 'No points passed cuts!'
1999
write(*,*) 'Loosen cuts or increase max_events',max_events
2000
c open(unit=22,file=result_file,status='old',access='append',
2001
c & err=23)
2002
c write(22,222) 'Iteration',0,'Mean: ',0d0,
2003
c & ' Fluctuation: ',0d0,
2004
c & 0d0,
2005
c & 0d0,'%'
2006
c 23 close(22)
2007
222 format(a10,I3,3x,a6,e10.4,a16,e10.3,e12.3,3x,f5.1,a1)
2008
2009
open(unit=66,file='results.dat',status='unknown')
2010
write(66,*) 0,0,0,0,0,0,0,1
2011
write(66,'(i4,4e15.5)') 1,0d0,0d0,0d0,0d0
2012
close(66)
2013
2014
c Remove file events.lhe (otherwise event combination gets screwed up)
2015
write(*,*) 'Deleting file events.lhe'
2016
open(unit=67,file='events.lhe',status='unknown')
2017
write(67,*)
2018
close(67)
2019
2020
stop
2021
end
2022
2023
subroutine average_grid(j,k,grid,grid2,x)
2024
c**************************************************************************
2025
c Special routine to deal with averaging over the grid bins
2026
c This routine starts averaging at bin k rather than bin 1 so that
2027
c one can accommodate cutoffs. With k=1 this should give the
2028
c standard sample/vegas/bases averaging results.
2029
c
2030
c Also stops averaging when reaches maximum value
2031
c
2032
c**************************************************************************
2033
implicit none
2034
c
2035
c Constants
2036
c
2037
include 'genps.inc'
2038
c
2039
c Arguments
2040
c
2041
integer j,k
2042
double precision grid(2,ng,0:maxinvar),grid2(0:ng,maxinvar)
2043
double precision x(maxinvar)
2044
c
2045
c Local
2046
c
2047
integer i,kmax
2048
double precision xo,xn
2049
c-----
2050
c Begin Code
2051
c-----
2052
kmax=k
2053
do i=k+1,ng-1
2054
if (grid(1,i,j) .gt. 0d0) kmax=i
2055
enddo
2056
xo = grid(1,k,j)
2057
xn = grid(1,k+1,j)
2058
grid(1,k,j) = (xo+xn)/2d0
2059
x(j) = grid(1,k,j)
2060
c do i=k+1,ng-1 !Original without kmax stuff
2061
do i=k+1,kmax-1
2062
grid(1, i, j) = xo + xn
2063
xo = xn
2064
xn = grid(1, i+1, j)
2065
grid(1, i, j) = (grid(1, i, j) + xn) / 3d0
2066
x(j) = x(j) + grid(1, i, j)
2067
end do
2068
c grid(1, ng, j) = (xn + xo) / 2d0 !Original without kmax stuff
2069
grid(1, kmax, j) = (xn + xo) / 2d0
2070
x(j) = x(j) + grid(1, ng, j)
2071
end
2072
2073
double precision function xbin(y,j)
2074
c**************************************************************************
2075
c Subroutine to determine which value y will map to give you the
2076
c value of x when put through grid j. That is what random number
2077
c do you need to be given to get the value x out of grid j and will be
2078
c between 0 < x < ng.
2079
c**************************************************************************
2080
implicit none
2081
c
2082
c Constants
2083
c
2084
include 'genps.inc'
2085
double precision tol
2086
parameter (tol=1d-12)
2087
c
2088
c Arguments
2089
c
2090
double precision y
2091
integer j
2092
c
2093
c Local
2094
c
2095
integer i,jl,ju
2096
double precision x,xo
2097
c
2098
c Global
2099
c
2100
double precision grid(2, ng, 0:maxinvar)
2101
common /data_grid/ grid
2102
double precision spole(maxinvar),swidth(maxinvar),bwjac
2103
common/to_brietwigner/spole ,swidth ,bwjac
2104
c
2105
c Data
2106
c
2107
data spole,swidth/maxinvar*0d0,maxinvar*0d0/
2108
c-----
2109
c Begin Code
2110
c-----
2111
bwjac = 1d0
2112
if (j .gt. 0) then
2113
if (swidth(j) .gt. 0d0) then
2114
call untranspole(spole(j),swidth(j),x,y,bwjac)
2115
else
2116
x=y
2117
endif
2118
else
2119
x=y
2120
endif
2121
if (x .eq. xgmax) then
2122
i=ng
2123
xbin = dble(ng)
2124
elseif (x .eq. xgmin) then
2125
xbin=0d0
2126
elseif(x .le. grid(2,1,j)) then
2127
i=1
2128
xo = grid(2,i,j)-xgmin
2129
xbin = dble(i)+(x-grid(2,i,j))/xo
2130
else
2131
jl = 1
2132
ju = ng
2133
do while (ju-jl .gt. 1) !Binary search
2134
i = (ju-jl)/2+jl
2135
if (grid(2,i,j) .le. x) then
2136
jl=i
2137
else
2138
ju=i
2139
endif
2140
enddo
2141
i=ju
2142
xo = grid(2,i,j)-grid(2,i-1,j)
2143
xbin = dble(i)+(x-grid(2,i,j))/xo
2144
endif
2145
c jbin=i
2146
c x =
2147
c if (x+tol .gt. grid(2,i,j) .and. i .ne. ng) then
2148
c write(*,'(a,2e23.16,e9.2)') 'Warning in DSAMPLE:JBIN ',
2149
c & x,grid(2,i,j),tol
2150
c x=2d0*grid(2,i,j)-x
2151
c jbin=i+1
2152
c endif
2153
end
2154
2155
2156
subroutine sample_write_g(idim,cpost)
2157
c**************************************************************************
2158
c Writes out grid in function form for dimension i with extension cpost
2159
c
2160
c**************************************************************************
2161
implicit none
2162
c
2163
c Constants
2164
c
2165
include 'genps.inc'
2166
c
2167
c Arguments
2168
c
2169
integer idim
2170
character*(*) cpost
2171
c
2172
c Local
2173
c
2174
character*60 fname
2175
integer i
2176
double precision xo,yo
2177
c
2178
c Global
2179
c
2180
double precision grid(2, ng, 0:maxinvar)
2181
common /data_grid/ grid
2182
2183
c-----
2184
c Begin Code
2185
c-----
2186
return
2187
if (idim .lt. 1 .or. idim .gt.maxinvar) then
2188
write(*,*) 'Error invalid dimension in sample_write_f',idim
2189
return
2190
endif
2191
if (idim .lt. 10) then
2192
write(fname,'(a,i1,a,a)') 'g_',idim,cpost,'.dat'
2193
elseif (idim .lt. 100) then
2194
write(fname,'(a,i2,a,a)') 'g_',idim,cpost,'.dat'
2195
endif
2196
open(unit=21,file=fname,status='unknown',err=99)
2197
do i=1,ng-1
2198
xo = (grid(2,i,idim)+grid(2,i+1,idim))/2d0
2199
yo =1d0/(-grid(2,i,idim)+grid(2,i+1,idim))
2200
write(21,*) xo,yo
2201
enddo
2202
close(21)
2203
return
2204
99 write(*,*) 'Error opening file ',fname
2205
end
2206
2207
function ran1(idum)
2208
dimension r(97)
2209
parameter (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6)
2210
parameter (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6)
2211
parameter (m3=243000,ia3=4561,ic3=51349)
2212
data iff /0/
2213
save r, ix1, ix2, ix3
2214
if (idum.lt.0.or.iff.eq.0) then
2215
iff=1
2216
ix1=mod(ic1-idum,m1)
2217
ix1=mod(ia1*ix1+ic1,m1)
2218
ix2=mod(ix1,m2)
2219
ix1=mod(ia1*ix1+ic1,m1)
2220
ix3=mod(ix1,m3)
2221
do 11 j=1,97
2222
ix1=mod(ia1*ix1+ic1,m1)
2223
ix2=mod(ia2*ix2+ic2,m2)
2224
r(j)=(float(ix1)+float(ix2)*rm2)*rm1
2225
11 continue
2226
idum=1
2227
endif
2228
ix1=mod(ia1*ix1+ic1,m1)
2229
ix2=mod(ia2*ix2+ic2,m2)
2230
ix3=mod(ia3*ix3+ic3,m3)
2231
j=1+(97*ix3)/m3
2232
if(j.gt.97.or.j.lt.1) stop
2233
ran1=r(j)
2234
r(j)=(float(ix1)+float(ix2)*rm2)*rm1
2235
return
2236
end
2237
2238
2239
Loggerhead is a web-based interface for Breezy
Version: 2.0.1