]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-par-129.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes needed by gfortran, it doesn't perform lazy evaluation (Piotr)
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-par-129.f
CommitLineData
9ef1c2d9 1c 15/02/2005 epos 1.03
2
3
4c----------------------------------------------------------------------
5 subroutine paramini(imod)
6c----------------------------------------------------------------------
7c Set parameters of the parametrisation of the eikonals.
8c
9c xDfit=Sum(i=0,1)(alpD(i)*xp**betDp(i)*xm**betDpp(i)*s**betD(i)
10c *xs**(gamD(i)*b2)*exp(-b2/delD(i))
11c
12c Parameters stored in /Dparam/ (epos.inc)
13c if imod=0, do settings only for iclpro, if imod=1, do settings
14c for iclpro=2 and iclpro
15c----------------------------------------------------------------------
16
17 include 'epos.inc'
18 include 'epos.incems'
19 include 'epos.incsem'
20 include 'epos.incpar'
21 double precision PhiExact,y
22 call utpri('parini',ish,ishini,3)
23
24c Initialisation of the variables
25
26 call Class('paramini ')
27
28c Variables used for xparg (only graphical application)
29
30 spmin=4.*q2min !Definition of spmin in psvin (epos-sha)
31 sfshlim=5.*spmin !transition energy for soft->hard
32 emaxDf=engy !energy for fit subroutines
33 smaxDf=emaxDf**2 !energy squared for fit subroutines
34c sfshlim=100.
35 nptf=10 !number of point for the fit
36 xmaxDf=1.d0 !maximum limit for the fit
37 xminDf=1.d0/dble(smaxDf) !minimum limit
38 xshmin=3.d-2 !minimum limit for sh variance fit
39 if(smaxDf.lt.sfshlim) xshmin=1.d0
40 xfitmin=0.1d0 !sh fitted between f(1) and f(1)*xfitmin
41 if(smaxDf.lt.sfshlim) xfitmin=1.d0
42
43 bmaxDf=2. !maximum b for variance fit
44
45 idxDmin=idxD0 !minimum indice for parameters
46 ntymin=ntymi !minimum indice for diagram type
47
48 ucfpro=utgam1(1.+alplea(iclpro))
49 ucftar=utgam1(1.+alplea(icltar))
50
51
52
53 iiiclegy=iclegy
54
55c for pi or K - p crosse section calculation, we need alpD, ... for
56c iclpro=1 or 3 and iclpro=2
57 iiiclpro1=iclpro
58 iiidirec=1
59 if(imod.eq.0)then
60 iiiclpro2=iclpro
61 else
62 iiiclpro2=2
63 if(iiiclpro1.lt.iiiclpro2)iiidirec=-1
64 endif
65 iclprosave=iclpro
66
67 do iiipro=iiiclpro2,iiiclpro1,iiidirec
68
69 iclpro=iiipro
70
71 if(ish.ge.4)write(ifch,*)'gamini'
72 & ,iscreen,iclpro,icltar,iclegy,engy,sfshlim
73
74 if(isetcs.le.1)then !if set mode, do fit
75
76c First try for fit parameters
77
78c linear fit of the 2 components of G as a function of x
79 call pompar(alpsf,betsf,0) !soft (taking into account hard)
80 call pompar(alpsh,betsh,1) !sh
81c Gaussian fit of the 2 components of G as a function of x and b
82 call variance(delsf,gamsf,0)
83 call variance(delsh,gamsh,1)
84 gamsf=max(0.,gamsf)
85 gamsh=max(0.,gamsh)
86
87
88c Fit GFF
89
90
91c fit parameters
92 numminDf=3 !minimum indice
93 numparDf=4 !maximum indice
94 betac=100. !temperature for chi2
95 betae=1. !temperature for error
96 fparDf=0.8 !variation amplitude for range
97
98 nmcxDf=20000 !number of try for x fit
99
100c starting values
101
102 parDf(1,3)=alpsf
103 parDf(2,3)=betsf
104 parDf(3,3)=alpsh
105 parDf(4,3)=betsh
106
107 if(smaxDf.ge.3.*sfshlim)then
108 call paramx !alpD and betD
109
110 call pompar(alpsf,betsf,-1) !soft (taking into account hard)
111 parDf(1,3)=alpsf
112 parDf(2,3)=max(-0.95+alppar,betsf)
113
114 endif
115
116 alpsf=parDf(1,3)
117 betsf=parDf(2,3)
118 alpsh=parDf(3,3)
119 betsh=parDf(4,3)
120
121
122 else !else parameters from table (inirj)
123 nbpsf=iDxD0
124 if(iclegy2.gt.1)then
125 al=1.+(log(engy)-log(egylow))/log(egyfac) !energy class
126 i2=min(iiiclegy+1,iclegy2)
127 i1=i2-1
128 else
129 i1=iclegy
130 i2=iclegy
131 al=float(iclegy)
132 endif
133 dl=al-i1
134 dl1=max(0.,1.-dl)
135 !linear interpolation
136 alpsf=alpDs(nbpsf,i2,iclpro,icltar)*dl
137 & +alpDs(nbpsf,i1,iclpro,icltar)*dl1
138 alpsh=alpDs(1,i2,iclpro,icltar)*dl
139 & +alpDs(1,i1,iclpro,icltar)*dl1
140 betsf=betDs(nbpsf,i2,iclpro,icltar)*dl
141 & +betDs(nbpsf,i1,iclpro,icltar)*dl1
142 betsh=betDs(1,i2,iclpro,icltar)*dl
143 & +betDs(1,i1,iclpro,icltar)*dl1
144 gamsf=gamDs(nbpsf,i2,iclpro,icltar)*dl
145 & +gamDs(nbpsf,i1,iclpro,icltar)*dl1
146 gamsh=gamDs(1,i2,iclpro,icltar)*dl
147 & +gamDs(1,i1,iclpro,icltar)*dl1
148 delsf=delDs(nbpsf,i2,iclpro,icltar)*dl
149 & +delDs(nbpsf,i1,iclpro,icltar)*dl1
150 delsh=delDs(1,i2,iclpro,icltar)*dl
151 & +delDs(1,i1,iclpro,icltar)*dl1
152
153
154c For the Plots
155 parDf(1,3)=alpsf
156 parDf(2,3)=betsf
157 parDf(3,3)=alpsh
158 parDf(4,3)=betsh
159
160 endif
161
162c if energy too small to have semi-hard interaction -> only soft diagram
163
164 if(smaxDf.lt.sfshlim.and.idxD0.eq.0)then !no hard: soft+hard=soft
165 alpsf=alpsf/2.
166 alpsh=alpsf
167 betsh=betsf
168 gamsh=gamsf
169 delsh=delsf
170 endif
171
172
173
174c Print results
175
176 if(ish.ge.4)then
177 write(ifch,*)"parameters for iclpro:",iclpro
178 write(ifch,*)"alp,bet,gam,del sf:",alpsf,betsf,gamsf,delsf
179 write(ifch,*)"alp,bet,gam,del sh:",alpsh,betsh,gamsh,delsh
180 endif
181
182
183c Record parameters
184
185
186 alpD(idxD0,iclpro,icltar)=alpsf
187 alpDp(idxD0,iclpro,icltar)=0.
188 alpDpp(idxD0,iclpro,icltar)=0.
189 betD(idxD0,iclpro,icltar)=betsf
190 betDp(idxD0,iclpro,icltar)=betsf
191 betDpp(idxD0,iclpro,icltar)=betsf
192 gamD(idxD0,iclpro,icltar)=gamsf
193 delD(idxD0,iclpro,icltar)=delsf
194
195 alpD(1,iclpro,icltar)=alpsh
196 alpDp(1,iclpro,icltar)=0.
197 alpDpp(1,iclpro,icltar)=0.
198 betD(1,iclpro,icltar)=betsh
199 betDp(1,iclpro,icltar)=betsh
200 betDpp(1,iclpro,icltar)=betsh
201 gamD(1,iclpro,icltar)=gamsh
202 delD(1,iclpro,icltar)=delsh
203
204 if(iomega.lt.2.and.alpdif.ne.1.)then
205 alpDp(2,iclpro,icltar)=0.
206 alpDpp(2,iclpro,icltar)=0.
207 betD(2,iclpro,icltar)=0. !max(0.,betsf)
208c alpdifs=alpdif
209c alpdif=0.99
210 betDp(2,iclpro,icltar)=-alpdif+alppar
211 betDpp(2,iclpro,icltar)=-alpdif+alppar
212c alpD(2,iclpro,icltar)=(alpsf+alpsh)*wdiff(iclpro)*wdiff(icltar)
213 alpD(2,iclpro,icltar)=wdiff(iclpro)*wdiff(icltar)
214 & /utgam1(1.-alpdif)**2
215 & *utgam1(2.-alpdif+alplea(iclpro))
216 & *utgam1(2.-alpdif+alplea(icltar))
217 & /chad(iclpro)/chad(icltar)
218c alpdif=alpdifs
219 gamD(2,iclpro,icltar)=0.
220 delD(2,iclpro,icltar)=4.*.0389*(gwidth*(r2had(iclpro)
221 & +r2had(icltar))+slopoms*log(smaxDf))
222 else
223 alpD(2,iclpro,icltar)=0.
224 alpDp(2,iclpro,icltar)=0.
225 alpDpp(2,iclpro,icltar)=0.
226 betD(2,iclpro,icltar)=0.
227 betDp(2,iclpro,icltar)=0.
228 betDpp(2,iclpro,icltar)=0.
229 gamD(2,iclpro,icltar)=0.
230 delD(2,iclpro,icltar)=1.
231 endif
232 if(ish.ge.4)write(ifch,*)"alp,bet,betp,del dif:"
233 & ,alpD(2,iclpro,icltar),betD(2,iclpro,icltar)
234 & ,betDp(2,iclpro,icltar),delD(2,iclpro,icltar)
235
236
237 bmxdif(iclpro,icltar)=conbmxdif() !important to do it before kfit, because it's used in.
238c call Kfit(-1) !xkappafit not used : if arg=-1, set xkappafit to 1
239 if(isetcs.le.1)then
240 if(isetcs.eq.0)then
241 call Kfit(-1) !xkappafit not used : if arg=-1, set xkappafit to 1)
242 else
243c call Kfit(-1) !xkappafit not used : if arg=-1, set xkappafit to 1)
244 call Kfit(iclegy)
245 endif
246c for plots record alpDs, betDs, etc ...
247 alpDs(idxD0,iclegy,iclpro,icltar)=alpsf
248 alpDps(idxD0,iclegy,iclpro,icltar)=0.
249 alpDpps(idxD0,iclegy,iclpro,icltar)=0.
250 betDs(idxD0,iclegy,iclpro,icltar)=betsf
251 betDps(idxD0,iclegy,iclpro,icltar)=betsf
252 betDpps(idxD0,iclegy,iclpro,icltar)=betsf
253 gamDs(idxD0,iclegy,iclpro,icltar)=gamsf
254 delDs(idxD0,iclegy,iclpro,icltar)=delsf
255
256 alpDs(1,iclegy,iclpro,icltar)=alpsh
257 alpDps(1,iclegy,iclpro,icltar)=0.
258 alpDpps(1,iclegy,iclpro,icltar)=0.
259 betDs(1,iclegy,iclpro,icltar)=betsh
260 betDps(1,iclegy,iclpro,icltar)=betsh
261 betDpps(1,iclegy,iclpro,icltar)=betsh
262 gamDs(1,iclegy,iclpro,icltar)=gamsh
263 delDs(1,iclegy,iclpro,icltar)=delsh
264 endif
265
266 enddo
267
268 if(iclpro.ne.iclprosave)stop'problem in parini with iclpro'
269
270 if(ish.ge.4)then !check PhiExact value for x=1
271 y=PhiExact(1.,1.d0,1.d0,smaxDf,0.)
272 write(ifch,*)'PhiExact=',y
273 endif
274
275 call utprix('parini',ish,ishini,3)
276
277 return
278 end
279
280c----------------------------------------------------------------------
281 subroutine Class(text)
282c----------------------------------------------------------------------
283 include 'epos.inc'
284 include 'epos.incpar'
285 parameter (eps=1.e-5) !to correct for precision problem)
286 character*10 text
287 if(iclegy1.eq.iclegy2)then
288 iclegy=iclegy1
289 else
290 iclegy=1+int( (log(engy)-log(egylow))/log(egyfac) + eps ) !energy class
291 if(iclegy.gt.iclegy2)then
292 write(ifch,*)'***********************************************'
293 write(ifch,*)'Warning in ',text
294 write(ifch,*)'Energy above the range used for the fit of D:'
295 write(ifch,*)egylow*egyfac**(iclegy1-1),egylow*egyfac**iclegy2
296 write(ifch,*)'***********************************************'
297 iclegy=iclegy2
298 endif
299 if(iclegy.lt.iclegy1)then
300 write(ifch,*)'***********************************************'
301 write(ifch,*)'Warning in ',text
302 write(ifch,*)'Energy below the range used for the fit of D:'
303 write(ifch,*)egylow*egyfac**(iclegy1-1),egylow*egyfac**iclegy2
304 write(ifch,*)'***********************************************'
305 iclegy=iclegy1
306 endif
307 endif
308 end
309
310c----------------------------------------------------------------------
311 subroutine param
312c----------------------------------------------------------------------
313c Set the parameter of the parametrisation of the eikonale.
314c We group the parameters into 4 array with a dimension of idxD1(=1)
315c to define xDfit (see below).
316c
317c xDfit=Sum(i,0,1)(alpD(i)*xp**betDp(i)*xm**betDpp(i)*s**betD(i)
318c *xs**(gamD(i)*b2)*exp(-b2/delD(i))
319c
320c subroutine used for tabulation.
321c----------------------------------------------------------------------
322 include 'epos.inc'
323 include 'epos.incsem'
324 include 'epos.incpar'
325
326c Initialisation of the variables
327
328 emaxDf=egyfac**(iclegy-1)*egylow
329 smaxDf=emaxDf**2
330
331 spmin=4.*q2min !Definition of spmin in psvin (epos-sha)
332 sfshlim=5.*spmin
333 nptf=10
334 xmaxDf=1.d0
335 xminDf=1d0/dble(smaxDf)
336 xshmin=3.d-2 !minimum limit for sh variance fit
337 if(smaxDf.lt.sfshlim) xshmin=1.d0
338 xfitmin=0.1d0 !sh fitted between f(1) and f(1)*xfitmin
339 if(smaxDf.lt.sfshlim) xfitmin=1.d0
340 bmaxDf=2.
341
342
343 if(idxD0.ne.0.and.idxD1.ne.1) stop "* idxD0/1 are not good! *"
344
345 engytmp=engy
346 engy=emaxDf
347
348c Initialisation of the parameters
349
350 do i=1,nbpf
351 do j=1,4
352 parDf(i,j)=1.
353 enddo
354 enddo
355
356c.......Calculations of the parameters
357
358c First try for fit parameters
359
360c linear fit of the 2 components of G as a function of x
361 call pompar(alpsf,betsf,0) !soft
362 call pompar(alpsh,betsh,1) !sh
363c Gaussian fit of the 2 components of G as a function of x and b
364 call variance(delsf,gamsf,0)
365 call variance(delsh,gamsh,1)
366 gamsf=max(0.,gamsf)
367 gamsh=max(0.,gamsh)
368
369c Fit GFF
370
371c fit parameters
372 numminDf=3 !minimum indice
373 numparDf=4 !maximum indice
374 betac=100. !temperature for chi2
375 betae=1. !temperature for error
376 fparDf=0.8 !variation amplitude for range
377
378 nmcxDf=20000 !number of try for x fit
379
380c starting values
381
382 parDf(1,3)=alpsf
383 parDf(2,3)=betsf
384 parDf(3,3)=alpsh
385 parDf(4,3)=betsh
386
387 if(smaxDf.ge.3.*sfshlim)then
388 call paramx !alpD and betD
389
390 call pompar(alpsf,betsf,-1) !soft (taking into account hard)
391 parDf(1,3)=alpsf
392 parDf(2,3)=max(-0.95+alppar,betsf)
393
394 endif
395
396 alpsf=parDf(1,3)
397 betsf=parDf(2,3)
398 alpsh=parDf(3,3)
399 betsh=parDf(4,3)
400
401 if(ish.ge.4)then
402 write(ifch,*)"param: fit parameters :",iscreen,iclpro,icltar
403 * ,iclegy,engy
404 write(ifch,*)"alp,bet,gam,del sf:",alpsf,betsf,gamsf,delsf
405 write(ifch,*)"alp,bet,gam,del sh:",alpsh,betsh,gamsh,delsh
406 endif
407
408 alpDs(idxD0,iclegy,iclpro,icltar)=alpsf
409 alpDps(idxD0,iclegy,iclpro,icltar)=betsf
410 alpDpps(idxD0,iclegy,iclpro,icltar)=0.
411 betDs(idxD0,iclegy,iclpro,icltar)=betsf
412 betDps(idxD0,iclegy,iclpro,icltar)=betsf
413 betDpps(idxD0,iclegy,iclpro,icltar)=betsf
414 gamDs(idxD0,iclegy,iclpro,icltar)=gamsf
415 delDs(idxD0,iclegy,iclpro,icltar)=delsf
416
417 alpDs(1,iclegy,iclpro,icltar)=alpsh
418 alpDps(1,iclegy,iclpro,icltar)=betsh
419 alpDpps(1,iclegy,iclpro,icltar)=0.
420 betDs(1,iclegy,iclpro,icltar)=betsh
421 betDps(1,iclegy,iclpro,icltar)=betsh
422 betDpps(1,iclegy,iclpro,icltar)=betsh
423 gamDs(1,iclegy,iclpro,icltar)=gamsh
424 delDs(1,iclegy,iclpro,icltar)=delsh
425
426 engy=engytmp
427
428 return
429 end
430
431
432c----------------------------------------------------------------------
433
434 subroutine pompar(alpha,beta,iqq)
435
436c----------------------------------------------------------------------
437c Return the power beta and the factor alpha of the fit of the eikonal
438c of a pomeron of type iqq : D(X)=alpha*(X)**beta
439c----------------------------------------------------------------------
440
441 include 'epos.inc'
442 include 'epos.incsem'
443 include 'epos.incpar'
444 double precision X,D1,D0,X0,D,droot
445 double precision Dsoftshval,xmax
446 dimension xlnXs(maxdataDf),xlnD(maxdataDf),sigma(maxdataDf)
447
448 do i=1,nptf
449 sigma(i)=1.e-2
450 enddo
451
452 if(iqq.le.0) then
453
454 iscr=iqq
455 xmax=min(0.1d0,10.d0*xminDf)
456 X0=xminDf
457 if(ish.ge.4)write (ifch,*) 'pompar (0) x0,xmax=',X0,xmax
458
459 do i=0,nptf-1
460 X=X0
461 if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
462 D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
463 if(D.eq.1.d-10)then
464 write(ifch,*)
465 & "Warning in pompar ! Dsoftshval(0) could be negative"
466 sigma(i+1)=1.e5
467 endif
468 xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
469 xlnD(i+1)=real(dlog(D))
470 enddo
471
472
473c Fit of D(X) between X0 and xmax
474
475 call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
476 if(beta.gt.10.)beta=10.
477
478 alpha=real(Dsoftshval(real(X0)*smaxDf,X0,0.d0,0.,iscr))
479 & *(real(X0)*smaxDf)**(-beta)
480
481
482 elseif(iqq.eq.1.and.xfitmin.ne.1.d0) then
483
484 iscr=2
485 xmax=1.d0
486
487c Definition of D0=D(X0)
488
489 D1=Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr)
490
491 D0=xfitmin*D1
492
493c Calculation of X0 and D(X)
494
495 X0=droot(D0,D1,xmax,iscr)
496 if(ish.ge.4)write (ifch,*) 'pompar (1) x0,xmax=',X0,xmax
497
498 do i=0,nptf-1
499 X=X0
500 if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
501 D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
502 if(D.eq.1.d-10)then
503 write(ifch,*)
504 & "Warning in pompar ! Dsoftshval(1) could be negative"
505 sigma(i+1)=1.e5
506 endif
507 xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
508 xlnD(i+1)=real(dlog(D))
509 enddo
510
511
512c Fit of D(X) between X0 and xmax
513
514 call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
515 if(beta.gt.10.)beta=10.
516
517 alpha=real(Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr))
518 & *(real(xmax)*smaxDf)**(-beta)
519
520
521 elseif(iqq.eq.10.and.xfitmin.ne.1.d0) then ! iqq=10
522
523 iscr=0
524 xmax=1.d0 !2.d0/max(2.d0,dlog(dble(smaxDf)/1.d3))
525
526c Definition of D0=D(X0)
527
528 D1=Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr)
529
530 D0=xfitmin*D1
531
532c Calculation of X0 and D(X)
533
534 X0=droot(D0,D1,xmax,iscr)
535 if(ish.ge.4)write (ifch,*) 'pompar (1) x0,xmax=',X0,xmax
536
537 do i=0,nptf-1
538 X=X0
539 if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
540 D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
541 if(D.eq.1.d-10)then
542 write(ifch,*)
543 & "Warning in pompar ! Dsoftshval(10) could be negative"
544 sigma(i+1)=1.e5
545 endif
546 xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
547 xlnD(i+1)=real(dlog(D))
548 enddo
549
550
551c Fit of D(X) between X0 and xmax
552
553 call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
554 if(beta.gt.10.)beta=10.
555
556 alpha=real(Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr))
557 & *(real(xmax)*smaxDf)**(-beta)
558
559
560
561 else !iqq=-1 or iqq=1 and xfitmin=1
562
563c Calculation of X0 and D(X)
564 iscr=0
565
566 X0=1.d0/dble(smaxDf)
567 xmax=max(2.d0/dble(smaxDf),
568 & min(max(0.03d0,dble(smaxDf)/2.d5),0.1d0))
569
570 if(ish.ge.4)write (ifch,*) 'pompar (-1) x0,xmax=',X0,xmax
571
572 do i=0,nptf-1
573 X=X0
574 if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
575 D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
576 if(D.eq.1.d-10)then
577 write(ifch,*)
578 & "Warning in pompar ! Dsoftshval(-1) could be negative"
579 sigma(i+1)=1.e5
580 endif
581 xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
582 xlnD(i+1)=real(dlog(D))
583 enddo
584
585c Fit of D(X) between X0 and xmax
586
587 call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
588 if(beta.gt.10.)beta=10.
589
590
591 alpha=real(Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr))
592 & *(real(xmax)*smaxDf)**(-beta)
593 endif
594
595 if(ish.ge.4)write(ifch,*) '%%%%%%%%%%%%% pompar %%%%%%%%%%%%%'
596 if(ish.ge.4)write(ifch,*) 'alpD ini =',alpha,' betD ini=',beta
597
598 return
599 end
600
601c----------------------------------------------------------------------
602
603 double precision function droot(d0,d1,xmax,iscr)
604
605c----------------------------------------------------------------------
606c Find x0 which gives f(x0)=D(x0*S)-d0=0
607c iqq=0 soft pomeron
608c iqq=1 semi-hard pomeron
609c----------------------------------------------------------------------
610 include 'epos.inc'
611 include 'epos.incsem'
612 include 'epos.incpar'
613 double precision Dsoftshval,d0,d1,d2,x0,x1,x2,f0,f1,f2,xmax
614 parameter (kmax=1000)
615
616
617 k=0
618 x0=min(xfitmin,100.d0*xminDf)
619 x1=xmax
620 5 d2=dabs(Dsoftshval(real(x0)*smaxDf,x0,0.d0,0.,iscr))
621 if(d2.lt.1.d-10.and.x0.lt.x1)then
622 x0=dsqrt(x0*x1)
623c write(ifch,*)"droot",x0,x1,d0,d1,d2
624 goto 5
625 elseif(d2.gt.d0)then
626 droot=max(x0,xfitmin)
627c write(ifch,*)"droot",x0,x1,d0,d1,d2
628 return
629 endif
630 f0=d2-d0
631 f1=d1-d0
632 if(f0*f1.lt.0.d0)then
633
634
635 10 x2=dsqrt(x0*x1)
636 d2=dabs(Dsoftshval(real(x2)*smaxDf,x2,0.d0,0.,iscr))
637 f2=d2-d0
638 k=k+1
639c write (ifch,*) '******************* droot **************'
640c write (ifch,*) x0,x1,x2,f0,f1,f2
641
642 if (f0*f2.lt.0.D0) then
643 x1=x2
644 f1=f2
645 else
646 x0=x2
647 f0=f2
648 endif
649
650 if (dabs((x1-x0)/x1).gt.(1.D-5).and.k.le.kmax.and.x1.ne.x0) then
651 goto 10
652 else
653 if (k.gt.kmax) then
654 write(ifch,*)'??? Warning in Droot: Delta=',dabs((x1-x0)/x1)
655c.........stop 'Error in Droot, too many steps'
656 endif
657 droot=dsqrt(x1*x0)
658 endif
659
660 else
661 droot=dsqrt(x1*x0)
662 endif
663
664 return
665 end
666
667c----------------------------------------------------------------------
668
669 double precision function drootom(d0,dmax,bmax,eps)
670
671c----------------------------------------------------------------------
672c Find b0 which gives f(b0)=(1-exp(-om(b0,iqq)))/dmax-d0=0
673 include 'epos.inc'
674 include 'epos.incsem'
675 include 'epos.incpar'
676 double precision om1intbc,d0,d1,d2,f0,f1,f2,dmax
677 parameter (kmax=1000)
678
679
680 k=0
681 b0=0.
682 b1=bmax
683 d2=(1.d0-exp(-om1intbc(b1)))/dmax
684 if(d2.gt.d0)then
685 drootom=b1
686c write(*,*)"drootom exit (1)",b0,b1,d0,d1,d2
687 return
688 endif
689 d1=(1.d0-exp(-om1intbc(b0)))/dmax
690 f0=d1-d0
691 f1=d2-d0
692 if(f0*f1.lt.0.d0)then
693
694
695 10 b2=0.5*(b0+b1)
696 d2=(1.d0-dexp(-om1intbc(b2)))/dmax
697 f2=d2-d0
698 k=k+1
699c write (*,*) '******************* drootom **************'
700c write (*,*) b0,b1,b2,f0,f1,f2
701
702 if (f1*f2.lt.0.D0) then
703 b0=b2
704 f0=f2
705 else
706 b1=b2
707 f1=f2
708 endif
709
710 if (abs(f2).gt.eps.and.k.le.kmax.and.b1.ne.b0) then
711 goto 10
712 else
713 if (k.gt.kmax) then
714 write(ifch,*)'??? Warning in Drootom: Delta=',abs((b1-b0)/b1)
715c.........stop 'Error in Droot, too many steps'
716 endif
717 drootom=0.5*(b1+b0)
718 endif
719
720 else
721c write(*,*)"drootom exit (2)",b0,b1,d0,d1,d2
722 drootom=0.5*(b1+b0)
723 endif
724
725 return
726 end
727
728c----------------------------------------------------------------------
729 subroutine variance(r2,alp,iqq)
730c----------------------------------------------------------------------
731c fit sigma2 into : 1/sigma2(x)=1/r2-alp*log(x*s)
732c iqq=0 -> soft pomeron
733c iqq=1 -> semi-hard pomeron
734c iqq=2 -> sum
735c----------------------------------------------------------------------
736 include 'epos.inc'
737 include 'epos.incsem'
738 include 'epos.incpar'
739 dimension Xs(maxdataDf),vari(maxdataDf),sigma(maxdataDf)
740 double precision X,X0,xmax
741
742 do i=1,nptf
743 sigma(i)=1.e-2
744 enddo
745
746 if(iqq.eq.0.or.xshmin.gt.0.95d0)then
747 X0=xminDf
748 xmax=xmaxDf
749 elseif(iqq.eq.2)then
750 X0=xshmin
751 xmax=xmaxDf
752 else
753 X0=.1d0/dlog(max(dexp(2.d0),dble(smaxDf)/1.d3))
754 if(smaxDf.lt.100.*q2min)X0=.95d0
755 xmax=xmaxDf
756 endif
757 if(iqq.ne.3.and.iqq.ne.4)then
758
759 do i=0,nptf-1
760 X=X0
761 if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
762 Xs(i+1)=log(real(X)*smaxDf)
763 sig2=sigma2(X,iqq)
764 if(sig2.le.0.) call utstop
765 & ('In variance, initial(1) sigma2 not def!&')
766 vari(i+1)=1./sig2
767 enddo
768
769c Fit of the variance of D(X,b) between X0 and xmaxDf
770
771 call fit(Xs,vari,nptf,sigma,0,tr2,talp)
772 r2=1./tr2
773 alp=-talp
774c in principle, the formula to convert 1/(del+eps*log(sy)) into
775c 1/del*(1-eps/del*log(sy)) is valid only if eps/del*log(sy)=alp*r2*log(sy)
776c is small. In practice, since the fit of G(x) being an approximation, each
777c component of the fit should not be taken separatly but we should consider
778c G as one function. Then it works even with large alp (gamD).
779c ttt=alp*r2*log(smaxDf)
780c if(ttt.gt.0.5)
781c & write(ifmt,*)'Warning, G(b) parametrization not optimal : ',
782c & 'gamD too large compared to delD !',ttt
783 else
784 if(iqq.eq.3)r2=sigma2(xmaxDf,3)
785 if(iqq.eq.4)r2=sigma2(xshmin,3)
786 if(r2.le.0.) call utstop
787 &('In variance, initial(2) sigma2 not def!&')
788 alp=0.
789 endif
790
791 if(ish.ge.4)then
792 write(ifch,*) '%%%%%%%%%% variance ini %%%%%%%%%%%%'
793 write(ifch,*) 'X0=',X0
794 write(ifch,*) 'delD ini=',r2
795 write(ifch,*) 'gamD ini=',alp
796 endif
797
798 return
799 end
800
801
802
803c----------------------------------------------------------------------
804
805 function sigma2(x,iqq)
806
807c----------------------------------------------------------------------
808c Return the variance for a given x of :
809c For G :
810c iqq=0 the soft pomeron
811c iqq=1 the semi-hard and valence quark pomeron
812c iqq=2 the sum
813c----------------------------------------------------------------------
814
815 include 'epos.inc'
816 include 'epos.incpar'
817 double precision x,Dsoftshval,sfsh,om51p,eps,range,sig2!,omNpuncut
818 external varifit
819 double precision varifit,Db(maxdataDf),bf(maxdataDf)
820
821 bmax=bmaxDf
822 sig2=bmax*0.5
823 bmin=-bmax
824 eps=1.d-10
825 ierro=0
826
827 if(iqq.eq.0)then
828 range=sig2
829 sfsh=om51p(real(x)*smaxDf,x,0.d0,0.,0)
830 if (dabs(sfsh).gt.eps) then
831 do i=0,nptf-1
832 bf(i+1)=dble(bmin+real(i)*(bmax-bmin)/real(nptf-1))
833 Db(i+1)=om51p(real(x)*smaxDf,x,0.d0,real(bf(i+1)),0)/sfsh
834 enddo
835 else
836 ierro=1
837 endif
838 elseif(iqq.eq.1.and.xshmin.lt..95d0)then
839 range=sig2
840 sfsh=0.d0
841 do j=1,4
842 sfsh=sfsh+om51p(real(x)*smaxDf,x,0.d0,0.,j)
843 enddo
844 if (dabs(sfsh).gt.eps) then
845 do i=0,nptf-1
846 bf(i+1)=dble(bmin+real(i)*(bmax-bmin)/real(nptf-1))
847 Db(i+1)=0.d0
848 do j=1,4
849 Db(i+1)=Db(i+1)+om51p(real(x)*smaxDf,x,0.d0,real(bf(i+1)),j)
850 enddo
851 Db(i+1)=Db(i+1)/sfsh
852 enddo
853 else
854 ierro=1
855 endif
856 else
857 sig2=2.d0*sig2
858 range=sig2
859 iscr=0
860 sfsh=Dsoftshval(real(x)*smaxDf,x,0.d0,0.,iscr)
861 if (dabs(sfsh).gt.eps) then
862 do i=0,nptf-1
863 bf(i+1)=dble(bmin+real(i)*(bmax-bmin)/real(nptf-1))
864 Db(i+1)=Dsoftshval(real(x)*smaxDf,x,0.d0,real(bf(i+1)),iscr)
865 & /sfsh
866 enddo
867 else
868 ierro=1
869 endif
870 endif
871
872c Fit of D(X,b) between -bmaxDf and bmaxDf
873
874 if(ierro.ne.1)then
875 nptft=nptf
876 call minfit(varifit,bf,Db,nptft,sig2,range)
877 sigma2=real(sig2)
878 else
879 sigma2=0.
880 endif
881
882 return
883 end
884
885c----------------------------------------------------------------------
886
887 subroutine paramx
888
889c----------------------------------------------------------------------
890c updates the 4 parameters alpsf,betsf,alpsh,betsh by fitting GFF
891c parDf(1,3) parDf(2,3) ... alp, bet soft
892c parDf(3,3) parDf(4,3) ... alp, bet semihard
893c----------------------------------------------------------------------
894
895 include 'epos.inc'
896 include 'epos.incpar'
897
898 double precision Dsoftshpar
899
900 external Dsoftshpar
901
902 dimension range(nbpf)
903
904 call givedatax
905
906 !determine parameter range
907 do i=numminDf,numparDf
908 range(i)=fparDf*parDf(i,3)
909 parDf(i,1)=parDf(i,3)-range(i)
910 parDf(i,2)=parDf(i,3)+range(i)
911 enddo
912
913
914 ! write(ifch,*) '%%%%%%%%%%%%%%%%%%% fitx %%%%%%%%%%%%%%%%%%%%%%%'
915
916 call fitx(Dsoftshpar,nmcxDf,chi2,err)
917
918 ! write(ifch,*) 'chi2=',chi2
919 ! write(ifch,*) 'err=',err
920 ! write(ifch,*) 'alpD(1)=',parDf(1,3),' betD(1)=',parDf(2,3)
921 ! write(ifch,*) 'alpD(2)=',parDf(3,3),' betD(2)=',parDf(4,3)
922
923 return
924 end
925
926
927c----------------------------------------------------------------------
928 subroutine givedatax
929c----------------------------------------------------------------------
930
931 include 'epos.inc'
932 include 'epos.incsem'
933 include 'epos.incpar'
934 double precision X,X0,X1,Dsoftshval,Xseuil
935
936 numdataDf=nptf
937
938 X0=max(xminDf,dble(sfshlim/smaxDf))
939 X1=xmaxDf
940 Xseuil=.0001d0
941
942c Fit of G(X) between X0 and X1
943 do i=0,nptf-1
944 X=X0
945 if (i.ne.0) X=X*(X1/X0)**(dble(i)/dble(nptf-1))
946 datafitD(i+1,2)=max(1.e-10,
947 & real(Dsoftshval(real(X)*smaxDf,X,0.d0,0.,1)))
948 datafitD(i+1,1)=real(X)
949 datafitD(i+1,3)=1.
950 if (X.gt.Xseuil) datafitD(i+1,3)=exp((Xseuil/X)-1.)
951 enddo
952
953 return
954 end
955
956
957
958
959
960c----------------------------------------------------------------------
961
962 function sigma1i(x)
963
964c----------------------------------------------------------------------
965c Return the variance of the sum of the soft pomeron and the semi-hard
966c pomeron for a given x.
967c----------------------------------------------------------------------
968
969 include 'epos.inc'
970 include 'epos.incpar'
971 double precision x,Dsoftshval,Dint
972
973
974 iscr=0
975 Dint=Dsoftshval(real(x)*smaxDf,x,0.d0,0.,iscr)
976
977 sigma1i=0.
978 if(Dint.ne.0.)
979 &sigma1i=real(-1.d0/dlog(Dsoftshval(real(x)*smaxDf,x,0.d0,1.,iscr)
980 & /Dint))
981
982
983 return
984 end
985
986
987c----------------------------------------------------------------------
988
989 SUBROUTINE minfit(func,x,y,ndata,a,range)
990
991c----------------------------------------------------------------------
992c Given a set of data points x(1:ndata),y(1:ndata), and the range of
993c the parameter a, fit it on function func by minimizing chi2.
994c In input a define the expected value of a, and on output they
995c correspond to the fited value.
996c ---------------------------------------------------------------------
997 include 'epos.inc'
998 double precision x(ndata),y(ndata),func,a,range,Smin,Som,a1,a2,eps
999 *,amin,rr,yp
1000 parameter (eps=1.d-5)
1001 external func
1002
1003
1004 Smin=1.d20
1005 amin=a
1006
1007
1008
1009 a1=a-range
1010 a2=a+range
1011
1012 do j=1,2000
1013 rr=dble(rangen())
1014 a=a1+(a2-a1)*rr
1015 k=0
1016
1017 10 if(a.lt.0.d0.and.k.lt.100) then
1018 rr=dble(rangen())
1019 a=a1+(a2-a1)*rr
1020 k=k+1
1021 goto 10
1022 endif
1023 if(k.ge.100) call utstop
1024 &('Always negative variance in minfit ...&')
1025
1026 Som=0.d0
1027 do k=1,ndata
1028 yp=min(1.d10,func(x(k),a)) !proposal function
1029 Som=Som+(yp-y(k))**2.d0
1030 enddo
1031 if(Som.lt.Smin)then
1032
1033 if(Smin.lt.1.)then
1034 if(a.gt.amin)then
1035 a1=amin
1036 else
1037 a2=amin
1038 endif
1039 endif
1040 amin=a
1041 Smin=Som
1042 endif
1043 if(Smin.lt.eps)goto 20
1044 enddo
1045
1046 20 continue
1047 a=amin
1048
1049 return
1050 end
1051
1052
1053
1054c----------------------------------------------------------------------
1055 subroutine fitx(func,nmc,chi2,err)
1056c----------------------------------------------------------------------
1057c Determines parameters of the funcion func
1058c representing the best fit of the data.
1059c At the end of the run, the "best" parameters are stored in parDf(n,3).
1060c The function func has to be defined via "function" using the parameters
1061c parDf(n,3), n=1,numparDf .
1062c Parameters as well as data are stored on /fitpar/:
1063c numparDf: number of parameters (input)
1064c parDf: array containing parameters:
1065c parDf(n,1): lower limit (input)
1066c parDf(n,2): upper limit (input)
1067c parDf(n,3): current parameter (internal and output = final result)
1068c parDf(n,4): previous parameter (internal)
1069c numdataDf: number of data points (input)
1070c datafitD: array containing data:
1071c datafitD(i,1): x value (input)
1072c datafitD(i,2): y value (input)
1073c datafitD(i,3): error (input)
1074c----------------------------------------------------------------------
1075
1076 include 'epos.inc'
1077 include 'epos.incpar'
1078 double precision func,x
1079 external func
1080
1081 ! write (ifch,*) 'numparDf,numminDf',numparDf,numminDf
1082
1083
1084c initial configuration (better if one start directly with initial one)
1085
1086c do n=numminDf,numparDf
1087c parDf(n,3)=parDf(n,1)+rangen()*(parDf(n,2)-parDf(n,1))
1088c enddo
1089
1090 chi2=0.
1091 err=0.
1092 do i=1,numdataDf
1093 x=dble(datafitD(i,1))
1094 fx=real(func(x))
1095 chi2=chi2+(log(fx)-log(datafitD(i,2)))**2/datafitD(i,3)**2
1096 err=err+(log(fx)-log(datafitD(i,2)))/datafitD(i,3)**2
1097 enddo
1098 err=abs(err)/real(numdataDf)
1099
1100c metropolis iteration
1101
1102 do i=1,nmc
1103c if(mod(i,int(real(nmc)/1000.)).eq.0)then
1104 betac=betac*(1.+1./real(nmc))!1.05
1105 betae=betae*(1.+1./real(nmc))!1.05
1106c endif
1107c if(mod(i,int(real(nmc)/20.)).eq.0)write(ifch,*)i,chi2,err
1108
1109 do n=numminDf,numparDf
1110 parDf(n,4)=parDf(n,3)
1111 enddo
1112 chi2x=chi2
1113 errx=err
1114
1115 n=numminDf+int(rangen()*(numparDf-numminDf+1))
1116 n=max(n,numminDf)
1117 n=min(n,numparDf)
1118c if(mod(i,int(real(nmc)/20.)).eq.0)write(ifch,*)n
1119
1120 10 parDf(n,3)=parDf(n,1)+rangen()*(parDf(n,2)-parDf(n,1))
1121 chi2=0
1122 err=0
1123 do j=1,numdataDf
1124 x=dble(datafitD(j,1))
1125 fx=real(func(x))
1126 chi2=chi2+(log(fx)-log(datafitD(j,2)))**2/datafitD(j,3)**2
1127 err=err+(log(fx)-log(datafitD(j,2)))/datafitD(j,3)**2
1128 enddo
1129 err=abs(err)/real(numdataDf)
1130
1131 if(chi2.gt.chi2x.and.rangen()
1132 $ .gt.exp(-min(50.,max(-50.,betac*(chi2-chi2x))))
1133 & .or.err.gt.errx.and.rangen()
1134 $ .gt.exp(-min(50.,max(-50.,betae*(err-errx))))
1135 & ) then
1136 do n=numminDf,numparDf
1137 parDf(n,3)=parDf(n,4)
1138 enddo
1139 chi2=chi2x
1140 err=errx
1141 endif
1142
1143 enddo
1144
1145 return
1146 end
1147
1148
1149c----------------------------------------------------------------------
1150
1151 SUBROUTINE fit(x,y,ndata,sig,mwt,a,b)
1152
1153c----------------------------------------------------------------------
1154c Given a set of data points x(1:ndata),y(1:ndata) with individual standard
1155c deviations sig(1:ndata), fit them to a straight line y = a + bx by
1156c minimizing chi2 .
1157c Returned are a,b and their respective probable uncertainties siga and sigb,
1158