Merge lp:~maddevelopers/mg5amcnlo/virtual_polyfit_EW into lp:~maddevelopers/mg5amcnlo/3.0.2

Proposed by Rikkert Frederix on 2019-09-12
Status: Merged
Merged at revision: 973
Proposed branch: lp:~maddevelopers/mg5amcnlo/virtual_polyfit_EW
Merge into: lp:~maddevelopers/mg5amcnlo/3.0.2
Diff against target: 1828 lines (+1113/-345)
15 files modified
Template/NLO/Cards/FKS_params.dat (+10/-6)
Template/NLO/SubProcesses/BinothLHA.f (+1/-2)
Template/NLO/SubProcesses/BinothLHA_OLP.f (+1/-1)
Template/NLO/SubProcesses/FKSParamReader.f (+0/-252)
Template/NLO/SubProcesses/FKSParams.f90 (+248/-0)
Template/NLO/SubProcesses/FKSParams.inc (+0/-31)
Template/NLO/SubProcesses/check_poles.f (+1/-1)
Template/NLO/SubProcesses/driver_mintFO.f (+1/-3)
Template/NLO/SubProcesses/driver_mintMC.f (+1/-3)
Template/NLO/SubProcesses/fks_singular.f (+31/-12)
Template/NLO/SubProcesses/makefile_fks_dir (+3/-3)
Template/NLO/SubProcesses/mint_module.f90 (+58/-20)
Template/NLO/SubProcesses/polfit.f (+720/-0)
madgraph/interface/amcatnlo_run_interface.py (+27/-2)
madgraph/iolibs/export_fks.py (+11/-9)
To merge this branch: bzr merge lp:~maddevelopers/mg5amcnlo/virtual_polyfit_EW
Reviewer Review Type Date Requested Status
marco zaro 2019-09-12 Approve on 2019-09-17
Review via email: mp+372681@code.launchpad.net

Description of the change

Replaced the grids to determine the approximate virtual by a polynomial fit using the polfit.f code (by Shampine, Davenport and Huddleston). It's still based on doing this one dimension at a time, but the polynomial fit is smoother than the grids, since the latter required rather coarse binning and had fixed bins which could not be adjusted to follow the importance sampling.

Also needed to replace the FKSParamReader with a fortran90 module to be able to use it correctly in the mint_module.f90

To post a comment you must log in.
marco zaro (marco-zaro) wrote :

Thanks Rikkert,
the code changes look correct, for what I can understand of the polyfit module
I have tried p p > t t~ [QCD QED], results are compatible within the uncertainty (also split-order by split-order). The Polyfit run is generally (much) faster:
If I require 1 per mil accuracy for that process, without polyfit it took 18 minutes on my macbook, while with polyfit less than 8, so great job!!

