LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf-5.9.1 / src / uncertainties.f
1 ! -*- F90 -*-
2
3 !--   31/05/2012 by Graeme Watt <Graeme.Watt(at)cern.ch>.
4 !--   Calculation of PDF uncertainties and PDF correlations.
5 !--
6 !--   Use formulae for PDF uncertainties and correlations in:
7 !--    G. Watt, JHEP 1109 (2011) 069 [arXiv:1106.5788 [hep-ph]].
8 !--   Code should distinguish between NNPDF (Monte Carlo approach),
9 !--   Alekhin02/ABKM09/ABM11 (symmetric Hessian approach).  Other
10 !--   PDF sets are assumed to use the asymmetric Hessian approach.
11 !--   List of subroutines in this file:
12 !--    GetMaxNumSets(MaxNumSets)
13 !--    GetPDFUncType(lMonteCarlo,lSymmetric)
14 !--    GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
15 !--    GetPDFuncertainty(values,central,errplus,errminus,errsym)
16 !--    GetPDFuncertaintyM(nset,values,central,errplus,errminus,errsym)
17 !--    GetPDFcorrelation(valuesA,valuesB,correlation)
18 !--    GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
19 !--   Input variables above are nset, values, valuesA, valuesB.
20 !--   Other arguments above are all output variables.
21
22
23 !-- Get flags indicating if Monte Carlo PDF set (NNPDF) and
24 !-- if should compute symmetric errors (NNPDF, Alekhin).
25
26 subroutine GetPDFUncType(lMonteCarlo,lSymmetric)
27   implicit none
28   logical lMonteCarlo,lSymmetric
29   integer nset
30   nset = 1
31   call GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
32   return
33 end subroutine GetPDFUncType
34
35 subroutine GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
36   implicit none
37   include 'parmsetup.inc'
38   logical lMonteCarlo,lSymmetric
39   character*16 name(nmxset)
40   integer nset,nmem(nmxset),ndef(nmxset),mem
41   common/NAME/name,nmem,ndef,mem
42   if ((name(nset).eq.'NNPDF').or.(name(nset).eq.'NNPDFint').or. &
43        (name(nset).eq.'NNPDF20int').or.(name(nset).eq.'NNPDF20intqed')) then ! NNPDF Monte Carlo PDF sets
44      lMonteCarlo = .true.
45      lSymmetric = .true.
46   else if ((name(nset).eq.'A02M').or.(name(nset).eq.'ABKM09').or. &
47        (name(nset).eq.'ABM11')) then ! symmetric eigenvector PDF sets
48      lMonteCarlo = .false.
49      lSymmetric = .true.
50   else ! default: assume asymmetric Hessian eigenvector PDF sets
51      lMonteCarlo = .false.
52      lSymmetric = .false.
53   end if
54 end subroutine GetPDFUncTypeM
55
56
57 !-- Calculate the PDF uncertainty using the appropriate formula for
58 !-- either the Hessian or Monte Carlo approach given an array
59 !-- "values(0:nmem)".  In the Monte Carlo approach, the uncertainty is
60 !-- given by the standard deviation, and the central (average) value
61 !-- is not necessarily "values(0)" for quantities with a non-linear
62 !-- dependence on PDFs.  In the Hessian approach, the central value is
63 !-- the best-fit "values(0)" and the uncertainty is given by either
64 !-- the symmetric or asymmetric formula using eigenvector PDF sets.
65
66 subroutine GetPDFuncertainty(values, &
67      &     central,errplus,errminus,errsym)
68   implicit none
69   integer nset
70   double precision values(0:*),central,errplus,errminus,errsym
71   nset = 1
72   call GetPDFuncertaintyM(nset,values, &
73        &     central,errplus,errminus,errsym)
74   return
75 end subroutine GetPDFuncertainty
76
77
78 subroutine GetPDFuncertaintyM(nset,values, &
79      &     central,errplus,errminus,errsym)
80   implicit none
81   integer nset,nmem,imem
82   double precision values(0:*),central,errplus,errminus,errsym
83   logical lMonteCarlo,lSymmetric
84
85   call numberPDFM(nset,nmem)
86   call GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
87
88   central = 0.D0 ! central value
89   errplus = 0.D0 ! positive uncertainty
90   errminus = 0.D0 ! negative uncertainty
91   errsym = 0.D0 ! symmetrised uncertainty
92
93   if (lMonteCarlo) then ! calculate average and standard deviation
94
95      do imem = 1, nmem
96         central = central + values(imem)
97         errsym = errsym + values(imem)**2
98      end do
99      central = central/nmem ! mean of values
100      errsym = errsym/nmem ! mean of squared values
101      errsym = nmem/(nmem-1.D0)*(errsym-central**2)
102      if (errsym.gt.0.D0) then
103         errsym = sqrt(errsym)
104      else
105         errsym = 0.D0
106      end if
107      errplus = errsym
108      errminus = errsym
109
110   else if (lSymmetric) then ! symmetric Hessian eigenvector PDF sets
111
112      do imem = 1, nmem
113         errsym = errsym + (values(imem)-values(0))**2
114      end do
115      errsym = sqrt(errsym)
116      errplus = errsym
117      errminus = errsym
118      central = values(0)
119
120   else ! default: assume asymmetric Hessian eigenvector PDF sets
121
122      !-- check that nmem is non-zero and even
123      if (nmem.ne.0.and.(mod(nmem,2).eq.0)) then
124         do imem = 1, nmem/2 ! sum over eigenvectors
125            errplus = errplus + max(0.D0, &
126                 &        values(2*imem-1)-values(0), &
127                 &        values(2*imem)-values(0))**2
128            errminus = errminus + max(0.D0, &
129                 &        values(0)-values(2*imem-1), &
130                 &        values(0)-values(2*imem))**2
131            errsym = errsym + (values(2*imem-1)-values(2*imem))**2
132         end do
133         errplus = sqrt(errplus)
134         errminus = sqrt(errminus)
135         errsym = 0.5D0*sqrt(errsym)
136      end if
137      central = values(0)
138
139   end if
140
141   return
142 end subroutine GetPDFuncertaintyM
143
144
145 !-- Calculate the PDF correlation using the appropriate formula for
146 !-- either the Hessian or Monte Carlo approach given two arrays
147 !-- "valuesA(0:nmem)" and "valuesB(0:nmem)".  The correlation can vary
148 !-- between -1 and +1 where values close to {-1,0,+1} mean that the two
149 !-- quantities A and B are {anticorrelated,uncorrelated,correlated}.
150
151 subroutine GetPDFcorrelation(valuesA,valuesB,correlation)
152   implicit none
153   integer nset
154   double precision valuesA(0:*),valuesB(0:*),correlation
155   nset = 1
156   call GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
157   return
158 end subroutine GetPDFcorrelation
159
160 subroutine GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
161   implicit none
162   integer nset,nmem,imem
163   double precision valuesA(0:*),valuesB(0:*),correlation
164   double precision A0,Ap,Am,As,B0,Bp,Bm,Bs
165   logical lMonteCarlo,lSymmetric
166
167   call numberPDFM(nset,nmem)
168   call GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
169
170   call GetPDFuncertaintyM(nset,valuesA,A0,Ap,Am,As)
171   call GetPDFuncertaintyM(nset,valuesB,B0,Bp,Bm,Bs)
172
173   correlation = 0.D0
174
175   if (lMonteCarlo) then ! calculate average and standard deviation
176
177      do imem = 1, nmem
178         correlation = correlation + valuesA(imem)*valuesB(imem)
179      end do
180      correlation = (correlation/nmem - A0*B0)/(As*Bs)*nmem/(nmem-1.D0)
181
182   else if (lSymmetric) then ! symmetric Hessian eigenvector PDF sets
183
184      do imem = 1, nmem
185         correlation = correlation + &
186              & (valuesA(imem)-A0)*(valuesB(imem)-B0)
187      end do
188      correlation = correlation/(As*Bs)
189
190   else ! default: assume asymmetric Hessian eigenvector PDF sets
191
192      !-- check that nmem is non-zero and even
193      if (nmem.ne.0.and.(mod(nmem,2).eq.0)) then
194         do imem = 1, nmem/2 ! sum over eigenvectors
195            correlation = correlation + &
196                 & (valuesA(2*imem-1)-valuesA(2*imem)) * &
197                 & (valuesB(2*imem-1)-valuesB(2*imem))
198         end do
199         correlation = correlation/(4.D0*As*Bs)
200      end if
201
202   end if
203
204   return
205 end subroutine GetPDFcorrelationM