Page MenuHomeHEPForge
Files F7879532

No OneTemporary
Actions

View Options

Index: trunk/share/tests/pdf_builtin.sin
===================================================================
--- trunk/share/tests/pdf_builtin.sin (revision 5839)
+++ trunk/share/tests/pdf_builtin.sin (revision 5840)
@@ -1,72 +1,84 @@
# $Id: pdf_builtin.sin 4932 2013-12-05 08:46:27Z jr_reuter $
# SINDARIN input for WHIZARD self-test
# Process p, p -> e2, E2
library = "pdf_builtin_1_lib"
model = "SM"
?logging = true
?openmp_logging = false
?vis_history = false
?integration_timer = false
?phs_s_mapping = false
alias q = u:d:s
alias Q = U:D:S
alias lepton = e1:e2:e3:E1:E2:E3
seed = 0
$method = "omega"
process pdf_builtin_1_p = q, Q => e2, E2
process pdf_builtin_2_p = A, A => e1, E1
compile ()
ms = 0
me = 0
mmu = 0
$phs_method = "wood"
$integration_method = "vamp"
cuts = all Pt > 100 GeV [lepton]
and all M > 10 GeV [lepton, lepton]
sqrts = 1000
iterations = 4:500, 2:500
!!! Tests should be run single-threaded
openmp_num_threads = 1
beams = u, U
integrate (pdf_builtin_1_p)
beams = p, p => pdf_builtin
$pdf_builtin_set = "cteq6l"
integrate (pdf_builtin_1_p)
beams = p, p => pdf_builtin
$pdf_builtin_set = "mstw2008lo"
integrate (pdf_builtin_1_p)
beams = p, p => pdf_builtin
$pdf_builtin_set = "mstw2008nlo"
integrate (pdf_builtin_1_p)
beams = p, p => pdf_builtin
$pdf_builtin_set = "mstw2008nnlo"
integrate (pdf_builtin_1_p)
beams = p, p => pdf_builtin
$pdf_builtin_set = "ct10"
integrate (pdf_builtin_1_p)
+beams = p, p => pdf_builtin
+$pdf_builtin_set = "CJ12_max"
+integrate (pdf_builtin_1_p)
+
+beams = p, p => pdf_builtin
+$pdf_builtin_set = "CJ12_mid"
+integrate (pdf_builtin_1_p)
+
+beams = p, p => pdf_builtin
+$pdf_builtin_set = "CJ12_min"
+integrate (pdf_builtin_1_p)
+
beams = A, A
integrate (pdf_builtin_2_p)
beams = p, p => pdf_builtin
$pdf_builtin_set = "mrst2004qedp"
integrate (pdf_builtin_2_p)
Index: trunk/share/tests/ref-output-double/pdf_builtin.ref
===================================================================
--- trunk/share/tests/ref-output-double/pdf_builtin.ref (revision 5839)
+++ trunk/share/tests/ref-output-double/pdf_builtin.ref (revision 5840)
@@ -1,407 +1,551 @@
?openmp_logging = false
?vis_history = false
?integration_timer = false
?phs_s_mapping = false
[user variable] q = PDG(2, 1, 3)
[user variable] Q = PDG(-2, -1, -3)
[user variable] lepton = PDG(11, 13, 15, -11, -13, -15)
seed = 0
$method = "omega"
| Process library 'pdf_builtin_1_lib': recorded process 'pdf_builtin_1_p'
| Process library 'pdf_builtin_1_lib': recorded process 'pdf_builtin_2_p'
| Process library 'pdf_builtin_lib': compiling ...
| Process library 'pdf_builtin_lib': ... success.
| Process library 'pdf_builtin_1_lib': compiling ...
| Process library 'pdf_builtin_1_lib': writing makefile
| Process library 'pdf_builtin_1_lib': removing old files
| Process library 'pdf_builtin_1_lib': writing driver
| Process library 'pdf_builtin_1_lib': creating source code
| Process library 'pdf_builtin_1_lib': compiling sources
| Process library 'pdf_builtin_1_lib': linking
| Process library 'pdf_builtin_1_lib': loading
| Process library 'pdf_builtin_1_lib': ... success.
SM.ms = 0.000000000000E+00
SM.me = 0.000000000000E+00
SM.mmu = 0.000000000000E+00
$phs_method = "wood"
$integration_method = "vamp"
sqrts = 1.000000000000E+03
openmp_num_threads = 1
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 0
| Initializing integration for process pdf_builtin_1_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_1_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 1
| Process components:
| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
| ------------------------------------------------------------------------
| Beam structure: u, ubar
| Beam data (collision):
| u (mass = 0.0000000E+00 GeV)
| ubar (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
| Phase space: 2 channels, 2 dimensions
| Phase space: found 2 channels, collected in 2 groves.
| Phase space: Using 2 equivalences between channels.
| Phase space: wood
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_1_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 2 chains, 2 channels, 2 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 400 1.8321997E+01 2.98E-01 1.63 0.33* 35.80
2 400 1.8828751E+01 1.70E-01 0.90 0.18* 40.76
3 400 1.8682775E+01 1.30E-01 0.70 0.14* 53.32
4 400 1.8636496E+01 1.06E-01 0.57 0.11* 40.51
|-----------------------------------------------------------------------------|
4 1600 1.8666882E+01 7.19E-02 0.39 0.15 40.51 0.78 4
|-----------------------------------------------------------------------------|
5 400 1.8881752E+01 1.31E-01 0.69 0.14 40.70
6 400 1.8558177E+01 1.27E-01 0.68 0.14* 39.83
|-----------------------------------------------------------------------------|
6 800 1.8715083E+01 9.12E-02 0.49 0.14 39.83 3.14 2
|=============================================================================|
$pdf_builtin_set = "cteq6l"
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 1
| Initializing integration for process pdf_builtin_1_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_1_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 1
| Process components:
| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
| ------------------------------------------------------------------------
| Beam structure: p, p => pdf_builtin
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Initialized builtin PDF CTEQ6L
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
| Phase space: 2 channels, 2 dimensions
| Phase space: found 2 channels, collected in 2 groves.
| Phase space: Using 2 equivalences between channels.
| Phase space: wood
| Beam structure: pdf_builtin, none => none, pdf_builtin
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_1_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 2 chains, 2 channels, 4 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 500 1.0289901E+00 2.67E-01 25.95 5.80* 1.50
2 499 1.5389851E+00 1.45E-01 9.41 2.10* 7.00
3 498 1.6947237E+00 8.32E-02 4.91 1.10* 14.36
4 497 1.8264951E+00 9.45E-02 5.17 1.15 11.07
|-----------------------------------------------------------------------------|
4 1994 1.6884026E+00 5.61E-02 3.32 1.48 11.07 3.10 4
|-----------------------------------------------------------------------------|
5 497 1.7911732E+00 9.61E-02 5.36 1.20 7.49
6 497 1.6593299E+00 7.87E-02 4.74 1.06* 6.94
|-----------------------------------------------------------------------------|
6 994 1.7122737E+00 6.09E-02 3.56 1.12 6.94 1.13 2
|=============================================================================|
$pdf_builtin_set = "mstw2008lo"
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 2
| Initializing integration for process pdf_builtin_1_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_1_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 1
| Process components:
| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
| ------------------------------------------------------------------------
| Beam structure: p, p => pdf_builtin
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Initialized builtin PDF MSTW2008LO
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
| Phase space: 2 channels, 2 dimensions
| Phase space: found 2 channels, collected in 2 groves.
| Phase space: Using 2 equivalences between channels.
| Phase space: wood
| Beam structure: pdf_builtin, none => none, pdf_builtin
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_1_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 2 chains, 2 channels, 4 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 500 1.1880877E+00 5.55E-01 46.70 10.44* 0.75
2 499 1.7294095E+00 1.56E-01 9.03 2.02* 5.67
3 498 1.7268042E+00 9.93E-02 5.75 1.28* 10.16
4 497 1.7156684E+00 8.15E-02 4.75 1.06* 16.39
|-----------------------------------------------------------------------------|
4 1994 1.7155968E+00 5.81E-02 3.39 1.51 16.39 0.31 4
|-----------------------------------------------------------------------------|
5 497 1.8369997E+00 9.31E-02 5.07 1.13 11.18
6 497 1.7518523E+00 9.22E-02 5.26 1.17 10.67
|-----------------------------------------------------------------------------|
6 994 1.7939895E+00 6.55E-02 3.65 1.15 10.67 0.42 2
|=============================================================================|
$pdf_builtin_set = "mstw2008nlo"
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 3
| Initializing integration for process pdf_builtin_1_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_1_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 1
| Process components:
| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
| ------------------------------------------------------------------------
| Beam structure: p, p => pdf_builtin
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Initialized builtin PDF MSTW2008NLO
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
| Phase space: 2 channels, 2 dimensions
| Phase space: found 2 channels, collected in 2 groves.
| Phase space: Using 2 equivalences between channels.
| Phase space: wood
| Beam structure: pdf_builtin, none => none, pdf_builtin
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_1_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 2 chains, 2 channels, 4 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 500 1.4521330E+00 3.77E-01 25.96 5.81* 1.79
2 499 1.6921136E+00 1.63E-01 9.66 2.16* 5.64
3 498 1.6653622E+00 9.62E-02 5.78 1.29* 7.33
4 497 1.6425744E+00 1.07E-01 6.53 1.46 8.71
|-----------------------------------------------------------------------------|
4 1994 1.6550091E+00 6.46E-02 3.91 1.74 8.71 0.12 4
|-----------------------------------------------------------------------------|
5 497 1.5526409E+00 7.47E-02 4.81 1.07* 8.24
6 497 1.5904078E+00 7.61E-02 4.78 1.07* 8.44
|-----------------------------------------------------------------------------|
6 994 1.5711731E+00 5.33E-02 3.39 1.07 8.44 0.13 2
|=============================================================================|
$pdf_builtin_set = "mstw2008nnlo"
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 4
| Initializing integration for process pdf_builtin_1_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_1_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 1
| Process components:
| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
| ------------------------------------------------------------------------
| Beam structure: p, p => pdf_builtin
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Initialized builtin PDF MSTW2008NNLO
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
| Phase space: 2 channels, 2 dimensions
| Phase space: found 2 channels, collected in 2 groves.
| Phase space: Using 2 equivalences between channels.
| Phase space: wood
| Beam structure: pdf_builtin, none => none, pdf_builtin
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_1_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 2 chains, 2 channels, 4 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 500 2.2132140E+00 8.52E-01 38.50 8.61* 1.07
2 499 1.8511677E+00 1.49E-01 8.07 1.80* 7.67
3 498 1.5390988E+00 8.55E-02 5.56 1.24* 13.79
4 497 1.6992156E+00 1.48E-01 8.70 1.94 5.13
|-----------------------------------------------------------------------------|
4 1994 1.6363962E+00 6.61E-02 4.04 1.80 5.13 1.33 4
|-----------------------------------------------------------------------------|
5 497 1.6264232E+00 8.16E-02 5.01 1.12* 4.91
6 497 1.7014429E+00 1.20E-01 7.06 1.57 3.55
|-----------------------------------------------------------------------------|
6 994 1.6501017E+00 6.75E-02 4.09 1.29 3.55 0.27 2
|=============================================================================|
$pdf_builtin_set = "ct10"
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 5
| Initializing integration for process pdf_builtin_1_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_1_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 1
| Process components:
| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
| ------------------------------------------------------------------------
| Beam structure: p, p => pdf_builtin
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Initialized builtin PDF CT10
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
| Phase space: 2 channels, 2 dimensions
| Phase space: found 2 channels, collected in 2 groves.
| Phase space: Using 2 equivalences between channels.
| Phase space: wood
| Beam structure: pdf_builtin, none => none, pdf_builtin
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_1_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 2 chains, 2 channels, 4 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
1 500 1.3862192E+00 3.57E-01 25.78 5.76* 2.38
2 499 1.3203714E+00 1.15E-01 8.72 1.95* 8.18
3 498 1.4726542E+00 7.75E-02 5.26 1.17* 13.45
4 497 1.4785343E+00 8.04E-02 5.44 1.21 10.57
|-----------------------------------------------------------------------------|
4 1994 1.4448402E+00 4.97E-02 3.44 1.54 10.57 0.50 4
|-----------------------------------------------------------------------------|
5 497 1.5258851E+00 8.23E-02 5.39 1.20* 9.58
6 497 1.5517191E+00 8.03E-02 5.17 1.15* 8.72
|-----------------------------------------------------------------------------|
6 994 1.5391235E+00 5.75E-02 3.73 1.18 8.72 0.05 2
|=============================================================================|
+$pdf_builtin_set = "CJ12_max"
| RNG: Initializing TAO random-number generator
| RNG: Setting seed for random-number generator to 6
+| Initializing integration for process pdf_builtin_1_p:
+| ------------------------------------------------------------------------
+| Process [scattering]: 'pdf_builtin_1_p'
+| Library name = 'pdf_builtin_1_lib'
+| Process index = 1
+| Process components:
+| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
+| ------------------------------------------------------------------------
+| Beam structure: p, p => pdf_builtin
+| Beam data (collision):
+| p (mass = 0.0000000E+00 GeV)
+| p (mass = 0.0000000E+00 GeV)
+| sqrts = 1.000000000000E+03 GeV
+| Initialized builtin PDF CJ12_max
+| Phase space: generating configuration ...
+| Phase space: ... success.
+| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
+| Phase space: 2 channels, 2 dimensions
+| Phase space: found 2 channels, collected in 2 groves.
+| Phase space: Using 2 equivalences between channels.
+| Phase space: wood
+| Beam structure: pdf_builtin, none => none, pdf_builtin
+| Beam structure: 1 channels, 2 dimensions
+| Applying user-defined cuts.
+| Starting integration for process 'pdf_builtin_1_p'
+| Integrate: iterations = 4:500:"gw", 2:500
+| Integrator: 2 chains, 2 channels, 4 dimensions
+| Integrator: Using VAMP channel equivalences
+| Integrator: 500 initial calls, 20 bins, stratified = T
+| Integrator: VAMP
+|=============================================================================|
+| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
+|=============================================================================|
+ 1 500 1.2718554E+00 3.39E-01 26.68 5.97* 1.36
+ 2 499 9.4516438E-01 8.84E-02 9.35 2.09* 5.81
+ 3 498 1.1890105E+00 6.73E-02 5.66 1.26* 15.09
+ 4 497 1.2755777E+00 6.19E-02 4.86 1.08* 14.48
+|-----------------------------------------------------------------------------|
+ 4 1994 1.1762030E+00 4.02E-02 3.42 1.53 14.48 3.17 4
+|-----------------------------------------------------------------------------|
+ 5 497 1.2817716E+00 5.97E-02 4.65 1.04* 14.20
+ 6 497 1.2209063E+00 6.01E-02 4.92 1.10 12.63
+|-----------------------------------------------------------------------------|
+ 6 994 1.2515528E+00 4.23E-02 3.38 1.07 12.63 0.52 2
+|=============================================================================|
+$pdf_builtin_set = "CJ12_mid"
+| RNG: Initializing TAO random-number generator
+| RNG: Setting seed for random-number generator to 7
+| Initializing integration for process pdf_builtin_1_p:
+| ------------------------------------------------------------------------
+| Process [scattering]: 'pdf_builtin_1_p'
+| Library name = 'pdf_builtin_1_lib'
+| Process index = 1
+| Process components:
+| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
+| ------------------------------------------------------------------------
+| Beam structure: p, p => pdf_builtin
+| Beam data (collision):
+| p (mass = 0.0000000E+00 GeV)
+| p (mass = 0.0000000E+00 GeV)
+| sqrts = 1.000000000000E+03 GeV
+| Initialized builtin PDF CJ12_mid
+| Phase space: generating configuration ...
+| Phase space: ... success.
+| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
+| Phase space: 2 channels, 2 dimensions
+| Phase space: found 2 channels, collected in 2 groves.
+| Phase space: Using 2 equivalences between channels.
+| Phase space: wood
+| Beam structure: pdf_builtin, none => none, pdf_builtin
+| Beam structure: 1 channels, 2 dimensions
+| Applying user-defined cuts.
+| Starting integration for process 'pdf_builtin_1_p'
+| Integrate: iterations = 4:500:"gw", 2:500
+| Integrator: 2 chains, 2 channels, 4 dimensions
+| Integrator: Using VAMP channel equivalences
+| Integrator: 500 initial calls, 20 bins, stratified = T
+| Integrator: VAMP
+|=============================================================================|
+| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
+|=============================================================================|
+ 1 500 1.4788065E+00 5.50E-01 37.19 8.32* 1.07
+ 2 499 1.2735090E+00 1.38E-01 10.83 2.42* 4.06
+ 3 498 1.1920139E+00 6.35E-02 5.33 1.19* 14.36
+ 4 497 1.1803539E+00 5.26E-02 4.46 0.99* 17.85
+|-----------------------------------------------------------------------------|
+ 4 1994 1.1935506E+00 3.88E-02 3.25 1.45 17.85 0.22 4
+|-----------------------------------------------------------------------------|
+ 5 497 1.0590799E+00 5.91E-02 5.58 1.24 9.23
+ 6 497 1.2928696E+00 1.07E-01 8.27 1.84 4.12
+|-----------------------------------------------------------------------------|
+ 6 994 1.1138863E+00 5.17E-02 4.65 1.46 4.12 3.66 2
+|=============================================================================|
+$pdf_builtin_set = "CJ12_min"
+| RNG: Initializing TAO random-number generator
+| RNG: Setting seed for random-number generator to 8
+| Initializing integration for process pdf_builtin_1_p:
+| ------------------------------------------------------------------------
+| Process [scattering]: 'pdf_builtin_1_p'
+| Library name = 'pdf_builtin_1_lib'
+| Process index = 1
+| Process components:
+| 1: 'pdf_builtin_1_p_i1': u:d:s, ubar:dbar:sbar => mu-, mu+ [omega]
+| ------------------------------------------------------------------------
+| Beam structure: p, p => pdf_builtin
+| Beam data (collision):
+| p (mass = 0.0000000E+00 GeV)
+| p (mass = 0.0000000E+00 GeV)
+| sqrts = 1.000000000000E+03 GeV
+| Initialized builtin PDF CJ12_min
+| Phase space: generating configuration ...
+| Phase space: ... success.
+| Phase space: writing configuration file 'pdf_builtin_1_p_i1.phs'
+| Phase space: 2 channels, 2 dimensions
+| Phase space: found 2 channels, collected in 2 groves.
+| Phase space: Using 2 equivalences between channels.
+| Phase space: wood
+| Beam structure: pdf_builtin, none => none, pdf_builtin
+| Beam structure: 1 channels, 2 dimensions
+| Applying user-defined cuts.
+| Starting integration for process 'pdf_builtin_1_p'
+| Integrate: iterations = 4:500:"gw", 2:500
+| Integrator: 2 chains, 2 channels, 4 dimensions
+| Integrator: Using VAMP channel equivalences
+| Integrator: 500 initial calls, 20 bins, stratified = T
+| Integrator: VAMP
+|=============================================================================|
+| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
+|=============================================================================|
+ 1 500 9.0561608E-01 2.51E-01 27.68 6.19* 1.43
+ 2 499 1.2607338E+00 1.25E-01 9.90 2.21* 5.67
+ 3 498 1.1710580E+00 6.02E-02 5.14 1.15* 14.78
+ 4 497 1.1381306E+00 5.39E-02 4.74 1.06* 12.59
+|-----------------------------------------------------------------------------|
+ 4 1994 1.1570742E+00 3.78E-02 3.27 1.46 12.59 0.62 4
+|-----------------------------------------------------------------------------|
+ 5 497 1.2424073E+00 6.24E-02 5.02 1.12 11.97
+ 6 497 1.1813185E+00 5.64E-02 4.78 1.06* 10.83
+|-----------------------------------------------------------------------------|
+ 6 994 1.2087820E+00 4.19E-02 3.46 1.09 10.83 0.53 2
+|=============================================================================|
+| RNG: Initializing TAO random-number generator
+| RNG: Setting seed for random-number generator to 9
| Initializing integration for process pdf_builtin_2_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_2_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 2
| Process components:
| 1: 'pdf_builtin_2_p_i1': A, A => e-, e+ [omega]
| ------------------------------------------------------------------------
| Beam structure: A, A
| Beam data (collision):
| A (mass = 0.0000000E+00 GeV)
| A (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_2_p_i1.phs'
| Phase space: 4 channels, 2 dimensions
| Phase space: found 4 channels, collected in 1 grove.
| Phase space: Using 8 equivalences between channels.
| Phase space: wood
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_2_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 1 chains, 4 channels, 2 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
- 1 500 9.4963594E+02 5.16E+01 5.43 1.21* 37.26
- 2 500 1.0083599E+03 2.17E+01 2.15 0.48* 56.50
- 3 500 9.9893714E+02 1.74E+01 1.75 0.39* 64.34
- 4 500 9.6584118E+02 1.83E+01 1.89 0.42 51.08
+ 1 500 1.0371009E+03 5.30E+01 5.11 1.14* 40.63
+ 2 500 1.0456202E+03 2.12E+01 2.03 0.45* 60.49
+ 3 500 1.0155966E+03 1.75E+01 1.72 0.39* 64.12
+ 4 500 1.0059769E+03 1.91E+01 1.90 0.43 40.10
|-----------------------------------------------------------------------------|
- 4 2000 9.8782302E+02 1.07E+01 1.08 0.48 51.08 1.10 4
+ 4 2000 1.0212371E+03 1.08E+01 1.06 0.47 40.10 0.72 4
|-----------------------------------------------------------------------------|
- 5 500 1.0292175E+03 1.63E+01 1.59 0.35* 54.40
- 6 500 1.0015150E+03 1.74E+01 1.74 0.39 52.91
+ 5 500 9.9383552E+02 1.80E+01 1.81 0.40* 39.61
+ 6 500 1.0009031E+03 2.09E+01 2.09 0.47 39.87
|-----------------------------------------------------------------------------|
- 6 1000 1.0162799E+03 1.19E+01 1.17 0.37 52.91 1.35 2
+ 6 1000 9.9683579E+02 1.36E+01 1.37 0.43 39.87 0.07 2
|=============================================================================|
$pdf_builtin_set = "mrst2004qedp"
| RNG: Initializing TAO random-number generator
-| RNG: Setting seed for random-number generator to 7
+| RNG: Setting seed for random-number generator to 10
| Initializing integration for process pdf_builtin_2_p:
| ------------------------------------------------------------------------
| Process [scattering]: 'pdf_builtin_2_p'
| Library name = 'pdf_builtin_1_lib'
| Process index = 2
| Process components:
| 1: 'pdf_builtin_2_p_i1': A, A => e-, e+ [omega]
| ------------------------------------------------------------------------
| Beam structure: p, p => pdf_builtin
| Beam data (collision):
| p (mass = 0.0000000E+00 GeV)
| p (mass = 0.0000000E+00 GeV)
| sqrts = 1.000000000000E+03 GeV
| Initialized builtin PDF MRST2004QEDp
| Phase space: generating configuration ...
| Phase space: ... success.
| Phase space: writing configuration file 'pdf_builtin_2_p_i1.phs'
| Phase space: 4 channels, 2 dimensions
| Phase space: found 4 channels, collected in 1 grove.
| Phase space: Using 8 equivalences between channels.
| Phase space: wood
| Beam structure: pdf_builtin, none => none, pdf_builtin
| Beam structure: 1 channels, 2 dimensions
| Applying user-defined cuts.
| Starting integration for process 'pdf_builtin_2_p'
| Integrate: iterations = 4:500:"gw", 2:500
| Integrator: 1 chains, 4 channels, 4 dimensions
| Integrator: Using VAMP channel equivalences
| Integrator: 500 initial calls, 20 bins, stratified = T
| Integrator: VAMP
|=============================================================================|
| It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] |
|=============================================================================|
- 1 500 7.3619633E-02 3.44E-02 46.68 10.44* 1.00
- 2 500 9.2741916E-02 1.15E-02 12.44 2.78* 4.65
- 3 500 1.5306946E-01 1.35E-02 8.80 1.97* 8.30
- 4 500 1.3542554E-01 7.38E-03 5.45 1.22* 14.94
+ 1 500 2.5151015E-01 1.05E-01 41.56 9.29* 1.10
+ 2 500 1.3019355E-01 1.33E-02 10.20 2.28* 6.31
+ 3 500 1.4584333E-01 1.30E-02 8.94 2.00* 6.30
+ 4 500 1.4904614E-01 9.95E-03 6.67 1.49* 11.22
|-----------------------------------------------------------------------------|
- 4 2000 1.2686716E-01 5.57E-03 4.39 1.96 14.94 5.42 4
+ 4 2000 1.4369827E-01 6.78E-03 4.72 2.11 11.22 0.80 4
|-----------------------------------------------------------------------------|
- 5 500 1.2749759E-01 7.70E-03 6.04 1.35 9.94
- 6 500 1.3511872E-01 8.79E-03 6.51 1.45 7.88
+ 5 500 1.5107240E-01 1.12E-02 7.41 1.66 8.55
+ 6 500 1.3786013E-01 9.41E-03 6.82 1.53* 7.78
|-----------------------------------------------------------------------------|
- 6 1000 1.3080485E-01 5.79E-03 4.43 1.40 7.88 0.43 2
+ 6 1000 1.4332518E-01 7.20E-03 5.03 1.59 7.78 0.82 2
|=============================================================================|
| WHIZARD run finished.
|=============================================================================|
Index: trunk/src/pdf_builtin/pdf_builtin.f90
===================================================================
--- trunk/src/pdf_builtin/pdf_builtin.f90 (revision 5839)
+++ trunk/src/pdf_builtin/pdf_builtin.f90 (revision 5840)
@@ -1,828 +1,868 @@
!$Id$
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Copyright (C) 1999-2014 by
! Wolfgang Kilian <kilian@physik.uni-siegen.de>
! Thorsten Ohl <ohl@physik.uni-wuerzburg.de>
! Juergen Reuter <juergen.reuter@desy.de>
! Christian Speckner <cnspeckn@googlemail.com>
!
! WHIZARD is free software; you can redistribute it and/or modify it
! under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2, or (at your option)
! any later version.
!
! WHIZARD is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Wrap a common interface around the different PDF sets.
module pdf_builtin
use kinds, only: default, double
use constants, only: PI
use diagnostics
use iso_varying_string, string_t => varying_string
use file_utils
use mrst2004qed
use cteq6pdf
use mstw2008
use ct10pdf
use CJ12
implicit none
save
private
! The available sets
- integer, parameter :: nsets = 10
+ integer, parameter :: nsets = 13
integer, parameter, public :: &
CTEQ6M = 1, CTEQ6D = 2, CTEQ6L = 3, CTEQ6L1 = 4, &
MRST2004QEDp = 5, MRST2004QEDn = 6, MSTW2008LO = 7, MSTW2008NLO = 8, &
MSTW2008NNLO = 9, CT10 = 10, CJ12_max = 11, CJ12_mid = 12, &
CJ12_min = 13
! Limits
real(kind=default), parameter :: &
cteq6_q_min = 1.3, cteq6_q_max = 10.E3, &
cteq6_x_min = 1.E-6, cteq6_x_max = 1., &
mrst2004qed_q_min = sqrt (1.26_default), mrst2004qed_q_max = sqrt (0.99E7_default), &
mrst2004qed_x_min = 1.01E-5, mrst2004qed_x_max = 1., &
mstw2008_q_min = 1.01, mstw2008_q_max = sqrt (0.99E9_default), &
mstw2008_x_min = 1.01E-6, mstw2008_x_max = 1., &
ct10_q_min = 1.3, ct10_q_max = 1.E5, &
ct10_x_min = 1.E-8, ct10_x_max = 1., &
cj12_q_min = 1.69, cj12_q_max = 3.E5, &
cj12_x_min = 2.E-5, cj12_x_max = 0.9
! Lambda_QCD and quark masses
real(kind=double), dimension(6), parameter, public :: alam_ct10 = &
- (/ 0.37219423568859566_double, 0.37219423568859566_double, &
+ [ 0.37219423568859566_double, 0.37219423568859566_double, &
0.37219423568859566_double, 0.37219423568859566_double, &
- 0.32560730624033124_double, 0.22600000000000001_double /)
+ 0.32560730624033124_double, 0.22600000000000001_double ]
real(kind=double), dimension(6), parameter, public :: alam_cteq6m = &
- (/ 0.37248648555333957_double, 0.37248648555333957_double, &
+ [ 0.37248648555333957_double, 0.37248648555333957_double, &
0.37248648555333957_double, 0.37248648555333957_double, &
- 0.32590307925177625_double, 0.22623045868581182_double /)
+ 0.32590307925177625_double, 0.22623045868581182_double ]
real(kind=double), dimension(6), parameter, public :: alam_cteq6l = &
- (/ 0.37218603304249098_double, 0.37218603304249098_double, &
+ [ 0.37218603304249098_double, 0.37218603304249098_double, &
0.37218603304249098_double, 0.37218603304249098_double, &
- 0.32559900534407232_double, 0.22599353258372407_double /)
+ 0.32559900534407232_double, 0.22599353258372407_double ]
real(kind=double), dimension(6), parameter, public :: alam_cteq6ll = &
- (/ 0.24560434095679567_double, 0.24560434095679567_double, &
+ [ 0.24560434095679567_double, 0.24560434095679567_double, &
0.24560434095679567_double, 0.24560434095679567_double, &
- 0.21495099184302788_double, 0.16499885346541945_double /)
+ 0.21495099184302788_double, 0.16499885346541945_double ]
real(kind=double), dimension(0:5), parameter :: amhat = &
- (/ 0.0000000000000000_double, 0.0000000000000000_double, &
+ [ 0.0000000000000000_double, 0.0000000000000000_double, &
0.0000000000000000_double, 0.0000000000000000_double, &
- 1.3000000000000000_double, 4.500000000000000_double /)
+ 1.3000000000000000_double, 4.500000000000000_double ]
!!! Global parameter for minimal value of evolution
real(kind=double), parameter, private :: amn = 0.375_double, &
tol = .0005_double
integer, parameter, private :: nhq = 2
real(kind=double), parameter :: &
ca = 3.D0, cf = 4./3.D0, tr = 0.5D0
real(kind=double), parameter :: &
b00 = 11./3.d0 * ca, b01 = -4./3.d0 * tr, b10 = 34./3.d0 * ca**2, &
b11 = -20./3.d0 * ca*tr - 4.* cf*tr
! Init flags
integer :: cteq6_initialized = -1
integer :: mstw2008_initialized = -1
logical :: ct10_initialized = .false.
- logical :: cj12_initialized = .false.
+ integer :: cj12_initialized = -1
logical :: &
mrst2004qedp_initialized = .false., &
mrst2004qedn_initialized = .false.
type(string_t) :: mrst2004qedp_prefix, mrst2004qedn_prefix, &
mstw2008_prefix
! 1 is for proton, 2 for pion, 3 for photon
type :: pdf_builtin_status_t
private
logical, dimension(3) :: initialized = .false.
end type pdf_builtin_status_t
! Public stuff
public :: pdf_init
public :: pdf_get_name
public :: pdf_evolve
public :: pdf_evolve_LHAPDF
public :: pdf_provides_photon
public :: pdf_get_id
public :: pdf_alphas
public :: pdf_builtin_status_reset
public :: pdf_builtin_status_is_initialized
public :: pdf_builtin_status_set_initialized
public :: pdf_builtin_status_t
contains
! Get PDF name
function pdf_get_name (pdftype) result (name)
integer, intent(in) :: pdftype
type(string_t) :: name
select case (pdftype)
case (CTEQ6M)
name = var_str ("CTEQ6M")
case (CTEQ6D)
name = var_str ("CTEQ6D")
case (CTEQ6L)
name = var_str ("CTEQ6L")
case (CTEQ6L1)
name = var_str ("CTEQ6L1")
case (MRST2004QEDp)
name = var_str ("MRST2004QEDp")
case (MRST2004QEDn)
name = var_str ("MRST2004QEDn")
case (MSTW2008LO)
name = var_str ("MSTW2008LO")
case (MSTW2008NLO)
name = var_str ("MSTW2008NLO")
case (MSTW2008NNLO)
name = var_str ("MSTW2008NNLO")
case (CT10)
name = var_str ("CT10")
case (CJ12_max)
name = var_str ("CJ12_max")
case (CJ12_mid)
name = var_str ("CJ12_mid")
case (CJ12_min)
name = var_str ("CJ12_min")
case default
call msg_fatal ("pdf_builtin: internal: invalid PDF set!")
end select
end function pdf_get_name
! Get the ID of a PDF set
function pdf_get_id (name) result (id)
type(string_t), intent(in) :: name
integer :: id
do id = 1, nsets
if (upper_case (pdf_get_name (id)) == upper_case (name)) return
end do
id = -1
end function pdf_get_id
! Query whether a PDF supplies a photon distribution
function pdf_provides_photon (pdftype) result (flag)
integer, intent(in) :: pdftype
logical :: flag
select case (pdftype)
case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1)
flag = .false.
case (MRST2004QEDp, MRST2004QEDn)
flag = .true.
case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO)
flag = .false.
case (CT10)
flag = .false.
case (CJ12_max, CJ12_mid, CJ12_min)
flag = .false.
case default
call msg_fatal ("pdf_builtin: internal: invalid PDF set!")
end select
end function pdf_provides_photon
! Initialize a PDF
subroutine pdf_init (pdf_status, pdftype, prefix, verbose)
integer, intent(in) :: pdftype
type(pdf_builtin_status_t), intent(inout) :: pdf_status
type(string_t), intent(in), optional :: prefix
type(string_t) :: mprefix
logical, intent(in), optional :: verbose
logical :: mverbose
if (any (pdf_status%initialized)) return
if (present (prefix)) then
mprefix = prefix
else
mprefix = ""
end if
if (present (verbose)) then
mverbose = verbose
else
mverbose = .true.
end if
select case (pdftype)
case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1)
if (cteq6_initialized == pdftype) return
call setctq6 (pdftype, char (mprefix))
cteq6_initialized = pdftype
case (MRST2004QEDp)
if (mrst2004qedp_initialized) return
mrst2004qedp_initialized = .true.
mrst2004qedp_prefix = mprefix
case (MRST2004QEDn)
if (mrst2004qedn_initialized) return
mrst2004qedn_initialized = .true.
mrst2004qedn_prefix = mprefix
case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO)
if (mstw2008_initialized == pdftype) return
mstw2008_initialized = pdftype
mstw2008_prefix = mprefix
case (CT10)
if (ct10_initialized) return
call setct10 (char (mprefix), 100)
ct10_initialized = .true.
+ case (CJ12_max)
+ if (cj12_initialized == pdftype) return
+ call setCJ (char (mprefix), 300)
+ cj12_initialized = pdftype
+ case (CJ12_mid)
+ if (cj12_initialized == pdftype) return
+ call setCJ (char (mprefix), 200)
+ cj12_initialized = pdftype
+ case (CJ12_min)
+ if (cj12_initialized == pdftype) return
+ call setCJ (char (mprefix), 100)
+ cj12_initialized = pdftype
case default
call msg_fatal ("pdf_builtin: internal: invalid PDF set!")
end select
if (mverbose) call msg_message ("Initialized builtin PDF " // &
char (pdf_get_name (pdftype)))
!!! Up to now only proton
call pdf_builtin_status_set_initialized (pdf_status, 1)
end subroutine pdf_init
! Evolve PDF
subroutine pdf_evolve (pdftype, x, q, f, fphoton)
integer, intent(in) :: pdftype
real(kind=default), intent(in) :: x, q
real(kind=default), intent(out), optional :: f(-6:6), fphoton
real(kind=double) :: mx, mq
real(kind=double) :: upv, dnv, ups, dns, str, chm, bot, glu, phot, &
sbar, bbar, cbar
type(string_t) :: setname, prefix
select case (pdftype)
case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1)
if (cteq6_initialized < 0) &
call msg_fatal ("pdf_builtin: internal: PDF set " // &
char (pdf_get_name (pdftype)) // " requested without initialization!")
if (cteq6_initialized /= pdftype) &
call msg_fatal ( &
"PDF sets " // char (pdf_get_name (pdftype)) // " and " // &
char (pdf_get_name (cteq6_initialized)) // &
" cannot be used simultaneously")
mx = max (min (x, cteq6_x_max), cteq6_x_min)
mq = max (min (q, cteq6_q_max), cteq6_q_min)
- if (present (f)) f = (/ 0._double, &
+ if (present (f)) f = [ 0._double, &
ctq6pdf (-5, mx, mq), ctq6pdf (-4, mx, mq), &
ctq6pdf (-3, mx, mq), ctq6pdf (-1, mx, mq), &
ctq6pdf (-2, mx, mq), ctq6pdf ( 0, mx, mq), &
ctq6pdf ( 2, mx, mq), ctq6pdf ( 1, mx, mq), &
ctq6pdf ( 3, mx, mq), ctq6pdf ( 4, mx, mq), &
- ctq6pdf ( 5, mx, mq), 0._double /)
+ ctq6pdf ( 5, mx, mq), 0._double ]
if (present (fphoton)) call msg_fatal ("photon pdf requested for " // &
char (pdf_get_name (pdftype)) // " which does not provide it!")
case (MRST2004QEDp)
if (.not. mrst2004qedp_initialized) call msg_fatal ( &
"pdf_builtin: internal: PDF set MRST2004QEDp requested without " // &
"initialization!")
mx = max (min (x, mrst2004qed_x_max), mrst2004qed_x_min)
mq = max (min (q, mrst2004qed_q_max), mrst2004qed_q_min)
call mrstqed (mx, mq, 1, upv, dnv, ups, dns, str, chm, bot, glu, phot, &
char (mrst2004qedp_prefix))
if (present (f)) f = &
- (/ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, &
- str, chm, bot, 0._double /) / mx
+ [ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, &
+ str, chm, bot, 0._double ] / mx
if (present (fphoton)) fphoton = phot / mx
case (MRST2004QEDn)
if (.not. mrst2004qedn_initialized) call msg_fatal ( &
"pdf_builtin: internal: PDF set MRST2004QEDn requested without " // &
"initialization!")
mx = max (min (x, mrst2004qed_x_max), mrst2004qed_x_min)
mq = max (min (q, mrst2004qed_q_max), mrst2004qed_q_min)
call mrstqed (mx, mq, 2, upv, dnv, ups, dns, str, chm, bot, glu, phot, &
char (mrst2004qedn_prefix))
if (present (f)) f = &
- (/ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, &
- str, chm, bot, 0._double /) / mx
+ [ 0._double, bot, chm, str, ups, dns, glu, dns + dnv, ups + upv, &
+ str, chm, bot, 0._double ] / mx
if (present (fphoton)) fphoton = phot / mx
case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO)
if (mstw2008_initialized < 0) &
call msg_fatal ("pdf_builtin: internal: PDF set " // &
char (pdf_get_name (pdftype)) // " requested without initialization!")
if (mstw2008_initialized /= pdftype) &
call msg_fatal ( &
"PDF sets " // char (pdf_get_name (pdftype)) // " and " // &
char (pdf_get_name (mstw2008_initialized)) // &
" cannot be used simultaneously")
select case (pdftype)
case (MSTW2008LO)
setname = var_str ("mstw2008lo")
case (MSTW2008NLO)
setname = var_str ("mstw2008nlo")
case (MSTW2008NNLO)
setname = var_str ("mstw2008nnlo")
end select
mx = max (min (x, mstw2008_x_max), mstw2008_x_min)
mq = max (min (q, mstw2008_q_max), mstw2008_q_min)
call getmstw2008 (char (mstw2008_prefix), char (setname), 0, mx, mq, &
upv, dnv, ups, dns, str, sbar, chm, cbar, bot, bbar, glu, phot)
if (present (f)) f = &
- (/ 0._double, bbar, cbar, sbar, ups, dns, glu, dns + dnv, ups + upv, &
- str, chm, bot, 0._double /) / mx
+ [ 0._double, bbar, cbar, sbar, ups, dns, glu, dns + dnv, ups + upv, &
+ str, chm, bot, 0._double ] / mx
if (present (fphoton)) call msg_fatal ("photon pdf requested for " // &
char (pdf_get_name (pdftype)) // " which does not provide it!")
case (CT10)
if (.not. ct10_initialized) call msg_fatal ( &
"pdf_builtin: internal: PDF set CT10 requested without " // &
"initialization!")
mx = max (min (x, ct10_x_max), ct10_x_min)
mq = max (min (q, ct10_q_max), ct10_q_min)
- if (present (f)) f = (/ 0._double, &
+ if (present (f)) f = [ 0._double, &
getct10pdf (-5, mx, mq), getct10pdf (-4, mx, mq), &
getct10pdf (-3, mx, mq), getct10pdf (-1, mx, mq), &
getct10pdf (-2, mx, mq), getct10pdf ( 0, mx, mq), &
getct10pdf ( 2, mx, mq), getct10pdf ( 1, mx, mq), &
getct10pdf ( 3, mx, mq), getct10pdf ( 4, mx, mq), &
- getct10pdf ( 5, mx, mq), 0._double /)
+ getct10pdf ( 5, mx, mq), 0._double ]
if (present (fphoton)) call msg_fatal ("photon pdf requested for " // &
"CT10 which does not provide it!")
+ case (CJ12_max, CJ12_mid, CJ12_min)
+ if (cj12_initialized < 0) &
+ call msg_fatal ("pdf_builtin: internal: PDF set " // &
+ char (pdf_get_name (pdftype)) // " requested without initialization!")
+ if (cj12_initialized /= pdftype) &
+ call msg_fatal ( &
+ "PDF sets " // char (pdf_get_name (pdftype)) // " and " // &
+ char (pdf_get_name (cj12_initialized)) // &
+ " cannot be used simultaneously")
+ mx = max (min (x, cj12_x_max), cj12_x_min)
+ mq = max (min (q, cj12_q_max), cj12_q_min)
+ if (present (f)) f = [ 0._double, &
+ CJpdf (-5, mx, mq), CJpdf (-4, mx, mq), CJpdf (-3, mx, mq), &
+ CJpdf (-2, mx, mq), CJpdf (-1, mx, mq), CJpdf ( 0, mx, mq), &
+ CJpdf ( 1, mx, mq), CJpdf ( 2, mx, mq), CJpdf ( 3, mx, mq), &
+ CJpdf ( 4, mx, mq), CJpdf ( 5, mx, mq), 0._double ]
+ if (present (fphoton)) call msg_fatal ("photon pdf requested for " // &
+ "CJ12 which does not provide it!")
case default
call msg_fatal ("pdf_builtin: internal: invalid PDF set!")
end select
end subroutine pdf_evolve
! included for compatibility with LHAPDF
! use a double precision array for the pdfs instead of a
! default precision
subroutine pdf_evolve_LHAPDF (set, x, q, ff)
integer, intent(in) :: set
double precision, intent(in) :: x, q
real(kind=default) :: dx, dq
double precision, dimension(-6:6), intent(out) :: ff
real(kind=default) :: f(-6:6)
dx = x
dq = q
call pdf_evolve(set, dx, dq, f)
ff = f
end subroutine pdf_evolve_LHAPDF
! PDF-specific running alphas
function pdf_alphas (pdftype, q) result (alphas)
real(kind=double) :: as
real(kind=default), intent(in) :: q
real(kind=double) :: qdummy
real(kind=default) :: alphas
integer, intent(in) :: pdftype
integer, parameter :: nset_dummy = 1
qdummy = q
select case (pdftype)
case (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1)
if (cteq6_initialized < 0) &
call msg_fatal ("pdf_builtin: internal: PDF set " // &
char (pdf_get_name (pdftype)) // " requested without initialization!")
if (cteq6_initialized /= pdftype) &
call msg_fatal ( &
"PDF sets " // char (pdf_get_name (pdftype)) // " and " // &
char (pdf_get_name (cteq6_initialized)) // &
" cannot be used simultaneously")
select case (pdftype)
case (CTEQ6M, CTEQ6D)
as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m)
case (CTEQ6L)
as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6l)
case (CTEQ6L1)
as = PI*pdf_builtin_alpi (qdummy,1,5,alam_cteq6ll)
case default
call msg_fatal ("Unrecognized PDF setup in evolution of alphas.")
end select
case (CT10)
if (.not. ct10_initialized) call msg_fatal ( &
"pdf_builtin: internal: PDF set CT10 requested without " // &
"initialization!")
as = PI*pdf_builtin_alpi (qdummy,2,5,alam_ct10)
case (MRST2004QEDp)
if (.not. mrst2004qedp_initialized) call msg_fatal ( &
"pdf_builtin: internal: PDF set MRST2004QEDp requested without " // &
"initialization!")
call mrst_alphas (as,qdummy)
case (MRST2004QEDn)
if (.not. mrst2004qedn_initialized) call msg_fatal ( &
"pdf_builtin: internal: PDF set MRST2004QEDn requested without " // &
"initialization!")
call mrst_alphas (as,qdummy)
case (MSTW2008LO, MSTW2008NLO, MSTW2008NNLO)
if (mstw2008_initialized < 0) &
call msg_fatal ("pdf_builtin: internal: PDF set " // &
char (pdf_get_name (pdftype)) // " requested without initialization!")
if (mstw2008_initialized /= pdftype) &
call msg_fatal ( &
"PDF sets " // char (pdf_get_name (pdftype)) // " and " // &
char (pdf_get_name (mstw2008_initialized)) // &
" cannot be used simultaneously")
select case (pdftype)
case (MSTW2008lo)
call pdf_builtin_alfa (as,qdummy,0)
case (MSTW2008nlo)
call pdf_builtin_alfa (as,qdummy,1)
case (MSTW2008nnlo)
call pdf_builtin_alfa (as,qdummy,2)
end select
+ case (CJ12_max, CJ12_mid, CJ12_min)
+ if (cj12_initialized < 0) &
+ call msg_fatal ("pdf_builtin: internal: PDF set " // &
+ char (pdf_get_name (pdftype)) // " requested without initialization!")
+ if (cj12_initialized /= pdftype) &
+ call msg_fatal ( &
+ "PDF sets " // char (pdf_get_name (pdftype)) // " and " // &
+ char (pdf_get_name (cteq6_initialized)) // &
+ " cannot be used simultaneously")
+ as = PI*pdf_builtin_alpi (qdummy,2,5,alam_cteq6m)
case default
call msg_fatal ("pdf_builtin: internal: invalid PDF set!")
end select
alphas = as
end function pdf_alphas
subroutine pdf_builtin_status_reset (pdf_builtin_status)
type(pdf_builtin_status_t), intent(inout) :: pdf_builtin_status
pdf_builtin_status%initialized = .false.
end subroutine pdf_builtin_status_reset
function pdf_builtin_status_is_initialized (pdf_builtin_status, set) result (flag)
logical :: flag
type(pdf_builtin_status_t), intent(in) :: pdf_builtin_status
integer, intent(in), optional :: set
if (present (set)) then
select case (set)
case (1:3); flag = pdf_builtin_status%initialized(set)
case default; flag = .false.
end select
else
flag = any (pdf_builtin_status%initialized)
end if
end function pdf_builtin_status_is_initialized
subroutine pdf_builtin_status_set_initialized (pdf_builtin_status, set)
type(pdf_builtin_status_t), intent(inout) :: pdf_builtin_status
integer, intent(in) :: set
pdf_builtin_status%initialized(set) = .true.
end subroutine pdf_builtin_status_set_initialized
function pdf_builtin_alpi (AMU,NORDER,NF,alam)
integer, intent(in) :: nf, norder
real(kind=double), intent(in) :: amu
real(kind=double), dimension(6), intent(in) :: alam
real(kind=double) :: alm
real(kind=double) :: pdf_builtin_alpi
integer :: irt, n_eff
n_eff = num_flavor (AMU,NF)
alm = alam (n_eff+1)
pdf_builtin_alpi = pdf_builtin_alpqcd (norder, n_eff, amu/alm, irt)
select case (irt)
case (1)
call msg_warning ("AMU < ALAM in pdf_builtin_alpi")
case (2)
call msg_warning ("pdf_builtin_alpi > 3; Be aware!")
end select
end function pdf_builtin_alpi
double precision function pdf_builtin_ALPQCD (IRDR, NF, RML, IRT)
real(kind=double) :: D0 = 0.D0, D1 = 1.D0, BIG = 1.0D15
real(kind=double) :: b0, b1
integer, intent(in) :: irdr, nf
integer :: irt
real(kind=double), intent(in) :: rml
real(kind=double) :: al, aln
real(kind=double) :: rm2
IRT = 0
if (IRDR .LT. 1 .OR. IRDR .GT. 2) THEN
print *, &
& 'Order out of range in pdf_builtin_ALPQCD: IRDR = ', IRDR
STOP
end if
b0 = (11.d0*CA - 2.* NF) / 3.d0
b1 = (34.d0*CA**2 - 10.d0*CA*NF - 6.d0*CF*NF) / 3.d0
rm2 = rml**2
if (RM2 .LE. 1.) THEN
IRT = 1
pdf_builtin_ALPQCD = 99.
RETURN
end if
aln = log (RM2)
al = 4.d0/ B0 / aln
IF (IRDR .GE. 2) AL = AL * (1.d0-B1*log(ALN) / ALN / B0**2)
if (AL .GE. 3.) THEN
IRT = 2
end if
pdf_builtin_ALPQCD = AL
end function pdf_builtin_ALPQCD
function num_flavor (amu,nf)
real(kind=double), intent(in) :: amu
integer, intent(in) :: nf
integer :: num_flavor, i
num_flavor = nf - nhq
IF ((num_flavor .EQ. nf) .OR. (amu .LE. amn)) GOTO 20
do 10 I = nf - NHQ + 1, nf
if (amu .GE. amhat(I)) THEN
num_flavor = I
else
goto 20
end if
10 continue
20 return
end function num_flavor
subroutine pdf_builtin_alfa (alfas,Qalfa,norder)
integer, intent(in) :: norder
real(kind=double), intent(in) :: Qalfa
real(kind=double), intent(out) :: alfas
real(kind=double) :: alphasq0
real(kind=double) :: mCharmSave, mBottomSave, mTopSave
mTopSave = 10000000000.0_double
mBottomSave = 4.75_double
mCharmSave = 1.4_double
select case (norder)
case (0)
alphasq0 = 0.68183_double
case (1)
alphasq0 = 0.491278_double
case (2)
alphasq0 = 0.45077_double
case default
call msg_fatal ("Wrong MSTW order!")
end select
if (Qalfa.GT.mTopSave) call msg_fatal &
("No MSTW alphas for such high scales.")
alfas = pdf_builtin_mstw_alphas (Qalfa,norder,mTopSave**2)
end subroutine pdf_builtin_alfa
function pdf_builtin_mstw_alphas (mur, norder, m2t) result (alfas)
real(kind=double), intent(in) :: mur
real(kind=double) :: alfas
integer, intent(in) :: norder
integer :: nf
real(kind=double), intent(in) :: m2t
real(kind=double) :: as0, asc, asb
real(kind=double) :: r2, m2, asi, asf, r20, r2b, r2c, r2t
real(kind=double), parameter :: m20 = 1.0_double, m2c = 1.96_double, &
m2b = 22.5625_double
real(kind=double) :: logfr
integer :: nff
real(kind=double) :: ast
nff = 4
select case (norder)
case (0)
as0 = 5.4258307424173556E-002_double
asc = 4.0838232991033577E-002_double
asb = 2.2297506503539639E-002_double
case (1)
as0 = 3.9094820221093209E-002_double
asc = 3.0202751103489092E-002_double
asb = 1.7751676069301441E-002_double
case (2)
as0 = 3.5871136848766867E-002_double
asc = 2.8092349518054758E-002_double
asb = 1.6936586710596203E-002_double
case default
call msg_fatal ("Wrong evolution order in MSTW alphas evolution.")
end select
r2 = mur**2
m2 = r2 * exp(+logfr)
if (m2 .gt. m2t) then
nf = 6
r2t = m2t * r2/m2
asi = ast
asf = mstw_alphas (r2,r2t,ast,nf,norder)
else
if (m2 .gt. m2b) then
nf = 5
r2b = m2b * r2/m2
asi = asb
asf = mstw_alphas (r2,r2b,asb,nf,norder)
else
if (m2 .gt. m2c) then
nf = 4
r2c = m2c * r2/m2
asi = asc
asf = mstw_alphas (r2,r2c,asc,nf,norder)
else
nf = 3
r20 = m20 * r2/m2
asi = as0
asf = mstw_alphas (r2,r20,as0,nf,norder)
end if
end if
end if
alfas = 4.D0*PI*ASF
end function pdf_builtin_mstw_alphas
function mstw_alphas (r2, r20, as0, nf, naord) result (as)
real(kind=double) :: as
real(kind=double), intent(in) :: r2, r20, as0
integer, intent(in) :: nf, naord
integer, parameter :: nastps = 20
integer :: k1
real(kind=double), parameter :: sxth = 1./6.
real(kind=double) :: lrrat, dlr, xk0, xk1, xk2, xk3
as = as0
lrrat = log (r2/r20)
dlr = lrrat / nastps
! ..Solution of the evolution equation depending on NAORD
! (fourth-order Runge-Kutta beyond the leading order)
select case (naord)
case (0)
as = as0 / (1.+ beta0(nf) * as0 * lrrat)
case (1)
do k1 = 1, nastps
xk0 = dlr * fbeta1 (as, nf)
xk1 = dlr * fbeta1 (as + 0.5 * xk0, nf)
xk2 = dlr * fbeta1 (as + 0.5 * xk1, nf)
xk3 = dlr * fbeta1 (as + xk2, nf)
as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3)
end do
case (2)
do k1 = 1, nastps
xk0 = dlr * fbeta2 (as, nf)
xk1 = dlr * fbeta2 (as + 0.5 * xk0, nf)
xk2 = dlr * fbeta2 (as + 0.5 * xk1, nf)
xk3 = dlr * fbeta2 (as + xk2, nf)
as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3)
end do
case (3)
do k1 = 1, nastps
xk0 = dlr * fbeta3 (as, nf)
xk1 = dlr * fbeta3 (as + 0.5 * xk0, nf)
xk2 = dlr * fbeta3 (as + 0.5 * xk1, nf)
xk3 = dlr * fbeta3 (as + xk2, nf)
as = as + sxth * (xk0 + 2.* xk1 + 2.* xk2 + xk3)
end do
end select
end function mstw_alphas
double precision function beta0 (nf)
integer, intent(in) :: nf
beta0 = b00 + b01 * nf
end function beta0
double precision function beta1 (nf)
integer, intent(in) :: nf
beta1 = b10 + b11 * nf
end function beta1
double precision function beta2 (nf)
integer, intent(in) :: nf
beta2 = 1428.50 - 279.611 * nf + 6.01852 * nf**2
end function beta2
double precision function beta3 (nf)
integer, intent(in) :: nf
beta3 = 29243.0 - 6946.30 * nf + 405.089 * nf**2 + &
1.49931 * nf**3
end function beta3
! ..The beta functions FBETAn at N^nLO for n = 1, 2, and 3
double precision function fbeta1 (a, nf)
real(kind=double), intent(in) :: a
integer, intent(in) :: nf
fbeta1 = - a**2 * ( beta0(nf) + a * beta1(nf) )
end function fbeta1
double precision function fbeta2 (a, nf)
real(kind=double), intent(in) :: a
integer, intent(in) :: nf
fbeta2 = - a**2 * ( beta0(nf) + a * ( beta1(nf) &
+ a * beta2(nf) ) )
end function fbeta2
double precision function fbeta3 (a, nf)
real(kind=double), intent(in) :: a
integer, intent(in) :: nf
fbeta3 = - a**2 * ( beta0(nf) + a * ( beta1(nf) &
+ a * ( beta2(nf) + a * beta3(nf)) ) )
end function fbeta3
subroutine mrst_alphas (alpha, q)
real(kind=double), intent(in) :: q
real(kind=double), intent(out) :: alpha
integer :: idir
integer, parameter :: nflav = 4
real(kind=double) :: alphas, qz2, qz, alambda, astep, &
tol2
alphas = 0.1205_double
qz2=8315._double
qz = sqrt(qz2)
alambda = 0.3000_double
astep = 0.010_double
tol2 = 0.0000001_double
idir = +1
100 continue
alambda = alambda + idir*astep
call mrst_lambda (nflav,alpha,qz,alambda)
if (idir*(alphas-alpha).gt.0.0) then
goto 200
else
astep = 0.5*astep
idir = -1*idir
goto 200
end if
200 continue
if(abs(alpha-alphas).gt.tol2) goto 100
! alambda found -- save it !!!!
! next call mrst_lambda to get alphas at q with the correct alambda
call mrst_lambda (nflav,alpha,q,alambda)
RETURN
end subroutine mrst_alphas
subroutine mrst_lambda (nflav,alpha,Q,alambda)
integer, intent(in) :: nflav
integer :: ith
real(kind=double), intent(in) :: q, alambda
real(kind=double), intent(out) :: alpha
real(kind=double) :: al, al2, alfinv, alfqc3, alfqc4, alfqc5, &
alfqs3, alfqs5
real(kind=double) :: x1, x2, t, tt
real(kind=double) :: qsdt, qsdtt, qsct, qsctt, qs, q2
real(kind=double) :: b0, b1, del, f, fp, as, as2
real(kind=double) :: flav
! The value of Lambda required corresponds to nflav=4
qsdt=8.18 !! This is the value of 4m_c^2
qsct=74.0 !! This is the value of 4m_b^2
al2=alambda*alambda
q2=q*q
t = log (q2/al2)
! CHECK: explicitly initialising ALFQC{3,4,5} (by AB)
alfqc3 = 0
alfqc4 = 0
alfqc5 = 0
ith=0
tt=t
qsdtt=qsdt/4.
qsctt=qsct/4.
al = alambda
al2=al*al
flav=4.
qs=al2*exp(T)
if (qs.lt.0.5d0) then !! running stops below 0.5
qs=0.5d0
t = log (qs/al2)
tt= t
end if
if (qs.lt.qsdtt .and. nflav.gt.3) GOTO 312
11 continue
b0 = 11._double - 2./3._double * flav
x1 = 4.*PI/b0
b1=102.-38.*FLAV/3.
X2=B1/B0**2
AS2=X1/T*(1.-X2* log(T)/T)
5 AS=AS2
F=-T+X1/AS-X2*log(X1/AS+X2)
FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
AS2=AS-F/FP
DEL=ABS(F/FP/AS)
IF((DEL-tol).GT.0.) goto 5
ALPHA=AS2
51 continue
IF(ITH.EQ.0) RETURN
GOTO (13,14,15) ITH
! GOTO 5
12 ITH=1
T=log(QSCTT/AL2)
GOTO 11
13 ALFQC4=ALPHA
FLAV=5.
ITH=2
GOTO 11
14 ALFQC5=ALPHA
ITH=3
T=TT
GOTO 11
15 ALFQS5=ALPHA
ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
ALPHA=1./ALFINV
RETURN
311 CONTINUE
B0=11-2.*FLAV/3.
X1=4.*PI/B0
B1=102.-38.*FLAV/3.
X2=B1/B0**2
AS2=X1/T*(1.-X2*log(T)/T)
35 AS=AS2
F=-T+X1/AS-X2*log(X1/AS+X2)
FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
AS2=AS-F/FP
DEL=ABS(F/FP/AS)
IF((DEL-tol).GT.0.) goto 35
ALPHA=AS2
351 continue
IF(ITH.EQ.0) RETURN
GOTO (313,314,315) ITH
312 ITH=1
T=log(QSDTT/AL2)
GOTO 311
313 ALFQC4=ALPHA
FLAV=3.
ITH=2
GOTO 311
314 ALFQC3=ALPHA
ITH=3
T=TT
GOTO 311
315 ALFQS3=ALPHA
ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
ALPHA=1./ALFINV
end subroutine mrst_lambda
end module pdf_builtin
Index: trunk/src/pdf_builtin/CJpdf.f
===================================================================
--- trunk/src/pdf_builtin/CJpdf.f (revision 5839)
+++ trunk/src/pdf_builtin/CJpdf.f (revision 5840)
@@ -1,362 +1,363 @@
c This is the set of user routines which can be used to read the CJ PDF
c tables and extract values of the CJ PDFs at specified values of
c momentum fraction x and factorization scale Q.
c
c
c Patterned after the CTEQ6 PDF routines
c
c J.F. Owens June 2012 - December 2012
c
c Three sets of PDFs are available: CJ12_min, CJ12_mid, and CJ12_max
c corresponding to the minimum, midddle, and maximum nuclear corrections.
c
c The following numbering scheme is used for the variable ISET
c
c ISET PDF
c
c 100 | CJ12_min central PDF
c 101-138 | CJ12_min error PDF sets
c |
c 200 | CJ12_mid central PDF
c 201-238 | CJ12_mid error PDF sets
c |
c 300 | CJ12_max central PDF
c 301-338 | CJ12_max error PDF sets
c
c The tables cover the x range 10^-6 < x < 1. and Q range 1.3 < Q < 10^5 GeV.
c Values outside these ranges must be considered as extrapolations.
c
c The user should initalize the PDF set by first calling SetCJ(ISET)
c
c A function call to CJpdf(Iptn,x,Q) returns the number density PDF of
c flavor Iptn at momentum fraction x and factorization scale Q
c
c A subroutine call to CJpdf_all(x,Q,bb,cb,sb,db,ub,g,u,d,s,c,b) returns all
c the PDFs with one call
c
c Flavor index: Iptn = -5:5 for bb, cb, sb, db, ub, g, u, d, s, c, b,
c
c For the CJ12 PDFs b=bb, c=cb, and s=sb
c
c An example of how to use the CJ PDFs is included in the program tst_CJpdf.f
c This will produce a table of PDFs for a specified value of ISET
c
module CJ12
use file_utils
contains
Function CJpdf(Iptn, x, Q)
Implicit double precision (A-H,O-Z)
Common/CJPAR2/Nx, Nq, Nfmx
Common/QCDtable/Alam, Nfl, Iord
If(x.lt.0.d0.or.x.gt.1.d0)then
Print*,'x out of range in CJPDF:',x
Stop
Endif
If(Q.lt.Alam)then
Print*,'Q out of range in CJpdf:',Q
Stop
Endif
If(Iptn.lt.-Nfmx.or.Iptn.gt.Nfmx)then
Print*,'Iptn out of range in CJpdf:',Iptn
Print*,'Maximum number of flavors is:',Nfmx
Stop
Endif
CJpdf=Getpdf(Iptn,x,Q)
c if(CJpdf.lt.0.d0)CJPDF=0.d0
Return
End
Subroutine CJpdf_all(x,Q,bb,cb,sb,db,ud,g,u,d,s,c,b)
Implicit Double Precision (a-h,o-z)
bb=CJpdf(-5,x,Q)
cb=CJpdf(-4,x,Q)
sb=CJpdf(-3,x,Q)
db=CJpdf(-2,x,Q)
ub=CJpdf(-1,x,Q)
g=CJpdf(0,x,Q)
u=CJpdf(1,x,Q)
d=CJpdf(2,x,Q)
s=sb
c=cb
b=bb
return
end
- Subroutine SetCJ(Iset)
+ Subroutine SetCJ (prefix, Iset)
Implicit Double Precision (A-H,O-Z)
Character Flnm(3)*9, nn*3, Tablefile*40
+ character (*) prefix
Data (Flnm(I), I=1,3)
> / 'CJ12_min_', 'CJ12_mid_', 'CJ12_max_' /
Data Isetold/-1/
save
C If data file not initialized, do so.
If(Iset.ne.Isetold) then
- IU= NextUn()
- If (Iset.ge.100 .and. Iset.le.140) Then
- write(nn,'(I3)') Iset
- Tablefile=Flnm(1)//nn(2:3)//'.tbl'
- Elseif (Iset.ge.200 .and. Iset.le.240) Then
- write(nn,'(I3)') Iset
- Tablefile=Flnm(2)//nn(2:3)//'.tbl'
- Elseif (Iset.ge.300 .and. Iset.le.340) Then
- write(nn,'(I3)') Iset
- Tablefile=Flnm(3)//nn(2:3)//'.tbl'
- Else
- Print *, 'Invalid Iset number in SetCJ :', Iset
- Stop
- Endif
- print*,'Opening ',Tablefile
- Open(IU, File=Tablefile, Status='OLD', Err=100)
+ IU= NextUn()
+ If (Iset.ge.100 .and. Iset.le.140) Then
+ write(nn,'(I3)') Iset
+ Tablefile=Flnm(1)//nn(2:3)//'.tbl'
+ Elseif (Iset.ge.200 .and. Iset.le.240) Then
+ write(nn,'(I3)') Iset
+ Tablefile=Flnm(2)//nn(2:3)//'.tbl'
+ Elseif (Iset.ge.300 .and. Iset.le.340) Then
+ write(nn,'(I3)') Iset
+ Tablefile=Flnm(3)//nn(2:3)//'.tbl'
+ Else
+ Print *, 'Invalid Iset number in SetCJ :', Iset
+ Stop
+ Endif
+C print*,'Opening ',Tablefile
+ Open(IU, File=prefix // "/" // Tablefile, Status='OLD', Err=100)
21 Call ReadTbl (IU)
Close (IU)
Isetold=Iset
Endif
Return
- 100 Print *, ' Data file ', Tablefile, ' cannot be opened '
- >//'in SetCJ!!'
+ 100 Print *, ' Data file ', prefix // "/" // trim (Tablefile),
+ >' cannot be opened in SetCJ!!'
Stop
C ********************
End
Subroutine ReadTbl (Nu)
Implicit Double Precision (A-H,O-Z)
Character Line*80
PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
PARAMETER (MXPQX=(MXF+3)*MXQ*MXX)
Common
> / CJPAR1 / Al, XV(MXX), QV(MXQ), SV(MXQ), UPD(MXPQX)
> / CJPAR2 / Nx, NQ, NfMx
> / XQrange / Qini, Qmax, Xmin
> / QCDtable / Alam, Nfl, Iorder
> / Masstbl / Amass(6)
Read (Nu, '(A)') Line
Read (Nu, '(A)') Line
Read (Nu, *) Dr, Fl, Al, (Amass(I),I=1,6)
Iorder = Nint(Dr)
Nfl = Nint(Fl)
Alam = Al
Read (Nu, '(A)') Line
Read (Nu, *) NX, NQ, NfMx
Read (Nu, '(A)') Line
Read (Nu, *) QINI, QMAX, (QV(I), I =1, NQ)
Read (Nu, '(A)') Line
Read (Nu, *) XMIN, (XV(I), I =1, NX)
Do 11 Iq = 1, NQ
SV(Iq) = Log(Log (QV(Iq) /Alam)/Log(Qini/Alam))
11 Continue
c
Nblk=Nx*NQ
Npts=Nblk*(NfMx+3)
Read(Nu,'(A)') Line
Read(Nu,25,err=999) (UPD(I), I=1,Npts)
25 format(5E14.4)
Return
999 print*,i,UPD(i-1),UPD(i)
stop
End
Function NextUn()
C Returns an unallocated FORTRAN i/o unit.
Logical EX
C
Do 10 N = 10, 300
INQUIRE (UNIT=N, OPENED=EX)
If (.NOT. EX) then
NextUn = N
Return
Endif
10 Continue
Stop ' There is no available I/O unit. '
C *************************
End
C
Function GetPDF(Iptn,x,Q)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
PARAMETER (MXPQX=(MXF+3)*MXQ*MXX)
Common
> / CJPAR1 / Al, XV(MXX), QV(MXQ), SV(MXQ), UPD(MXPQX)
> / CJPAR2 / Nx, NQ, NfMx
> / XQrange / Qini, Qmax, Xmin
> / QCDtable / Alam, Nfl, Iorder
> / Masstbl / Amass(6)
c
c Linear interpolation in s
c
s=log(log(Q/Alam)/log(Qini/Alam))
do 2 j=1,NQ
if(s.lt.SV(j))then
J2=J
if(J2.eq.1)J2=2
J1=J2-1
S2=SV(j2)
S1=SV(j2-1)
goto 1
endif
2 continue
1 continue
c
c In CJ12 s=sbar, c=cbar, and b=bbar
c
if(Iptn.gt.2)then
Iptn_temp=-Iptn
else
Iptn_temp=Iptn
endif
A1=GetFV(Iptn_temp,x,J1)
A2=GetFV(Iptn_temp,x,J2)
ANS=A1*(S-S2)/(S1-S2)+A2*(S-S1)/(S2-S1)
if(ans.lt.0.d0)ans=0.d0
GetPDF=ANS
RETURN
END
c
FUNCTION GetFV(Iptn,x,J)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
PARAMETER (MXPQX=(MXF+3)*MXQ*MXX)
Dimension xx(4),fx(4)
Common
> / CJPAR1 / Al, XV(MXX), QV(MXQ), SV(MXQ), UPD(MXPQX)
> / CJPAR2 / Nx, NQ, NfMx
> / XQrange / Qini, Qmax, Xmin
> / QCDtable / Alam, Nfl, Iorder
> / Masstbl / Amass(6)
c
c Interpolate in x using the form A*x**alpha*(1-x)**beta
c which is locally valid for each PDF
c
DO 1 I=1,Nx
IF(X.LT.XV(I)) GO TO 2
1 CONTINUE
c 2 I=I-1
2 I=I-2
If(I.le.0.d0) I=2
if(i.gt.(nx-3))i=nx-3
In1=I+Nx*(J-1)+Nx*NQ*(Iptn+5)
xx(1)=xv(I)
xx(2)=xv(i+1)
xx(3)=xv(I+2)
xx(4)=xv(I+3)
fx(1)=UPD(In1)
fx(2)=UPD(In1+1)
fx(3)=UPD(In1+2)
fx(4)=UPD(In1+3)
c call polint(xx,fx,4,x,ans,dy)
call polint4(xx,fx,x,ans)
getfv=ans
END
SUBROUTINE POLINT4 (XA,YA,X,Y)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C The POLINT4 routine is based on the POLINT routine from "Numerical Recipes",
C but assuming N=4, and ignoring the error estimation.
C suggested by Z. Sullivan.
DIMENSION XA(*),YA(*)
H1=XA(1)-X
H2=XA(2)-X
H3=XA(3)-X
H4=XA(4)-X
W=YA(2)-YA(1)
DEN=W/(H1-H2)
D1=H2*DEN
C1=H1*DEN
W=YA(3)-YA(2)
DEN=W/(H2-H3)
D2=H3*DEN
C2=H2*DEN
W=YA(4)-YA(3)
DEN=W/(H3-H4)
D3=H4*DEN
C3=H3*DEN
W=C2-D1
DEN=W/(H1-H3)
CD1=H3*DEN
CC1=H1*DEN
W=C3-D2
DEN=W/(H2-H4)
CD2=H4*DEN
CC2=H2*DEN
W=CC2-CD1
DEN=W/(H1-H4)
DD1=H4*DEN
DC1=H1*DEN
If((H3+H4).lt.0D0) Then
Y=YA(4)+D3+CD2+DD1
Elseif((H2+H3).lt.0D0) Then
Y=YA(3)+D2+CD1+DC1
Elseif((H1+H2).lt.0D0) Then
Y=YA(2)+C2+CD1+DC1
ELSE
Y=YA(1)+C1+CC1+DC1
ENDIF
RETURN
END
SUBROUTINE POLINT(XA,YA,N,X,Y,DY)
IMPLICIT REAL*8 (A-H,O-Z)
PARAMETER (NMAX=25)
DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
NS=1
DIF=ABS(X-XA(1))
DO 11 I=1,N
DIFT=ABS(X-XA(I))
IF (DIFT.LT.DIF) THEN
NS=I
DIF=DIFT
ENDIF
C(I)=YA(I)
D(I)=YA(I)
11 CONTINUE
Y=YA(NS)
NS=NS-1
DO 13 M=1,N-1
DO 12 I=1,N-M
HO=XA(I)-X
HP=XA(I+M)-X
W=C(I+1)-D(I)
DEN=HO-HP
IF(DEN.EQ.0.)then
print*,xa(i),xa(i+m),x
print*, 'Enter <CR> to continue'
read(*,*)
endif
DEN=W/DEN
D(I)=HP*DEN
C(I)=HO*DEN
12 CONTINUE
IF (2*NS.LT.N-M)THEN
DY=C(NS+1)
ELSE
DY=D(NS)
NS=NS-1
ENDIF
Y=Y+DY
13 CONTINUE
RETURN
END
end module CJ12

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:35 PM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806068
Default Alt Text
(77 KB)

Event Timeline