review: Approve

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'Template/NLO/Cards/FKS_params.dat'
2--- Template/NLO/Cards/FKS_params.dat 2017-08-04 19:52:15 +0000
3+++ Template/NLO/Cards/FKS_params.dat 2019-09-12 11:14:33 +0000
4@@ -153,6 +153,16 @@
5 0.005d0
6 ! Default :: 0.005d0
7 !
8+! Use a polynomial fit for the virtual corrections. Gets a better accuracy,
9+! but uses more memmory and disk space.
10+#UsePolyVirtual
11+.True.
12+! Default :: .True.
13+!
14+! ==========================================================================
15+! Combination of matrix elements
16+! ==========================================================================
17+!
18 ! For fixed order calculations only (ignored for (N)LO+PS runs): This
19 ! parameter determines if parton flavour configurations with idential
20 ! matrix elements should be separated at the level of the analysis or
21@@ -170,11 +180,5 @@
22 ! Default :: .false.
23 !
24 ! ==========================================================================
25-! Arbitrary numerical parameters used in the FKS formalism
26-! ==========================================================================
27-!
28-! To be implemented by the FKS authors
29-!
30-! ==========================================================================
31 ! End of FKS_params.dat file
32 ! ==========================================================================
33
34=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
35--- Template/NLO/SubProcesses/BinothLHA.f 2017-06-23 14:14:03 +0000
36+++ Template/NLO/SubProcesses/BinothLHA.f 2019-09-12 11:14:33 +0000
37@@ -4,12 +4,11 @@
38 c that calls the OLP and returns the virtual weights. For convenience
39 c also the born_wgt is passed to this subroutine.
40 c
41+ use FKSParams
42 implicit none
43 include "nexternal.inc"
44 include "coupl.inc"
45 include 'born_nhel.inc'
46-c general MadFKS parameters
47- include 'FKSParams.inc'
48 double precision pi, zero,mone
49 parameter (pi=3.1415926535897932385d0)
50 parameter (zero=0d0)
51
52=== modified file 'Template/NLO/SubProcesses/BinothLHA_OLP.f'
53--- Template/NLO/SubProcesses/BinothLHA_OLP.f 2015-08-13 09:26:37 +0000
54+++ Template/NLO/SubProcesses/BinothLHA_OLP.f 2019-09-12 11:14:33 +0000
55@@ -6,11 +6,11 @@
56 c
57 C************************************************************************
58 c
59+ use FKSParams
60 implicit none
61 include "nexternal.inc"
62 include "coupl.inc"
63 include "Binoth_proc.inc"
64- include "FKSParams.inc"
65 double precision pi
66 parameter (pi=3.1415926535897932385d0)
67 double precision pin(0:3,nexternal-1),p(0:4,nexternal-1)
68
69=== removed file 'Template/NLO/SubProcesses/FKSParamReader.f'
70--- Template/NLO/SubProcesses/FKSParamReader.f 2017-09-01 10:34:24 +0000
71+++ Template/NLO/SubProcesses/FKSParamReader.f 1970-01-01 00:00:00 +0000
72@@ -1,252 +0,0 @@
73- subroutine FKSParamReader(filename, printParam, force)
74-
75- implicit none
76-
77- logical HasReadOnce, force
78- data HasReadOnce/.FALSE./
79-
80-
81- character(*) filename
82- CHARACTER*64 buff, buff2, mode
83- include "FKSParams.inc"
84- include "orders.inc"
85-
86- integer i,j
87-
88- logical printParam, couldRead, paramPrinted
89- data paramPrinted/.FALSE./
90- couldRead=.False.
91-! Default parameters
92-
93- if (HasReadOnce.and..not.force) then
94- goto 901
95- endif
96-
97-! Make sure to have default parameters if not set
98-! in the FKSParams.dat card (if it is an old one for instance)
99- call DefaultFKSParam()
100-
101- open(68, file=fileName, err=676, action='READ')
102- do
103- read(68,*,end=999) buff
104- if(index(buff,'#').eq.1) then
105-
106- if (buff .eq. '#IRPoleCheckThreshold') then
107- read(68,*,end=999) IRPoleCheckThreshold
108- if (IRPoleCheckThreshold .lt. -1.01d0 ) then
109- stop 'IRPoleCheckThreshold must be >= -1.0d0.'
110- endif
111-
112- elseif (buff .eq. '#PrecisionVirtualAtRunTime') then
113- read(68,*,end=999) PrecisionVirtualAtRunTime
114- if (IRPoleCheckThreshold .lt. -1.01d0 ) then
115- stop 'PrecisionVirtualAtRunTime must be >= -1.0d0.'
116- endif
117-
118- else if (buff .eq. '#NHelForMCoverHels') then
119- read(68,*,end=999) NHelForMCoverHels
120- if (NHelForMCoverHels .lt. -1) then
121- stop 'NHelForMCoverHels must be >= -1.'
122- endif
123- else if (buff .eq. '#QCD^2==') then
124- read(68,*,end=999) QCD_squared_selected
125- if (QCD_squared_selected .lt. -1) then
126- stop 'QCD_squared_selected must be >= -1.'
127- endif
128- else if (buff .eq. '#QED^2==') then
129- read(68,*,end=999) QED_squared_selected
130- if (QED_squared_selected .lt. -1) then
131- stop 'QED_squared_selected must be >= -1.'
132- endif
133- else if (buff .eq. '#VirtualFraction') then
134- read(68,*,end=999) Virt_fraction
135- if (Virt_fraction .lt. 0 .or. virt_fraction .gt.1) then
136- stop 'VirtualFraction should be a fraction'/
137- $ /' between 0 and 1'
138- endif
139- else if (buff .eq. '#MinVirtualFraction') then
140- read(68,*,end=999) Min_Virt_fraction
141- if (min_virt_fraction .lt. 0 .or. min_virt_fraction .gt.1)
142- $ then
143- stop 'VirtualFraction should be a fraction'/
144- $ /' between 0 and 1'
145- endif
146- else if (buff .eq. '#SeparateFlavourConfigurations') then
147- read(68,*,end=999) separate_flavour_configs
148-
149- else if (buff .eq. '#VetoedContributionTypes') then
150- read(68,*,end=999) VetoedContributionTypes(0)
151- if (VetoedContributionTypes(0) .lt. 0.or.
152- & VetoedContributionTypes(0) .gt. maxContribsSelected) then
153- write(*,*) 'VetoedContributionTypes length should be '/
154- & /'>= 0 and <=',maxContribsSelected
155- stop 'Format error in FKS_params.dat.'
156- endif
157- read(68,*,end=999) (VetoedContributionTypes(I),I=1,
158- $ VetoedContributionTypes(0))
159- do I=1,VetoedContributionTypes(0)
160- if (VetoedContributionTypes(I).lt.1.or.
161- $ VetoedContributionTypes(I).gt.maxContribType) then
162- write(*,*) 'VetoedContributionTypes must be >=1 and '/
163- & /'<=',maxContribType
164- stop 'Format error in FKS_params.dat.'
165- endif
166- enddo
167- do I=VetoedContributionTypes(0)+1,maxContribsSelected
168- VetoedContributionTypes(I)=-1
169- enddo
170-
171- else if (buff .eq. '#SelectedContributionTypes') then
172- read(68,*,end=999) SelectedContributionTypes(0)
173- if (SelectedContributionTypes(0) .lt. 0 .or.
174- & SelectedContributionTypes(0) .gt. maxContribsSelected) then
175- write(*,*) 'SelectedContributionTypes length should be '/
176- & /'>= 0 and <=',maxContribsSelected
177- stop 'Format error in FKS_params.dat.'
178- endif
179- read(68,*,end=999) (SelectedContributionTypes(I),I=1,
180- $ SelectedContributionTypes(0))
181- do I=1,SelectedContributionTypes(0)
182- if (SelectedContributionTypes(I).lt.1.or.
183- $ SelectedContributionTypes(I).gt.maxContribType) then
184- write(*,*) 'SelectedContributionTypes must be >=1 and '/
185- & /'<=',maxContribType
186- stop 'Format error in FKS_params.dat.'
187- endif
188- enddo
189- do I=SelectedContributionTypes(0)+1,maxContribsSelected
190- SelectedContributionTypes(I)=-1
191- enddo
192-
193- else if (buff .eq. '#SelectedCouplingOrders') then
194- read(68,*,end=999) SelectedCouplingOrders(1,0)
195- if (SelectedCouplingOrders(1,0) .lt. 0 .or.
196- & SelectedCouplingOrders(1,0) .gt. maxCouplingsSelected) then
197- write(*,*) 'SelectedCouplingOrders length should be >='/
198- & /' 0 and <=',maxCouplingsSelected
199- stop 'Format error in FKS_params.dat.'
200- endif
201- do j = 2, maxCouplingTypes
202- SelectedCouplingOrders(j,0) = SelectedCouplingOrders(1,0)
203- enddo
204- do j = 1, SelectedCouplingOrders(1,0)
205- read(68,*,end=999) (SelectedCouplingOrders(i,j),i=1,
206- $ nsplitorders)
207- do i=nsplitorders+1,maxCouplingTypes
208- SelectedCouplingOrders(i,j)=-1
209- enddo
210- enddo
211-
212- else
213- write(*,*) 'The parameter name ',buff(2:),
214- &' is not reckognized.'
215- stop 'Format error in FKS_params.dat.'
216- endif
217- endif
218- enddo
219- 999 continue
220- couldRead=.True.
221- goto 998
222-
223- 676 continue
224- write(*,*) 'ERROR :: MadFKS parameter file ',fileName,
225- &' could not be found or is malformed. Please specify it.'
226- stop
227-C Below is the code if one desires to let the code continue with
228-C a non existing or malformed parameter file
229- write(*,*) 'WARNING :: The file ',fileName,' could not be ',
230- & ' open or did not contain the necessary information. The ',
231- & ' default MadFKS parameters will be used.'
232- call DefaultFKSParam()
233- goto 998
234-
235- 998 continue
236-
237- if(printParam.and..not.paramPrinted) then
238- write(*,*)
239- & '==============================================================='
240- if (couldRead) then
241- write(*,*) 'INFO: MadFKS read these parameters from '
242- &,filename
243- else
244- write(*,*) 'INFO: MadFKS used the default parameters.'
245- endif
246- write(*,*)
247- & '==============================================================='
248- write(*,*) ' > IRPoleCheckThreshold = ',IRPoleCheckThreshold
249- write(*,*) ' > PrecisionVirtualAtRunTime = '
250- $ ,PrecisionVirtualAtRunTime
251- if (SelectedContributionTypes(0).gt.0) then
252- write(*,*) ' > SelectedContributionTypes = ',
253- & (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
254- else
255- write(*,*) ' > SelectedContributionTypes = All'
256- endif
257- if (VetoedContributionTypes(0).gt.0) then
258- write(*,*) ' > VetoedContributionTypes = ',
259- & (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
260- else
261- write(*,*) ' > VetoedContributionTypes = None'
262- endif
263- if (QCD_squared_selected.eq.-1) then
264- write(*,*) ' > QCD_squared_selected = All'
265- else
266- write(*,*) ' > QCD_squared_selected = ',QCD_squared_selected
267- endif
268- if (QED_squared_selected.eq.-1) then
269- write(*,*) ' > QED_squared_selected = All'
270- else
271- write(*,*) ' > QED_squared_selected = ',QED_squared_selected
272- endif
273- if (SelectedCouplingOrders(1,0).gt.0) then
274- do j=1,SelectedCouplingOrders(1,0)
275- write(*,*) ' > SelectedCouplingOrders(',j,') = ',
276- & (SelectedCouplingOrders(i,j),i=1,nsplitorders)
277- enddo
278- else
279- write(*,*) ' > SelectedCouplingOrders = All'
280- endif
281- write(*,*) ' > NHelForMCoverHels = ',NHelForMCoverHels
282- write(*,*) ' > VirtualFraction = ',Virt_fraction
283- write(*,*) ' > MinVirtualFraction = ',Min_virt_fraction
284- write(*,*) ' > SeparateFlavourConfigs = ',separate_flavour_configs
285- write(*,*)
286- & '==============================================================='
287- paramPrinted=.TRUE.
288- endif
289-
290- close(68)
291- HasReadOnce=.TRUE.
292- 901 continue
293- end
294-
295- subroutine DefaultFKSParam()
296-
297- implicit none
298- integer i,j
299- include "FKSParams.inc"
300-
301- IRPoleCheckThreshold=1.0d-5
302- NHelForMCoverHels=5
303- PrecisionVirtualAtRunTime=1d-3
304- Virt_fraction=1d0
305- QED_squared_selected=-1
306- QCD_squared_selected=-1
307- Min_virt_fraction=0.005d0
308- separate_flavour_configs=.false.
309- IncludeBornContributions=.true.
310- SelectedContributionTypes(0)=0
311- VetoedContributionTypes(0)=0
312- do i=1, maxContribsSelected
313- SelectedContributionTypes(I)=-1
314- VetoedContributionTypes(I)=-1
315- enddo
316- do j=1,maxCouplingTypes
317- SelectedCouplingOrders(j,0) = 0
318- enddo
319- do j=1,maxCouplingsSelected
320- do i=1,maxCouplingTypes
321- SelectedCouplingOrders(i,j) = -1
322- enddo
323- enddo
324- end
325
326=== added file 'Template/NLO/SubProcesses/FKSParams.f90'
327--- Template/NLO/SubProcesses/FKSParams.f90 1970-01-01 00:00:00 +0000
328+++ Template/NLO/SubProcesses/FKSParams.f90 2019-09-12 11:14:33 +0000
329@@ -0,0 +1,248 @@
330+!====================================================================
331+!
332+! Define common block with all general parameters used by MadFKS
333+! See their definitions in the file FKS_params.dat.
334+!
335+!====================================================================
336+module FKSParams
337+ character(len=64), parameter :: paramFileName='FKS_params.dat'
338+ integer,parameter :: maxContribsSelected=100, &
339+ maxCouplingsSelected=100, &
340+ maxContribType=15, &
341+ maxCouplingTypes=20
342+ real*8 :: IRPoleCheckThreshold,Virt_fraction, PrecisionVirtualAtRunTime,Min_virt_fraction
343+ integer :: NHelForMCoverHels,VetoedContributionTypes(0:maxContribsSelected), &
344+ SelectedContributionTypes(0:maxContribsSelected),QED_squared_selected, &
345+ SelectedCouplingOrders(maxCouplingTypes,0:maxCouplingsSelected), &
346+ QCD_squared_selected
347+ logical :: separate_flavour_configs,IncludeBornContributions,use_poly_virtual
348+
349+contains
350+
351+ subroutine FKSParamReader(filename, printParam, force)
352+ ! Reads the file 'filename' and sets the parameters found in that file.
353+ implicit none
354+ logical, save :: HasReadOnce=.False.,paramPrinted=.false.
355+ logical :: force,couldRead,printParam
356+ character(*) :: filename
357+ CHARACTER(len=64) :: buff, buff2, mode
358+ include "orders.inc"
359+ integer :: i,j
360+ couldRead=.False.
361+ if (HasReadOnce.and..not.force) then
362+ goto 901
363+ endif
364+! Make sure to have default parameters if not set in the FKSParams.dat card
365+! (if it is an old one for instance)
366+ call DefaultFKSParam()
367+! Overwrite the default parameters from file:
368+ open(68, file=fileName, err=676, action='READ')
369+ do
370+ read(68,*,end=999) buff
371+ if(index(buff,'#').eq.1) then
372+ if (buff .eq. '#IRPoleCheckThreshold') then
373+ read(68,*,end=999) IRPoleCheckThreshold
374+ if (IRPoleCheckThreshold .lt. -1.01d0 ) then
375+ stop 'IRPoleCheckThreshold must be >= -1.0d0.'
376+ endif
377+ elseif (buff .eq. '#PrecisionVirtualAtRunTime') then
378+ read(68,*,end=999) PrecisionVirtualAtRunTime
379+ if (IRPoleCheckThreshold .lt. -1.01d0 ) then
380+ stop 'PrecisionVirtualAtRunTime must be >= -1.0d0.'
381+ endif
382+ else if (buff .eq. '#NHelForMCoverHels') then
383+ read(68,*,end=999) NHelForMCoverHels
384+ if (NHelForMCoverHels .lt. -1) then
385+ stop 'NHelForMCoverHels must be >= -1.'
386+ endif
387+ else if (buff .eq. '#QCD^2==') then
388+ read(68,*,end=999) QCD_squared_selected
389+ if (QCD_squared_selected .lt. -1) then
390+ stop 'QCD_squared_selected must be >= -1.'
391+ endif
392+ else if (buff .eq. '#QED^2==') then
393+ read(68,*,end=999) QED_squared_selected
394+ if (QED_squared_selected .lt. -1) then
395+ stop 'QED_squared_selected must be >= -1.'
396+ endif
397+ else if (buff .eq. '#VirtualFraction') then
398+ read(68,*,end=999) Virt_fraction
399+ if (Virt_fraction .lt. 0 .or. virt_fraction .gt.1) then
400+ stop 'VirtualFraction should be a fraction between 0 and 1'
401+ endif
402+ else if (buff .eq. '#MinVirtualFraction') then
403+ read(68,*,end=999) Min_Virt_fraction
404+ if (min_virt_fraction .lt. 0 .or. min_virt_fraction .gt.1) then
405+ stop 'VirtualFraction should be a fraction between 0 and 1'
406+ endif
407+ else if (buff .eq. '#SeparateFlavourConfigurations') then
408+ read(68,*,end=999) separate_flavour_configs
409+ else if (buff .eq. '#UsePolyVirtual') then
410+ read(68,*,end=999) use_poly_virtual
411+ else if (buff .eq. '#VetoedContributionTypes') then
412+ read(68,*,end=999) VetoedContributionTypes(0)
413+ if (VetoedContributionTypes(0) .lt. 0.or. &
414+ VetoedContributionTypes(0) .gt. maxContribsSelected) then
415+ write(*,*) 'VetoedContributionTypes length should be >= 0 and <=', &
416+ maxContribsSelected
417+ stop 'Format error in FKS_params.dat.'
418+ endif
419+ read(68,*,end=999) (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
420+ do I=1,VetoedContributionTypes(0)
421+ if (VetoedContributionTypes(I).lt.1.or. &
422+ VetoedContributionTypes(I).gt.maxContribType) then
423+ write(*,*) 'VetoedContributionTypes must be >=1 and <=',maxContribType
424+ stop 'Format error in FKS_params.dat.'
425+ endif
426+ enddo
427+ do I=VetoedContributionTypes(0)+1,maxContribsSelected
428+ VetoedContributionTypes(I)=-1
429+ enddo
430+
431+ else if (buff .eq. '#SelectedContributionTypes') then
432+ read(68,*,end=999) SelectedContributionTypes(0)
433+ if (SelectedContributionTypes(0) .lt. 0 .or. &
434+ SelectedContributionTypes(0) .gt. maxContribsSelected) then
435+ write(*,*) 'SelectedContributionTypes length should be >= 0 and <=', &
436+ maxContribsSelected
437+ stop 'Format error in FKS_params.dat.'
438+ endif
439+ read(68,*,end=999) (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
440+ do I=1,SelectedContributionTypes(0)
441+ if (SelectedContributionTypes(I).lt.1.or. &
442+ SelectedContributionTypes(I).gt.maxContribType) then
443+ write(*,*) 'SelectedContributionTypes must be >=1 and <=',maxContribType
444+ stop 'Format error in FKS_params.dat.'
445+ endif
446+ enddo
447+ do I=SelectedContributionTypes(0)+1,maxContribsSelected
448+ SelectedContributionTypes(I)=-1
449+ enddo
450+ else if (buff .eq. '#SelectedCouplingOrders') then
451+ read(68,*,end=999) SelectedCouplingOrders(1,0)
452+ if (SelectedCouplingOrders(1,0) .lt. 0 .or. &
453+ SelectedCouplingOrders(1,0) .gt. maxCouplingsSelected) then
454+ write(*,*) 'SelectedCouplingOrders length should be >= 0 and <=', &
455+ maxCouplingsSelected
456+ stop 'Format error in FKS_params.dat.'
457+ endif
458+ do j = 2, maxCouplingTypes
459+ SelectedCouplingOrders(j,0) = SelectedCouplingOrders(1,0)
460+ enddo
461+ do j = 1, SelectedCouplingOrders(1,0)
462+ read(68,*,end=999) (SelectedCouplingOrders(i,j),i=1,nsplitorders)
463+ do i=nsplitorders+1,maxCouplingTypes
464+ SelectedCouplingOrders(i,j)=-1
465+ enddo
466+ enddo
467+ else
468+ write(*,*) 'The parameter name ',buff(2:),'is not reckognized.'
469+ stop 'Format error in FKS_params.dat.'
470+ endif
471+ endif
472+ enddo
473+999 continue
474+ couldRead=.True.
475+ goto 998
476+
477+676 continue
478+ write(*,*) 'ERROR :: MadFKS parameter file ',fileName, &
479+ ' could not be found or is malformed. Please specify it.'
480+ stop
481+ ! Below is the code if one desires to let the code continue with
482+ ! a non existing or malformed parameter file
483+ write(*,*) 'WARNING :: The file ',fileName,' could not be ', &
484+ ' open or did not contain the necessary information. The ', &
485+ ' default MadFKS parameters will be used.'
486+ call DefaultFKSParam()
487+ goto 998
488+998 continue
489+
490+ if(printParam.and..not.paramPrinted) then
491+ write(*,*) &
492+ '==============================================================='
493+ if (couldRead) then
494+ write(*,*) 'INFO: MadFKS read these parameters from ',filename
495+ else
496+ write(*,*) 'INFO: MadFKS used the default parameters.'
497+ endif
498+ write(*,*) &
499+ '==============================================================='
500+ write(*,*) ' > IRPoleCheckThreshold = ',IRPoleCheckThreshold
501+ write(*,*) ' > PrecisionVirtualAtRunTime = ',PrecisionVirtualAtRunTime
502+ if (SelectedContributionTypes(0).gt.0) then
503+ write(*,*) ' > SelectedContributionTypes = ', &
504+ (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
505+ else
506+ write(*,*) ' > SelectedContributionTypes = All'
507+ endif
508+ if (VetoedContributionTypes(0).gt.0) then
509+ write(*,*) ' > VetoedContributionTypes = ', &
510+ (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
511+ else
512+ write(*,*) ' > VetoedContributionTypes = None'
513+ endif
514+ if (QCD_squared_selected.eq.-1) then
515+ write(*,*) ' > QCD_squared_selected = All'
516+ else
517+ write(*,*) ' > QCD_squared_selected = ',QCD_squared_selected
518+ endif
519+ if (QED_squared_selected.eq.-1) then
520+ write(*,*) ' > QED_squared_selected = All'
521+ else
522+ write(*,*) ' > QED_squared_selected = ',QED_squared_selected
523+ endif
524+ if (SelectedCouplingOrders(1,0).gt.0) then
525+ do j=1,SelectedCouplingOrders(1,0)
526+ write(*,*) ' > SelectedCouplingOrders(',j,') = ', &
527+ (SelectedCouplingOrders(i,j),i=1,nsplitorders)
528+ enddo
529+ else
530+ write(*,*) ' > SelectedCouplingOrders = All'
531+ endif
532+ write(*,*) ' > NHelForMCoverHels = ',NHelForMCoverHels
533+ write(*,*) ' > VirtualFraction = ',Virt_fraction
534+ write(*,*) ' > MinVirtualFraction = ',Min_virt_fraction
535+ write(*,*) ' > SeparateFlavourConfigs = ',separate_flavour_configs
536+ write(*,*) ' > UsePolyVirtual = ',use_poly_virtual
537+ write(*,*) &
538+ '==============================================================='
539+ paramPrinted=.TRUE.
540+ endif
541+
542+ close(68)
543+ HasReadOnce=.TRUE.
544+901 continue
545+ end subroutine FKSParamReader
546+
547+ subroutine DefaultFKSParam()
548+ ! Sets the default parameters
549+ implicit none
550+ integer i,j
551+ IRPoleCheckThreshold=1.0d-5
552+ NHelForMCoverHels=5
553+ PrecisionVirtualAtRunTime=1d-3
554+ Virt_fraction=1d0
555+ QED_squared_selected=-1
556+ QCD_squared_selected=-1
557+ Min_virt_fraction=0.005d0
558+ separate_flavour_configs=.false.
559+ use_poly_virtual=.true.
560+ IncludeBornContributions=.true.
561+ SelectedContributionTypes(0)=0
562+ VetoedContributionTypes(0)=0
563+ do i=1, maxContribsSelected
564+ SelectedContributionTypes(I)=-1
565+ VetoedContributionTypes(I)=-1
566+ enddo
567+ do j=1,maxCouplingTypes
568+ SelectedCouplingOrders(j,0) = 0
569+ enddo
570+ do j=1,maxCouplingsSelected
571+ do i=1,maxCouplingTypes
572+ SelectedCouplingOrders(i,j) = -1
573+ enddo
574+ enddo
575+ end subroutine DefaultFKSParam
576+
577+end module FKSParams
578
579=== removed file 'Template/NLO/SubProcesses/FKSParams.inc'
580--- Template/NLO/SubProcesses/FKSParams.inc 2017-08-04 19:52:15 +0000
581+++ Template/NLO/SubProcesses/FKSParams.inc 1970-01-01 00:00:00 +0000
582@@ -1,31 +0,0 @@
583-!====================================================================
584-!
585-! Define common block with all general parameters used by MadFKS
586-! See their definitions in the file FKS_params.dat.
587-!
588-!====================================================================
589-
590- character*64 paramFileName
591- parameter ( paramFileName='FKS_params.dat')
592-
593- integer maxContribsSelected, maxCouplingsSelected
594- parameter( maxContribsSelected=100, maxCouplingsSelected=100 )
595- integer maxContribType, maxCouplingTypes
596- parameter ( maxContribType=15, maxCouplingTypes=20 )
597-
598- real*8 IRPoleCheckThreshold,Virt_fraction,
599- & PrecisionVirtualAtRunTime,Min_virt_fraction
600-
601- integer NHelForMCoverHels, VetoedContributionTypes(0:maxContribsSelected),
602- & SelectedContributionTypes(0:maxContribsSelected),
603- & QED_squared_selected, QCD_squared_selected,
604- & SelectedCouplingOrders(maxCouplingTypes,0:maxCouplingsSelected)
605-
606- logical separate_flavour_configs, IncludeBornContributions
607-
608- common /FKSPARAMS/IRPoleCheckThreshold,virt_fraction,
609- & PrecisionVirtualAtRunTime,min_virt_fraction,
610- & NHelForMCoverHels,separate_flavour_configs,
611- & QCD_squared_selected, QED_squared_selected,
612- & VetoedContributionTypes, SelectedContributionTypes,
613- & SelectedCouplingOrders
614
615=== modified file 'Template/NLO/SubProcesses/check_poles.f'
616--- Template/NLO/SubProcesses/check_poles.f 2018-12-21 15:27:32 +0000
617+++ Template/NLO/SubProcesses/check_poles.f 2019-09-12 11:14:33 +0000
618@@ -3,6 +3,7 @@
619 c This is the driver for the whole calulation
620 c**************************************************************************
621 use mint_module
622+ use FKSParams
623 implicit none
624 C
625 C CONSTANTS
626@@ -64,7 +65,6 @@
627 common /to_polecheck/force_polecheck, polecheck_passed
628 integer ret_code_ml
629 common /to_ret_code/ret_code_ml
630- include 'FKSParams.inc'
631
632 C-----
633 C BEGIN CODE
634
635=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
636--- Template/NLO/SubProcesses/driver_mintFO.f 2019-09-12 09:03:04 +0000
637+++ Template/NLO/SubProcesses/driver_mintFO.f 2019-09-12 11:14:33 +0000
638@@ -4,6 +4,7 @@
639 c**************************************************************************
640 use extra_weights
641 use mint_module
642+ use FKSParams
643 implicit none
644 C
645 C CONSTANTS
646@@ -66,9 +67,6 @@
647 include "timing_variables.inc"
648 real*4 tOther, tTot
649
650-c general MadFKS parameters
651- include "FKSParams.inc"
652-
653 c applgrid
654 integer iappl
655 common /for_applgrid/ iappl
656
657=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
658--- Template/NLO/SubProcesses/driver_mintMC.f 2019-09-12 09:03:04 +0000
659+++ Template/NLO/SubProcesses/driver_mintMC.f 2019-09-12 11:14:33 +0000
660@@ -4,6 +4,7 @@
661 c**************************************************************************
662 use extra_weights
663 use mint_module
664+ use FKSParams
665 implicit none
666 C
667 C CONSTANTS
668@@ -74,9 +75,6 @@
669 include "timing_variables.inc"
670 real*4 tOther, tTot
671
672-c general MadFKS parameters
673- include "FKSParams.inc"
674-
675 double precision deravg,derstd,dermax,xi_i_fks_ev_der_max
676 & ,y_ij_fks_ev_der_max
677 integer ntot_granny,derntot,ncase(0:6)
678
679=== modified file 'Template/NLO/SubProcesses/fks_singular.f'
680--- Template/NLO/SubProcesses/fks_singular.f 2019-09-12 09:03:04 +0000
681+++ Template/NLO/SubProcesses/fks_singular.f 2019-09-12 11:14:33 +0000
682@@ -1456,6 +1456,7 @@
683 c contribution
684 use weight_lines
685 use extra_weights
686+ use FKSParams
687 implicit none
688 include 'nexternal.inc'
689 include 'run.inc'
690@@ -1463,7 +1464,6 @@
691 include 'coupl.inc'
692 include 'fks_info.inc'
693 include 'q_es.inc'
694- include 'FKSParams.inc'
695 include 'orders.inc'
696 integer type,i,j
697 logical foundIt,foundOrders
698@@ -1651,6 +1651,7 @@
699 use weight_lines
700 use extra_weights
701 use mint_module
702+ use FKSParams
703 implicit none
704 include 'nexternal.inc'
705 include 'run.inc'
706@@ -1658,7 +1659,6 @@
707 include 'timing_variables.inc'
708 include 'genps.inc'
709 include 'orders.inc'
710- include 'FKSParams.inc'
711 integer orders(nsplitorders)
712 integer i,j,k,iamp,icontr_orig
713 logical virt_found
714@@ -2011,11 +2011,11 @@
715 c wgts() array to include the weights.
716 use weight_lines
717 use extra_weights
718+ use FKSParams
719 implicit none
720 include 'nexternal.inc'
721 include 'run.inc'
722 include 'timing_variables.inc'
723- include 'FKSParams.inc'
724 include 'genps.inc'
725 integer i,kr,kf,iwgt_save,dd
726 double precision xlum(maxscales),dlum,pi,mu2_r(maxscales),c_mu2_r
727@@ -2094,11 +2094,11 @@
728 c computations (ickkw.eq.-1).
729 use weight_lines
730 use extra_weights
731+ use FKSParams
732 implicit none
733 include 'nexternal.inc'
734 include 'run.inc'
735 include 'timing_variables.inc'
736- include 'FKSParams.inc'
737 include 'genps.inc'
738 integer i,ks,kh,iwgt_save
739 double precision xlum(maxscales),dlum,pi,mu2_r(maxscales)
740@@ -2200,11 +2200,11 @@
741 c wgts() array to include the weights.
742 use weight_lines
743 use extra_weights
744+ use FKSParams
745 implicit none
746 include 'nexternal.inc'
747 include 'run.inc'
748 include 'timing_variables.inc'
749- include 'FKSParams.inc'
750 include 'genps.inc'
751 integer n,izero,i,nn
752 parameter (izero=0)
753@@ -5877,12 +5877,22 @@
754 tOLP=tOLP+(tAfter-tBefore)
755 virtual_over_born=virt_wgt/born_wgt
756 if (ickkw.ne.-1) then
757- virt_wgt=virt_wgt-average_virtual(0,ichan)*born_wgt
758+ if (use_poly_virtual) then
759+ virt_wgt=virt_wgt-polyfit(0)*born_wgt
760+ else
761+ virt_wgt=virt_wgt-average_virtual(0,ichan)*born_wgt
762+ endif
763 do iamp=1,amp_split_size
764 if (amp_split_virt(iamp).eq.0d0) cycle
765+ if (use_poly_virtual) then
766+ amp_split_virt(iamp)=amp_split_virt(iamp)-
767+ $ polyfit(iamp)
768+ $ *amp_split_born_for_virt(iamp)
769+ else
770 amp_split_virt(iamp)=amp_split_virt(iamp)-
771 $ average_virtual(iamp,ichan)
772 $ *amp_split_born_for_virt(iamp)
773+ endif
774 enddo
775 endif
776 if (abrv.ne.'virt') then
777@@ -5902,12 +5912,21 @@
778 $ amp_split_virt_save(1:amp_split_size)
779 endif
780 if (abrv(1:4).ne.'virt' .and. ickkw.ne.-1) then
781- avv_wgt=average_virtual(0,ichan)*born_wgt
782- do iamp=1, amp_split_size
783- if (amp_split_born_for_virt(iamp).eq.0d0) cycle
784- amp_split_avv(iamp)= average_virtual(iamp,ichan)
785- $ *amp_split_born_for_virt(iamp)
786- enddo
787+ if (use_poly_virtual) then
788+ avv_wgt=polyfit(0)*born_wgt
789+ do iamp=1, amp_split_size
790+ if (amp_split_born_for_virt(iamp).eq.0d0) cycle
791+ amp_split_avv(iamp)= polyfit(iamp)
792+ $ *amp_split_born_for_virt(iamp)
793+ enddo
794+ else
795+ avv_wgt=average_virtual(0,ichan)*born_wgt
796+ do iamp=1, amp_split_size
797+ if (amp_split_born_for_virt(iamp).eq.0d0) cycle
798+ amp_split_avv(iamp)= average_virtual(iamp,ichan)
799+ $ *amp_split_born_for_virt(iamp)
800+ enddo
801+ endif
802 endif
803
804 c eq.(MadFKS.C.13)
805
806=== modified file 'Template/NLO/SubProcesses/makefile_fks_dir'
807--- Template/NLO/SubProcesses/makefile_fks_dir 2018-12-21 15:27:32 +0000
808+++ Template/NLO/SubProcesses/makefile_fks_dir 2019-09-12 11:14:33 +0000
809@@ -26,16 +26,16 @@
810
811 # Files for all executables
812 FILES= $(patsubst %.f,%.o,$(wildcard parton_lum_*.f)) $(patsubst \
813- %.f,%.o,$(wildcard matrix_*.f)) real_me_chooser.o \
814+ %.f,%.o,$(wildcard matrix_*.f)) FKSParams.o real_me_chooser.o \
815 chooser_functions.o genps_fks.o setcuts.o setscales.o \
816 veto_xsec.o $(patsubst %.f,%.o,$(wildcard b_sf_???.f)) \
817 born.o sborn_sf.o extra_cnt_wrapper.o $(patsubst \
818 %.f,%.o,$(wildcard born_cnt_*.f)) fks_Sij.o \
819 $(fastjetfortran_madfks) fks_singular.o montecarlocounter.o \
820 reweight_xsec.o boostwdir2.o initcluster.o cluster.o \
821- splitorders_stuff.o reweight.o get_color.o FKSParamReader.o \
822+ splitorders_stuff.o reweight.o get_color.o \
823 iproc_map.o MC_integer.o $(reweight_xsec_events_pdf_dummy) \
824- $(applgrid_interface) weight_lines.o mint_module.o
825+ $(applgrid_interface) weight_lines.o mint_module.o polfit.o
826
827 # Files needed for mintFO & mintMC
828 RUN= $(FO_ANALYSE) $(FILES) cuts.o pythia_unlops.o recluster.o \
829
830=== modified file 'Template/NLO/SubProcesses/mint_module.f90'
831--- Template/NLO/SubProcesses/mint_module.f90 2019-08-15 13:57:46 +0000
832+++ Template/NLO/SubProcesses/mint_module.f90 2019-09-12 11:14:33 +0000
833@@ -64,8 +64,8 @@
834 ! nintegrals>6 : virtual and born order by order
835 !
836
837-
838 module mint_module
839+ use FKSParams ! contains use_poly_virtual
840 implicit none
841 integer, parameter, private :: nintervals=32 ! max number of intervals in the integration grids
842 integer, parameter, public :: ndimmax=60 ! max number of dimensions of the integral
843@@ -92,7 +92,7 @@
844 integer, dimension(maxchannels), public :: iconfigs
845 double precision, public :: accuracy,min_virt_fraction_mint,wgt_mult
846 double precision, dimension(0:n_ave_virt,maxchannels), public :: average_virtual
847- double precision, dimension(0:n_ave_virt), public :: virt_wgt_mint,born_wgt_mint
848+ double precision, dimension(0:n_ave_virt), public :: virt_wgt_mint,born_wgt_mint,polyfit
849 double precision, dimension(maxchannels), public :: virtual_fraction
850 double precision, dimension(nintegrals,0:maxchannels), public :: ans,unc
851 logical :: only_virt,new_point,pass_cuts_check
852@@ -340,9 +340,13 @@
853 ! overwrite xgrid with the new xgrid
854 if (regridded(kchan)) xgrid(1:nint_used,1:ndim,kchan)=xgrid_new(1:nint_used,1:ndim)
855 enddo
856- do k_ord_virt=0,n_ord_virt
857- call regrid_ave_virt(k_ord_virt)
858- enddo
859+ if (use_poly_virtual) then
860+ call do_polyfit()
861+ else
862+ do k_ord_virt=0,n_ord_virt
863+ call regrid_ave_virt(k_ord_virt)
864+ enddo
865+ endif
866 ! Regrid the MC over integers (used for the MC over FKS dirs)
867 call regrid_MC_integer
868 end subroutine update_integration_grids
869@@ -664,10 +668,18 @@
870 isix=2*k_ord_virt+6
871 endif
872 if (f(ithree).ne.0d0) then
873- f(isix)=f(isix)/virtual_fraction(ichan)
874- virtual=(f(ithree)+average_virtual(k_ord_virt,ichan)*f(isix))*virtual_fraction(ichan)
875- born=f(isix)*virtual_fraction(ichan)
876- call fill_ave_virt(x,k_ord_virt,virtual,born)
877+ born=f(isix)
878+ ! virt_wgt_mint=(virtual-average_virtual*born)/virtual_fraction. Compensate:
879+ if (use_poly_virtual) then
880+ virtual=f(ithree)*virtual_fraction(ichan)+ &
881+ polyfit(k_ord_virt)*f(isix)
882+ call add_point_polyfit(ichan,k_ord_virt,x(1:ndim-3), &
883+ virtual/born,born/wgt_mult)
884+ else
885+ virtual=f(ithree)*virtual_fraction(ichan)+ &
886+ average_virtual(k_ord_virt,ichan)*f(isix)
887+ call fill_ave_virt(x,k_ord_virt,virtual,born)
888+ endif
889 else
890 f(isix)=0d0
891 endif
892@@ -801,7 +813,11 @@
893 if(imode.eq.0) nhits(icell(kdim),kdim,ichan)=nhits(icell(kdim),kdim,ichan)+1
894 enddo
895 do k_ord_virt=0,n_ord_virt
896- call get_ave_virt(x,k_ord_virt)
897+ if (use_poly_virtual) then
898+ call get_polyfit(ichan,k_ord_virt,x(1:ndim-3),polyfit(k_ord_virt))
899+ else
900+ call get_ave_virt(x,k_ord_virt)
901+ endif
902 enddo
903 end subroutine get_random_x
904
905@@ -988,14 +1004,18 @@
906
907 subroutine reset_mint_grids
908 implicit none
909- integer :: kdim,kchan,kint,k_ord_virt
910+ integer :: kdim,kchan,kint
911 do kint=0,nint_used
912 xgrid(kint,1:ndim,1:nchans)=dble(kint)/nint_used
913 enddo
914 nhits(1:nint_used,1:ndim,1:nchans)=0
915 regridded(1:nchans)=.true.
916 nhits_in_grids(1:nchans)=0
917- call init_ave_virt
918+ if (use_poly_virtual) then
919+ call init_polyfit(ndim-3,nchans,n_ord_virt,1000)
920+ else
921+ call init_ave_virt
922+ endif
923 ans_chan(0:nchans)=0d0
924 if (double_events) then
925 ! when double events, start with the very first channel only. For the
926@@ -1032,11 +1052,13 @@
927 write (12,*) 'MAX',(ymax(j,i,kchan),i=1,ndim)
928 enddo
929 endif
930- do j=1,nintervals_virt
931- do k=0,n_ord_virt
932- write (12,*) 'AVE',(ave_virt(j,i,k,kchan),i=1,ndim)
933+ if (.not.use_poly_virtual) then
934+ do j=1,nintervals_virt
935+ do k=0,n_ord_virt
936+ write (12,*) 'AVE',(ave_virt(j,i,k,kchan),i=1,ndim)
937+ enddo
938 enddo
939- enddo
940+ endif
941 if (imode.ge.1) then
942 write (12,*) 'MAX',ymax_virt(kchan)
943 endif
944@@ -1046,6 +1068,7 @@
945 write (12,*) 'AVE',virtual_fraction(kchan),average_virtual(0,kchan)
946 enddo
947 write (12,*) 'IDE',(ifold(i),i=1,ndim)
948+ if (use_poly_virtual) call save_polyfit(12)
949 close (12)
950 end subroutine write_grids_to_file
951
952@@ -1053,6 +1076,7 @@
953 ! Read the MINT integration grids from file
954 implicit none
955 integer :: i,j,k,kchan,idum
956+ integer,dimension(maxchannels) :: points
957 character(len=3) :: dummy
958 open (unit=12, file='mint_grids',status='old')
959 ans(1,0)=0d0
960@@ -1066,11 +1090,13 @@
961 read (12,*) dummy,(ymax(j,i,kchan),i=1,ndim)
962 enddo
963 endif
964- do j=1,nintervals_virt
965- do k=0,n_ord_virt
966- read (12,*) dummy,(ave_virt(j,i,k,kchan),i=1,ndim)
967+ if (.not.use_poly_virtual) then
968+ do j=1,nintervals_virt
969+ do k=0,n_ord_virt
970+ read (12,*) dummy,(ave_virt(j,i,k,kchan),i=1,ndim)
971+ enddo
972 enddo
973- enddo
974+ endif
975 if (imode.ge.2) then
976 read (12,*) dummy,ymax_virt(kchan)
977 endif
978@@ -1083,6 +1109,18 @@
979 enddo
980 read (12,*) dummy,(ifold(i),i=1,ndim)
981 unc(1,0)=sqrt(unc(1,0))
982+ ! polyfit stuff:
983+ if (use_poly_virtual) then
984+ do kchan=1,nchans
985+ read (12,*) dummy,points(kchan)
986+ enddo
987+ do kchan=1,nchans
988+ backspace(12)
989+ enddo
990+ call init_polyfit(ndim-3,nchans,n_ord_virt,maxval(points(1:nchans)))
991+ call restore_polyfit(12)
992+ call do_polyfit()
993+ endif
994 close (12)
995 ! check for zero cross-section: if restoring grids corresponding to
996 ! sigma=0, just terminate the run
997
998=== added file 'Template/NLO/SubProcesses/polfit.f'
999--- Template/NLO/SubProcesses/polfit.f 1970-01-01 00:00:00 +0000
1000+++ Template/NLO/SubProcesses/polfit.f 2019-09-12 11:14:33 +0000
1001@@ -0,0 +1,720 @@
1002+ subroutine init_polyfit(ndims,nchans,k_ord_virt,npoints)
1003+! wrapper subroutine that saves all information that's needed for calls
1004+! to the polfit() and pvalue() subroutines. Since I was too lazy to
1005+! write a proper module, use fortran entry statements to jump to the
1006+! right locations in this subroutine.
1007+ implicit none
1008+ integer ndims,nchans,ndim,nchan,maxpoint,maxdeg,ierr,i,ichan,ic
1009+ $ ,iunit,npoints,n_ord_virt,k_ord_virt,ko
1010+ integer,allocatable,dimension(:) :: n
1011+ integer,allocatable,dimension(:,:,:) :: ndeg
1012+ real*8,allocatable,dimension(:) :: r,yp,a
1013+ real*8,allocatable,dimension(:,:,:) :: y,w,x2d,temp3
1014+ real*8,allocatable,dimension(:,:,:,:) :: a2d
1015+ logical,allocatable,dimension(:) :: valid_ord_virt
1016+ parameter (maxdeg=20)
1017+ real*8 absXS,virt,ave_fun,fun_at_x,eps,sum_w,yfit,x(*)
1018+ logical fit_done,verbose
1019+ parameter (verbose=.false.)
1020+ character(len=3) :: dummy
1021+ save
1022+ maxpoint=npoints
1023+ nchan=nchans
1024+ ndim=ndims
1025+ n_ord_virt=k_ord_virt
1026+ if (allocated(ndeg)) deallocate(ndeg)
1027+ allocate(ndeg(ndim,0:n_ord_virt,nchan))
1028+ if (allocated(n)) deallocate(n)
1029+ allocate(n(nchan))
1030+ if (allocated(a)) deallocate(a)
1031+ allocate(a(maxpoint*3+maxdeg*3+3))
1032+ if (allocated(y)) deallocate(y)
1033+ allocate(y(maxpoint,0:n_ord_virt,nchan))
1034+ if (allocated(w)) deallocate(w)
1035+ allocate(w(maxpoint,0:n_ord_virt,nchan))
1036+ if (allocated(r)) deallocate(r)
1037+ allocate(r(maxpoint))
1038+ if (allocated(yp)) deallocate(yp)
1039+ allocate(yp(maxdeg))
1040+ if (allocated(x2d)) deallocate(x2d)
1041+ allocate(x2d(maxpoint,ndim,nchan))
1042+ if (allocated(a2d)) deallocate(a2d)
1043+ allocate(a2d(maxdeg*3+3,ndim,0:n_ord_virt,nchan))
1044+ if (allocated(valid_ord_virt)) deallocate(valid_ord_virt)
1045+ allocate(valid_ord_virt(0:n_ord_virt))
1046+ n(1:nchan)=0
1047+ fit_done=.false.
1048+ valid_ord_virt(0:n_ord_virt)=.false.
1049+ return
1050+
1051+ entry add_point_polyfit(ichan,k_ord_virt,x,virt,absXS)
1052+ ! only increment 'n' and 'x2d' if it's the zeroth "order" in the
1053+ ! virtual/born ratio
1054+ valid_ord_virt(k_ord_virt)=.true.
1055+ if (k_ord_virt.eq.0) then
1056+ if (n(ichan).eq.maxpoint) then
1057+! increase the size of the allocated variables:
1058+ maxpoint=maxpoint+1000
1059+ allocate(temp3(maxpoint,0:n_ord_virt,nchan))
1060+ do i=1,nchan
1061+ temp3(1:n(i),0:n_ord_virt,i)=y(1:n(i),0:n_ord_virt,i)
1062+ enddo
1063+ call move_alloc(temp3,y)
1064+ allocate(temp3(maxpoint,0:n_ord_virt,nchan))
1065+ do i=1,nchan
1066+ temp3(1:n(i),0:n_ord_virt,i)=w(1:n(i),0:n_ord_virt,i)
1067+ enddo
1068+ call move_alloc(temp3,w)
1069+ deallocate(r)
1070+ allocate(r(maxpoint))
1071+ allocate(temp3(maxpoint,ndim,nchan))
1072+ do i=1,nchan
1073+ temp3(1:n(i),1:ndim,i)=x2d(1:n(i),1:ndim,i)
1074+ enddo
1075+ call move_alloc(temp3,x2d)
1076+ deallocate(a)
1077+ allocate(a(maxpoint*3+maxdeg*3+3))
1078+ endif
1079+ n(ichan)=n(ichan)+1
1080+! add the point to the saved information
1081+ do i=1,ndim
1082+ x2d(n(ichan),i,ichan)=x(i)
1083+ enddo
1084+! initialise all to zero
1085+ y(n(ichan),0:n_ord_virt,ichan)=0d0
1086+ w(n(ichan),0:n_ord_virt,ichan)=0d0
1087+ endif
1088+ y(n(ichan),k_ord_virt,ichan)=virt
1089+ w(n(ichan),k_ord_virt,ichan)=sqrt(abs(absXS))
1090+ return
1091+
1092+ entry do_polyfit()
1093+ ! perform the fit for each dimension and each channel
1094+ do ic=1,nchan
1095+ do ko=0,n_ord_virt
1096+ if (.not.valid_ord_virt(ko)) cycle
1097+ do i=1,ndim
1098+ eps=-1d0
1099+ call polfit(n(ic),x2d(1,i,ic),y(1,ko,ic),w(1,ko,ic)
1100+ $ ,maxdeg,ndeg(i,ko,ic),eps,r,ierr,a)
1101+ ! this is the information used by pvalue() to interpolate
1102+ ! later:
1103+ a2d(1:maxdeg*3+3,i,ko,ic)=a(1:maxdeg*3+3)
1104+ ! print something in the logs:
1105+ if (verbose) write (*,'(a,i3,i3,i3,i8,i3,i3,f10.4)')
1106+ $ 'polyfit -',ic,ko,i,n(ic),ndeg(i,ko,ic),ierr,eps
1107+ enddo
1108+ enddo
1109+ enddo
1110+ fit_done=.true.
1111+ return
1112+
1113+ entry get_polyfit(ichan,k_ord_virt,x,fun_at_x)
1114+ ! if we haven't do a fit to the data, just return zero. This is
1115+ ! typically the case only during the very first iteration.
1116+ if ((.not.fit_done) .or. (.not.valid_ord_virt(k_ord_virt))) then
1117+ fun_at_x=0d0
1118+ return
1119+ endif
1120+ ! Do the interpolation. Since we have several dimensions, we use
1121+ ! the following:
1122+ ! f(x) = average + sum_i (f_i(x_i)-average)
1123+ ! where 'i' loops over all dimensions. In other words, the
1124+ ! information from each of the separate dimensions should only add
1125+ ! (or subtract) something from the average value. To obtain the
1126+ ! average, simply call pvalue() to return the zeroth-order
1127+ ! polynomial.
1128+ fun_at_x=0d0
1129+ do i=1,ndim
1130+ a(1:maxdeg*3+3)=a2d(1:maxdeg*3+3,i,k_ord_virt,ichan)
1131+ call pvalue(ndeg(i,k_ord_virt,ichan),0,x(i),yfit,yp,a)
1132+ call pvalue(0,0,x(i),ave_fun,yp,a)
1133+ fun_at_x=fun_at_x+(yfit-ave_fun)+ave_fun/dble(ndim)
1134+ enddo
1135+ return
1136+
1137+ entry save_polyfit(iunit)
1138+ do ic=1,nchan
1139+ write (iunit,*) 'POL',n(ic)
1140+ enddo
1141+ do ic=1,nchan
1142+ do i=1,n(ic)
1143+ write (iunit,*) 'POL',x2d(i,1:ndim,ic),
1144+ & y(i,0:n_ord_virt,ic),w(i,0:n_ord_virt,ic)
1145+ enddo
1146+ enddo
1147+ return
1148+
1149+ entry restore_polyfit(iunit)
1150+ do ic=1,nchan
1151+ read (iunit,*) dummy,n(ic)
1152+ enddo
1153+ do ic=1,nchan
1154+ do i=1,n(ic)
1155+ read (iunit,*) dummy,x2d(i,1:ndim,ic)
1156+ & ,y(i,0:n_ord_virt,ic),w(i,0:n_ord_virt,ic)
1157+ enddo
1158+ enddo
1159+ do ic=1,nchan
1160+ do ko=0,n_ord_virt
1161+ do i=1,n(ic)
1162+ if (y(i,ko,ic).ne.0d0) valid_ord_virt(ko)=.true.
1163+ enddo
1164+ enddo
1165+ enddo
1166+ if (verbose) write (*,*) 'polyfit - valid orders virt'
1167+ $ ,valid_ord_virt(0:n_ord_virt)
1168+ return
1169+ end
1170+
1171+
1172+
1173+
1174+
1175+*DECK POLFIT
1176+ SUBROUTINE POLFIT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A)
1177+c
1178+c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179+cRF: make double precision and replace DO-loops and arirthmic
1180+c IF-statements to comply more modern standards. Since we made it
1181+c double precision, the 'extended, partial-double' precision part has
1182+c become useless. Could remove it, but there is no real reason for
1183+c that (should only be marginally faster).
1184+c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1185+c
1186+C***BEGIN PROLOGUE POLFIT
1187+C***PURPOSE Fit discrete data in a least squares sense by polynomials
1188+C in one variable.
1189+C***LIBRARY SLATEC
1190+C***CATEGORY K1A1A2
1191+C***TYPE SINGLE PRECISION (POLFIT-S, DPOLFT-D)
1192+C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
1193+C***AUTHOR Shampine, L. F., (SNLA)
1194+C Davenport, S. M., (SNLA)
1195+C Huddleston, R. E., (SNLL)
1196+C***DESCRIPTION
1197+C
1198+C Abstract
1199+C
1200+C Given a collection of points X(I) and a set of values Y(I) which
1201+C correspond to some function or measurement at each of the X(I),
1202+C subroutine POLFIT computes the weighted least-squares polynomial
1203+C fits of all degrees up to some degree either specified by the user
1204+C or determined by the routine. The fits thus obtained are in
1205+C orthogonal polynomial form. Subroutine PVALUE may then be
1206+C called to evaluate the fitted polynomials and any of their
1207+C derivatives at any point. The subroutine PCOEF may be used to
1208+C express the polynomial fits as powers of (X-C) for any specified
1209+C point C.
1210+C
1211+C The parameters for POLFIT are
1212+C
1213+C Input --
1214+C N - the number of data points. The arrays X, Y and W
1215+C must be dimensioned at least N (N .GE. 1).
1216+C X - array of values of the independent variable. These
1217+C values may appear in any order and need not all be
1218+C distinct.
1219+C Y - array of corresponding function values.
1220+C W - array of positive values to be used as weights. If
1221+C W(1) is negative, POLFIT will set all the weights
1222+C to 1.0, which means unweighted least squares error
1223+C will be minimized. To minimize relative error, the
1224+C user should set the weights to: W(I) = 1.0/Y(I)**2,
1225+C I = 1,...,N .
1226+C MAXDEG - maximum degree to be allowed for polynomial fit.
1227+C MAXDEG may be any non-negative integer less than N.
1228+C Note -- MAXDEG cannot be equal to N-1 when a
1229+C statistical test is to be used for degree selection,
1230+C i.e., when input value of EPS is negative.
1231+C EPS - specifies the criterion to be used in determining
1232+C the degree of fit to be computed.
1233+C (1) If EPS is input negative, POLFIT chooses the
1234+C degree based on a statistical F test of
1235+C significance. One of three possible
1236+C significance levels will be used: .01, .05 or
1237+C .10. If EPS=-1.0 , the routine will
1238+C automatically select one of these levels based
1239+C on the number of data points and the maximum
1240+C degree to be considered. If EPS is input as
1241+C -.01, -.05, or -.10, a significance level of
1242+C .01, .05, or .10, respectively, will be used.
1243+C (2) If EPS is set to 0., POLFIT computes the
1244+C polynomials of degrees 0 through MAXDEG .
1245+C (3) If EPS is input positive, EPS is the RMS
1246+C error tolerance which must be satisfied by the
1247+C fitted polynomial. POLFIT will increase the
1248+C degree of fit until this criterion is met or
1249+C until the maximum degree is reached.
1250+C
1251+C Output --
1252+C NDEG - degree of the highest degree fit computed.
1253+C EPS - RMS error of the polynomial of degree NDEG .
1254+C R - vector of dimension at least NDEG containing values
1255+C of the fit of degree NDEG at each of the X(I) .
1256+C Except when the statistical test is used, these
1257+C values are more accurate than results from subroutine
1258+C PVALUE normally are.
1259+C IERR - error flag with the following possible values.
1260+C 1 -- indicates normal execution, i.e., either
1261+C (1) the input value of EPS was negative, and the
1262+C computed polynomial fit of degree NDEG
1263+C satisfies the specified F test, or
1264+C (2) the input value of EPS was 0., and the fits of
1265+C all degrees up to MAXDEG are complete, or
1266+C (3) the input value of EPS was positive, and the
1267+C polynomial of degree NDEG satisfies the RMS
1268+C error requirement.
1269+C 2 -- invalid input parameter. At least one of the input
1270+C parameters has an illegal value and must be corrected
1271+C before POLFIT can proceed. Valid input results
1272+C when the following restrictions are observed
1273+C N .GE. 1
1274+C 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0.
1275+C 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0.
1276+C W(1)=-1.0 or W(I) .GT. 0., I=1,...,N .
1277+C 3 -- cannot satisfy the RMS error requirement with a
1278+C polynomial of degree no greater than MAXDEG . Best
1279+C fit found is of degree MAXDEG .
1280+C 4 -- cannot satisfy the test for significance using
1281+C current value of MAXDEG . Statistically, the
1282+C best fit found is of order NORD . (In this case,
1283+C NDEG will have one of the values: MAXDEG-2,
1284+C MAXDEG-1, or MAXDEG). Using a higher value of
1285+C MAXDEG may result in passing the test.
1286+C A - work and output array having at least 3N+3MAXDEG+3
1287+C locations
1288+C
1289+C Note - POLFIT calculates all fits of degrees up to and including
1290+C NDEG . Any or all of these fits can be evaluated or
1291+C expressed as powers of (X-C) using PVALUE and PCOEF
1292+C after just one call to POLFIT .
1293+C
1294+C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston,
1295+C Curve fitting by polynomials in one variable, Report
1296+C SLA-74-0270, Sandia Laboratories, June 1974.
1297+C***ROUTINES CALLED PVALUE, XERMSG
1298+C***REVISION HISTORY (YYMMDD)
1299+C 740601 DATE WRITTEN
1300+C 890531 Changed all specific intrinsics to generic. (WRB)
1301+C 890531 REVISION DATE from Version 3.2
1302+C 891214 Prologue converted to Version 4.0 format. (BAB)
1303+C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
1304+C 920501 Reformatted the REFERENCES section. (WRB)
1305+C 920527 Corrected erroneous statements in DESCRIPTION. (WRB)
1306+C*** END PROLOGUE POLFIT
1307+ IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1308+ DOUBLE PRECISION TEMD1,TEMD2
1309+ DIMENSION X(*), Y(*), W(*), R(*), A(*), YP(0)
1310+ DIMENSION CO(4,3)
1311+ SAVE CO
1312+ DATA CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2),
1313+ 1 CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3),
1314+ 2 CO(4,3)/-13.086850,-2.4648165,-3.3846535,-1.2973162,
1315+ 3 -3.3381146,-1.7812271,-3.2578406,-1.6589279,
1316+ 4 -1.6282703,-1.3152745,-3.2640179,-1.9829776/
1317+C*** FIRST EXECUTABLE STATEMENT POLFIT
1318+ M = ABS(N)
1319+ IF (M .EQ. 0) GO TO 30
1320+ IF (MAXDEG .LT. 0) GO TO 30
1321+ A(1) = MAXDEG
1322+ MOP1 = MAXDEG + 1
1323+ IF (M .LT. MOP1) GO TO 30
1324+ IF (EPS .LT. 0.0 .AND. M .EQ. MOP1) GO TO 30
1325+ XM = M
1326+ ETST = EPS*EPS*XM
1327+ IF (W(1) .LT. 0.0) GO TO 2
1328+ DO I = 1,M
1329+ IF (W(I) .LE. 0.0) GO TO 30
1330+ enddo
1331+ GO TO 4
1332+ 2 DO I = 1,M
1333+ W(I) = 1.0
1334+ enddo
1335+ 4 IF (EPS .GE. 0.0) GO TO 8
1336+C
1337+C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR
1338+C CHOOSING DEGREE OF POLYNOMIAL FIT
1339+C
1340+ IF (EPS .GT. (-.55)) GO TO 5
1341+ IDEGF = M - MAXDEG - 1
1342+ KSIG = 1
1343+ IF (IDEGF .LT. 10) KSIG = 2
1344+ IF (IDEGF .LT. 5) KSIG = 3
1345+ GO TO 8
1346+ 5 KSIG = 1
1347+ IF (EPS .LT. (-.03)) KSIG = 2
1348+ IF (EPS .LT. (-.07)) KSIG = 3
1349+C
1350+C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING
1351+C
1352+ 8 K1 = MAXDEG + 1
1353+ K2 = K1 + MAXDEG
1354+ K3 = K2 + MAXDEG + 2
1355+ K4 = K3 + M
1356+ K5 = K4 + M
1357+ DO I = 2,K4
1358+ A(I) = 0.0
1359+ enddo
1360+ W11 = 0.0
1361+ IF (N .LT. 0) GO TO 11
1362+C
1363+C UNCONSTRAINED CASE
1364+C
1365+ DO I = 1,M
1366+ K4PI = K4 + I
1367+ A(K4PI) = 1.0
1368+ W11 = W11 + W(I)
1369+ enddo
1370+ GO TO 13
1371+C
1372+C CONSTRAINED CASE
1373+C
1374+ 11 DO I = 1,M
1375+ K4PI = K4 + I
1376+ W11 = W11 + W(I)*A(K4PI)**2
1377+ enddo
1378+C
1379+C COMPUTE FIT OF DEGREE ZERO
1380+C
1381+ 13 TEMD1 = 0.0D0
1382+ DO I = 1,M
1383+ K4PI = K4 + I
1384+ TEMD1 = TEMD1 + DBLE(W(I))*DBLE(Y(I))*DBLE(A(K4PI))
1385+ enddo
1386+ TEMD1 = TEMD1/DBLE(W11)
1387+ A(K2+1) = TEMD1
1388+ SIGJ = 0.0
1389+ DO I = 1,M
1390+ K4PI = K4 + I
1391+ K5PI = K5 + I
1392+ TEMD2 = TEMD1*DBLE(A(K4PI))
1393+ R(I) = TEMD2
1394+ A(K5PI) = TEMD2 - DBLE(R(I))
1395+ SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
1396+ enddo
1397+ J = 0
1398+C
1399+C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION
1400+C
1401+ if (eps.lt.0d0) then
1402+ goto 24
1403+ elseif (eps.eq.0d0) then
1404+ goto 26
1405+ else
1406+ goto 27
1407+ endif
1408+
1409+C
1410+C INCREMENT DEGREE
1411+C
1412+ 16 J = J + 1
1413+ JP1 = J + 1
1414+ K1PJ = K1 + J
1415+ K2PJ = K2 + J
1416+ SIGJM1 = SIGJ
1417+C
1418+C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1
1419+C
1420+ IF (J .GT. 1) A(K1PJ) = W11/W1
1421+C
1422+C COMPUTE NEW A COEFFICIENT
1423+C
1424+ TEMD1 = 0.0D0
1425+ DO I = 1,M
1426+ K4PI = K4 + I
1427+ TEMD2 = A(K4PI)
1428+ TEMD1 = TEMD1 + DBLE(X(I))*DBLE(W(I))*TEMD2*TEMD2
1429+ enddo
1430+ A(JP1) = TEMD1/DBLE(W11)
1431+C
1432+C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS
1433+C
1434+ W1 = W11
1435+ W11 = 0.0
1436+ DO I = 1,M
1437+ K3PI = K3 + I
1438+ K4PI = K4 + I
1439+ TEMP = A(K3PI)
1440+ A(K3PI) = A(K4PI)
1441+ A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP
1442+ W11 = W11 + W(I)*A(K4PI)**2
1443+ enddo
1444+C
1445+C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE
1446+C PRECISION
1447+C
1448+ TEMD1 = 0.0D0
1449+ DO I = 1,M
1450+ K4PI = K4 + I
1451+ K5PI = K5 + I
1452+ TEMD2 = DBLE(W(I))*DBLE((Y(I)-R(I))-A(K5PI))*DBLE(A(K4PI))
1453+ TEMD1 = TEMD1 + TEMD2
1454+ enddo
1455+ TEMD1 = TEMD1/DBLE(W11)
1456+ A(K2PJ+1) = TEMD1
1457+C
1458+C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND
1459+C ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE
1460+C COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT,
1461+C THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST
1462+C SIGNIFICANT BITS ARE IN A(K5PI) .
1463+C
1464+ SIGJ = 0.0
1465+ DO I = 1,M
1466+ K4PI = K4 + I
1467+ K5PI = K5 + I
1468+ TEMD2 = DBLE(R(I)) + DBLE(A(K5PI)) + TEMD1*DBLE(A(K4PI))
1469+ R(I) = TEMD2
1470+ A(K5PI) = TEMD2 - DBLE(R(I))
1471+ SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
1472+ enddo
1473+C
1474+C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE
1475+C MAXDEG HAS BEEN REACHED
1476+C
1477+ if (eps.lt.0d0) then
1478+ goto 23
1479+ elseif (eps.eq.0d0) then
1480+ goto 26
1481+ else
1482+ goto 27
1483+ endif
1484+C
1485+C COMPUTE F STATISTICS (INPUT EPS .LT. 0.)
1486+C
1487+ 23 IF (SIGJ .EQ. 0.0) GO TO 29
1488+ DEGF = M - J - 1
1489+ DEN = (CO(4,KSIG)*DEGF + 1.0)*DEGF
1490+ FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN
1491+ FCRIT = FCRIT*FCRIT
1492+ F = (SIGJM1 - SIGJ)*DEGF/SIGJ
1493+ IF (F .LT. FCRIT) GO TO 25
1494+C
1495+C POLYNOMIAL OF DEGREE J SATISFIES F TEST
1496+C
1497+ 24 SIGPAS = SIGJ
1498+ JPAS = J
1499+ NFAIL = 0
1500+ IF (MAXDEG .EQ. J) GO TO 32
1501+ GO TO 16
1502+C
1503+C POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE
1504+C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND.
1505+C
1506+ 25 NFAIL = NFAIL + 1
1507+ IF (NFAIL .GE. 3) GO TO 29
1508+ IF (MAXDEG .EQ. J) GO TO 32
1509+ GO TO 16
1510+C
1511+C RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT
1512+C EPS = 0.)
1513+C
1514+ 26 IF (MAXDEG .EQ. J) GO TO 28
1515+ GO TO 16
1516+C
1517+C SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.)
1518+C
1519+ 27 IF (SIGJ .LE. ETST) GO TO 28
1520+ IF (MAXDEG .EQ. J) GO TO 31
1521+ GO TO 16
1522+C
1523+C RETURNS
1524+C
1525+ 28 IERR = 1
1526+ NDEG = J
1527+ SIG = SIGJ
1528+ GO TO 33
1529+ 29 IERR = 1
1530+ NDEG = JPAS
1531+ SIG = SIGPAS
1532+ GO TO 33
1533+ 30 IERR = 2
1534+ write (*,*) 'SLATEC', 'POLFIT', 'INVALID INPUT PARAMETER.', 2,
1535+ + 1
1536+ stop 1
1537+ GO TO 37
1538+ 31 IERR = 3
1539+ NDEG = MAXDEG
1540+ SIG = SIGJ
1541+ GO TO 33
1542+ 32 IERR = 4
1543+ NDEG = JPAS
1544+ SIG = SIGPAS
1545+C
1546+ 33 A(K3) = NDEG
1547+C
1548+C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT
1549+C ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES
1550+C
1551+ IF(EPS .GE. 0.0 .OR. NDEG .EQ. MAXDEG) GO TO 36
1552+ NDER = 0
1553+ DO I = 1,M
1554+ CALL PVALUE (NDEG,NDER,X(I),R(I),YP,A)
1555+ enddo
1556+ 36 EPS = SQRT(SIG/XM)
1557+ 37 RETURN
1558+ END
1559+
1560+
1561+*DECK PVALUE
1562+ SUBROUTINE PVALUE (L, NDER, X, YFIT, YP, A)
1563+c
1564+cRF: make double precision and replace DO-loops to comply with more
1565+c modern standards
1566+c
1567+C***BEGIN PROLOGUE PVALUE
1568+C***PURPOSE Use the coefficients generated by POLFIT to evaluate the
1569+C polynomial fit of degree L, along with the first NDER of
1570+C its derivatives, at a specified point.
1571+C***LIBRARY SLATEC
1572+C***CATEGORY K6
1573+C***TYPE SINGLE PRECISION (PVALUE-S, DP1VLU-D)
1574+C***KEYWORDS CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION
1575+C***AUTHOR Shampine, L. F., (SNLA)
1576+C Davenport, S. M., (SNLA)
1577+C***DESCRIPTION
1578+C
1579+C Written by L. F. Shampine and S. M. Davenport.
1580+C
1581+C Abstract
1582+C
1583+C The subroutine PVALUE uses the coefficients generated by POLFIT
1584+C to evaluate the polynomial fit of degree L , along with the first
1585+C NDER of its derivatives, at a specified point. Computationally
1586+C stable recurrence relations are used to perform this task.
1587+C
1588+C The parameters for PVALUE are
1589+C
1590+C Input --
1591+C L - the degree of polynomial to be evaluated. L may be
1592+C any non-negative integer which is less than or equal
1593+C to NDEG , the highest degree polynomial provided
1594+C by POLFIT .
1595+C NDER - the number of derivatives to be evaluated. NDER
1596+C may be 0 or any positive value. If NDER is less
1597+C than 0, it will be treated as 0.
1598+C X - the argument at which the polynomial and its
1599+C derivatives are to be evaluated.
1600+C A - work and output array containing values from last
1601+C call to POLFIT .
1602+C
1603+C Output --
1604+C YFIT - value of the fitting polynomial of degree L at X
1605+C YP - array containing the first through NDER derivatives
1606+C of the polynomial of degree L . YP must be
1607+C dimensioned at least NDER in the calling program.
1608+C
1609+C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston,
1610+C Curve fitting by polynomials in one variable, Report
1611+C SLA-74-0270, Sandia Laboratories, June 1974.
1612+C***ROUTINES CALLED XERMSG
1613+C***REVISION HISTORY (YYMMDD)
1614+C 740601 DATE WRITTEN
1615+C 890531 Changed all specific intrinsics to generic. (WRB)
1616+C 890531 REVISION DATE from Version 3.2
1617+C 891214 Prologue converted to Version 4.0 format. (BAB)
1618+C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
1619+C 900510 Convert XERRWV calls to XERMSG calls. (RWC)
1620+C 920501 Reformatted the REFERENCES section. (WRB)
1621+C*** END PROLOGUE PVALUE
1622+ IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1623+ DIMENSION YP(*),A(*)
1624+ CHARACTER*8 XERN1, XERN2
1625+C***FIRST EXECUTABLE STATEMENT PVALUE
1626+ IF (L .LT. 0) GO TO 12
1627+ NDO = MAX(NDER,0)
1628+ NDO = MIN(NDO,L)
1629+ MAXORD = A(1) + 0.5
1630+ K1 = MAXORD + 1
1631+ K2 = K1 + MAXORD
1632+ K3 = K2 + MAXORD + 2
1633+ NORD = A(K3) + 0.5
1634+ IF (L .GT. NORD) GO TO 11
1635+ K4 = K3 + L + 1
1636+ IF (NDER .LT. 1) GO TO 2
1637+ DO I = 1,NDER
1638+ YP(I) = 0.0
1639+ enddo
1640+ 2 IF (L .GE. 2) GO TO 4
1641+ IF (L .EQ. 1) GO TO 3
1642+C
1643+C L IS 0
1644+C
1645+ VAL = A(K2+1)
1646+ GO TO 10
1647+C
1648+C L IS 1
1649+C
1650+ 3 CC = A(K2+2)
1651+ VAL = A(K2+1) + (X-A(2))*CC
1652+ IF (NDER .GE. 1) YP(1) = CC
1653+ GO TO 10
1654+C
1655+C L IS GREATER THAN 1
1656+C
1657+ 4 NDP1 = NDO + 1
1658+ K3P1 = K3 + 1
1659+ K4P1 = K4 + 1
1660+ LP1 = L + 1
1661+ LM1 = L - 1
1662+ ILO = K3 + 3
1663+ IUP = K4 + NDP1
1664+ DO I = ILO,IUP
1665+ A(I) = 0.0
1666+ enddo
1667+ DIF = X - A(LP1)
1668+ KC = K2 + LP1
1669+ A(K4P1) = A(KC)
1670+ A(K3P1) = A(KC-1) + DIF*A(K4P1)
1671+ A(K3+2) = A(K4P1)
1672+C
1673+C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES
1674+C
1675+ DO I = 1,LM1
1676+ IN = L - I
1677+ INP1 = IN + 1
1678+ K1I = K1 + INP1
1679+ IC = K2 + IN
1680+ DIF = X - A(INP1)
1681+ VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1)
1682+ IF (NDO .LE. 0) GO TO 8
1683+ DO N = 1,NDO
1684+ K3PN = K3P1 + N
1685+ K4PN = K4P1 + N
1686+ YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN)
1687+ enddo
1688+C
1689+C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS
1690+C
1691+ DO N = 1,NDO
1692+ K3PN = K3P1 + N
1693+ K4PN = K4P1 + N
1694+ A(K4PN) = A(K3PN)
1695+ A(K3PN) = YP(N)
1696+ enddo
1697+ 8 A(K4P1) = A(K3P1)
1698+ A(K3P1) = VAL
1699+ enddo
1700+C
1701+C NORMAL RETURN OR ABORT DUE TO ERROR
1702+C
1703+ 10 YFIT = VAL
1704+ RETURN
1705+C
1706+ 11 WRITE (XERN1, '(I8)') L
1707+ WRITE (XERN2, '(I8)') NORD
1708+ write (*,*) 'SLATEC', 'PVALUE',
1709+ * 'THE ORDER OF POLYNOMIAL EVALUATION, L = ' // XERN1 //
1710+ * ' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 //
1711+ * ', COMPUTED BY POLFIT -- EXECUTION TERMINATED.', 8, 2
1712+ stop 1
1713+ RETURN
1714+C
1715+ 12 write (*,*) 'SLATEC', 'PVALUE',
1716+ + 'INVALID INPUT PARAMETER. ORDER OF POLYNOMIAL EVALUATION ' //
1717+ + 'REQUESTED IS NEGATIVE -- EXECUTION TERMINATED.', 2, 2
1718+ stop 1
1719+ RETURN
1720+ END
1721+
1722
1723=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
1724--- madgraph/interface/amcatnlo_run_interface.py 2019-08-28 08:02:11 +0000
1725+++ madgraph/interface/amcatnlo_run_interface.py 2019-09-12 11:14:33 +0000
1726@@ -2376,9 +2376,12 @@
1727 # latter is only the only line that contains integers.
1728 for j,fs in enumerate([files_mint_grids,files_MC_integer]):
1729 linesoffiles=[]
1730+ polyfit_data=[]
1731 for f in fs:
1732 with open(f,'r+') as fi:
1733- linesoffiles.append(fi.readlines())
1734+ data=fi.readlines()
1735+ linesoffiles.append([ dat for dat in data if 'POL' not in dat.split()[0] ])
1736+ polyfit_data.append([ dat for dat in data if 'POL' in dat.split()[0] ])
1737 to_write=[]
1738 for rowgrp in zip(*linesoffiles):
1739 action=list(set([row.strip().split()[0] for row in rowgrp])) # list(set()) structure to remove duplicants
1740@@ -2414,8 +2417,30 @@
1741 write_string.append(int(sum(floatgrp)/len(floatgrp)))
1742 else:
1743 raise aMCatNLOError('Unknown action for combining grids: %s' % action[0])
1744-
1745 to_write.append(action[0] + " " + (" ".join(str(ws) for ws in write_string)) + "\n")
1746+
1747+ if polyfit_data:
1748+ # special for data regarding virtuals. Need to simply append
1749+ # all the data, but skipping doubles.
1750+ for i,onefile in enumerate(polyfit_data):
1751+ # Get the number of channels, and the number of PS points per channel
1752+ data_amount_in_file=[int(oneline.split()[1]) for oneline in onefile if len(oneline.split())==2]
1753+ if i==0:
1754+ filtered_list=[ [] for i in range(len(data_amount_in_file)) ]
1755+ start=len(data_amount_in_file)
1756+ for channel,channel_size in enumerate(data_amount_in_file):
1757+ end=start+channel_size
1758+ for data_channel in onefile[start:end]:
1759+ if data_channel not in filtered_list[channel]:
1760+ filtered_list[channel].append(data_channel)
1761+ start=end
1762+ # The amount of data in each file (per channel):
1763+ for channel in filtered_list:
1764+ to_write.append("POL " + str(len(channel)) + "\n")
1765+ # All the data:
1766+ for ch in filtered_list:
1767+ for dat in ch:
1768+ to_write.append(dat)
1769 # write the data over the master location
1770 if j==0:
1771 with open(pjoin(location,'mint_grids'),'w') as f:
1772
1773=== modified file 'madgraph/iolibs/export_fks.py'
1774--- madgraph/iolibs/export_fks.py 2019-08-15 13:09:08 +0000
1775+++ madgraph/iolibs/export_fks.py 2019-09-12 11:14:33 +0000
1776@@ -624,8 +624,7 @@
1777 'FKS_params.dat',
1778 'initial_states_map.dat',
1779 'OLE_order.olc',
1780- 'FKSParams.inc',
1781- 'FKSParamReader.f',
1782+ 'FKSParams.f90',
1783 'cuts.inc',
1784 'unlops.inc',
1785 'pythia_unlops.f',
1786@@ -688,7 +687,8 @@
1787 'randinit',
1788 'sudakov.inc',
1789 'maxconfigs.inc',
1790- 'timing_variables.inc']
1791+ 'timing_variables.inc',
1792+ 'polfit.f']
1793
1794 for file in linkfiles:
1795 ln('../' + file , '.')
1796@@ -1149,24 +1149,26 @@
1797
1798 amp_split_size=len(amp_split_orders)
1799
1800- text = 'C The orders to be integrated for the Born and at NLO\n'
1801+ text = '! The orders to be integrated for the Born and at NLO\n'
1802 text += 'integer nsplitorders\n'
1803 text += 'parameter (nsplitorders=%d)\n' % len(split_orders)
1804 text += 'character*3 ordernames(nsplitorders)\n'
1805 text += 'data ordernames / %s /\n' % ', '.join(['"%3s"' % o for o in split_orders])
1806 text += 'integer born_orders(nsplitorders), nlo_orders(nsplitorders)\n'
1807- text += 'C the order of the coupling orders is %s\n' % ', '.join(split_orders)
1808+ text += '! the order of the coupling orders is %s\n' % ', '.join(split_orders)
1809 text += 'data born_orders / %s /\n' % ', '.join([str(max_born_orders[o]) for o in split_orders])
1810 text += 'data nlo_orders / %s /\n' % ', '.join([str(max_nlo_orders[o]) for o in split_orders])
1811- text += 'C The position of the QCD /QED orders in the array\n'
1812+ text += '! The position of the QCD /QED orders in the array\n'
1813 text += 'integer qcd_pos, qed_pos\n'
1814- text += 'C if = -1, then it is not in the split_orders\n'
1815+ text += '! if = -1, then it is not in the split_orders\n'
1816 text += 'parameter (qcd_pos = %d)\n' % qcd_pos
1817 text += 'parameter (qed_pos = %d)\n' % qed_pos
1818- text += 'C this is to keep track of the various coupling combinations entering each ME\n'
1819+ text += '! this is to keep track of the various \n'
1820+ text += '! coupling combinations entering each ME\n'
1821 text += 'integer amp_split_size, amp_split_size_born\n'
1822 text += 'parameter (amp_split_size = %d)\n' % amp_split_size
1823- text += 'parameter (amp_split_size_born = %d) ! the first entries in amp_split are for the born\n' % amp_split_size_born
1824+ text += '! the first entries in the next line in amp_split are for the born \n'
1825+ text += 'parameter (amp_split_size_born = %d)\n' % amp_split_size_born
1826 text += 'double precision amp_split(amp_split_size)\n'
1827 text += 'double complex amp_split_cnt(amp_split_size,2,nsplitorders)\n'
1828 text += 'common /to_amp_split/amp_split, amp_split_cnt\n'

Subscribers

People subscribed via source and target branches

to all changes: