COSMOS v7.655  COSMOSv7655
(AirShowerMC)
mnparm.f
Go to the documentation of this file.
1 *
2 * $Id: mnparm.F,v 1.2 1996/03/15 18:02:50 james Exp $
3 *
4 * $Log: mnparm.F,v $
5 * Revision 1.2 1996/03/15 18:02:50 james
6 * Modified Files:
7 * mnderi.F eliminate possible division by zero
8 * mnexcm.F suppress print on STOP when print flag=-1
9 * set FVAL3 to flag if FCN already called with IFLAG=3
10 * mninit.F set version 96.03
11 * mnlims.F remove arguments, not needed
12 * mnmigr.F VLEN -> LENV in debug print statement
13 * mnparm.F move call to MNRSET to after NPAR redefined, to zero all
14 * mnpsdf.F eliminate possible division by zero
15 * mnscan.F suppress printout when print flag =-1
16 * mnset.F remove arguments in call to MNLIMS
17 * mnsimp.F fix CSTATU so status is PROGRESS only if new minimum
18 * mnvert.F eliminate possible division by zero
19 *
20 * Revision 1.1.1.1 1996/03/07 14:31:31 mclareni
21 * Minuit
22 *
23 *
24  SUBROUTINE mnparm(K,CNAMJ,UK,WK,A,B,IERFLG)
25 *
26 * $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
27 *
28 * $Log: d506dp.inc,v $
29 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
30 * Minuit
31 *
32 *
33 *
34 *
35 * d506dp.inc
36 *
37 C ************ DOUBLE PRECISION VERSION *************
38  IMPLICIT DOUBLE PRECISION (a-h,o-z)
39 CC Called from MNPARS and user-callable
40 CC Implements one parameter definition, that is:
41 CC K (external) parameter number
42 CC CNAMK parameter name
43 CC UK starting value
44 CC WK starting step size or uncertainty
45 CC A, B lower and upper physical parameter limits
46 CC and sets up (updates) the parameter lists.
47 CC Output: IERFLG=0 if no problems
48 CC >0 if MNPARM unable to implement definition
49 CC
50 *
51 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
52 *
53 * $Log: d506cm.inc,v $
54 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
55 * Minuit
56 *
57 *
58 *
59 *
60 * d506cm.inc
61 *
62  parameter(mne=100 , mni=50)
63  parameter(mnihl=mni*(mni+1)/2)
64  CHARACTER*10 CPNAM
65  COMMON
66  1/mn7nam/ cpnam(mne)
67  2/mn7ext/ u(mne) ,alim(mne) ,blim(mne)
68  3/mn7err/ erp(mni) ,ern(mni) ,werr(mni) ,globcc(mni)
69  4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
70  5/mn7int/ x(mni) ,xt(mni) ,dirin(mni)
71  6/mn7fx2/ xs(mni) ,xts(mni) ,dirins(mni)
72  7/mn7der/ grd(mni) ,g2(mni) ,gstep(mni) ,gin(mne) ,dgrd(mni)
73  8/mn7fx3/ grds(mni) ,g2s(mni) ,gsteps(mni)
74  9/mn7fx1/ ipfix(mni) ,npfix
75  a/mn7var/ vhmat(mnihl)
76  b/mn7vat/ vthmat(mnihl)
77  c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
78 C
79  parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
80  parameter(zero=0.0, one=1.0, half=0.5)
81  COMMON
82  d/mn7npr/ maxint ,npar ,maxext ,nu
83  e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
84  e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
85  f/mn7tit/ cfrom ,cstatu ,ctitl ,cword ,cundef ,cvrsn ,covmes
86  g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
87  h/mn7min/ amin ,up ,edm ,fval3 ,epsi ,apsi ,dcovar
88  i/mn7cnv/ nfcn ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
89  j/mn7arg/ word7(maxp)
90  k/mn7log/ lwarn ,lrepor ,limset ,lnolim ,lnewmn ,lphead
91  l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
92  m/mn7rpt/ xpt(maxcpt) ,ypt(maxcpt)
93  n/mn7cpt/ chpt(maxcpt)
94  o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr ,ke2cr
95  CHARACTER CTITL*50, CWORD*(maxcwd), CUNDEF*10, CFROM*8,
96  + cvrsn*6, covmes(0:3)*22, cstatu*10, chpt*1
97  LOGICAL LWARN, LREPOR, LIMSET, LNOLIM, LNEWMN, LPHEAD
98  CHARACTER*(*) CNAMJ
99  CHARACTER CNAMK*10, CHBUFI*4
100 C
101  cnamk = cnamj
102  kint = npar
103  IF (k.LT.1 .OR. k.GT.maxext) THEN
104 C parameter number exceeds allowed maximum value
105  WRITE (isyswr,9) k,maxext
106  9 FORMAT (/' MINUIT USER ERROR. PARAMETER NUMBER IS',i11/
107  + ', ALLOWED RANGE IS ONE TO',i4/)
108  GO TO 800
109  ENDIF
110 C normal parameter request
111  ktofix = 0
112  IF (nvarl(k) .LT. 0) GO TO 50
113 C previously defined parameter is being redefined
114 C find if parameter was fixed
115  DO 40 ix= 1, npfix
116  IF (ipfix(ix) .EQ. k) ktofix = k
117  40 CONTINUE
118  IF (ktofix .GT. 0) THEN
119  CALL mnwarn('W','PARAM DEF','REDEFINING A FIXED PARAMETER.')
120  IF (kint .GE. maxint) THEN
121  WRITE (isyswr,'(A)') ' CANNOT RELEASE. MAX NPAR EXCEEDED.'
122  GO TO 800
123  ENDIF
124  CALL mnfree(-k)
125  ENDIF
126 C if redefining previously variable parameter
127  IF(niofex(k) .GT. 0) kint = npar-1
128  50 CONTINUE
129 C
130 C . . .print heading
131  IF (lphead .AND. isw(5).GE.0) THEN
132  WRITE (isyswr,61)
133  lphead = .false.
134  ENDIF
135  61 FORMAT(/' PARAMETER DEFINITIONS:'/
136  + ' NO. NAME VALUE STEP SIZE LIMITS')
137  IF (wk .GT. zero) GO TO 122
138 C . . .constant parameter . . . .
139  IF (isw(5) .GE. 0) WRITE (isyswr, 82) k,cnamk,uk
140  82 FORMAT (1x,i5,1x,1h',A10,1H',1x,g13.5, ' constant')
141  nvl = 0
142  GO TO 200
143  122 IF (a.EQ.zero .AND. b.EQ.zero) THEN
144 C variable parameter without limits
145  nvl = 1
146  IF (isw(5) .GE. 0) WRITE (isyswr, 127) k,cnamk,uk,wk
147  127 FORMAT (1x,i5,1x,1h',A10,1H',1x,2g13.5, ' no limits')
148  ELSE
149 C variable parameter with limits
150  nvl = 4
151  lnolim = .false.
152  IF (isw(5) .GE. 0) WRITE (isyswr, 132) k,cnamk,uk,wk,a,b
153  132 FORMAT(1x,i5,1x,1h',A10,1H',1x,2g13.5,2x,2g13.5)
154  ENDIF
155 C . . request for another variable parameter
156  kint = kint + 1
157  IF (kint .GT. maxint) THEN
158  WRITE (isyswr,135) maxint
159  135 FORMAT (/' MINUIT USER ERROR. TOO MANY VARIABLE PARAMETERS.'/
160  + ' THIS VERSION OF MINUIT DIMENSIONED FOR',i4//)
161  GO TO 800
162  ENDIF
163  IF (nvl .EQ. 1) GO TO 200
164  IF (a .EQ. b) THEN
165  WRITE (isyswr,'(/A,A/A/)') ' USER ERROR IN MINUIT PARAMETER',
166  + ' DEFINITION',' UPPER AND LOWER LIMITS EQUAL.'
167  GO TO 800
168  ENDIF
169  IF (b .LT. a) THEN
170  sav = b
171  b = a
172  a = sav
173  CALL mnwarn('W','PARAM DEF','PARAMETER LIMITS WERE REVERSED.')
174  IF (lwarn) lphead=.true.
175  ENDIF
176  IF ((b-a) .GT. 1.0e7) THEN
177  WRITE (chbufi,'(I4)') k
178  CALL mnwarn('W','PARAM DEF',
179  + 'LIMITS ON PARAM'//chbufi//' TOO FAR APART.')
180  IF (lwarn) lphead=.true.
181  ENDIF
182  danger = (b-uk)*(uk-a)
183  IF (danger .LT. 0.)
184  + CALL mnwarn('W','PARAM DEF','STARTING VALUE OUTSIDE LIMITS.')
185  IF (danger .EQ. 0.)
186  + CALL mnwarn('W','PARAM DEF','STARTING VALUE IS AT LIMIT.')
187  200 CONTINUE
188 C . . . input OK, set values, arrange lists,
189 C calculate step sizes GSTEP, DIRIN
190  cfrom = 'PARAMETR'
191  nfcnfr = nfcn
192  cstatu= 'NEW VALUES'
193  nu = max(nu,k)
194  cpnam(k) = cnamk
195  u(k) = uk
196  alim(k) = a
197  blim(k) = b
198  nvarl(k) = nvl
199 C K is external number of new parameter
200 C LASTIN is the number of var. params with ext. param. no.< K
201  lastin = 0
202  DO 240 ix= 1, k-1
203  IF (niofex(ix) .GT. 0) lastin=lastin+1
204  240 CONTINUE
205 C KINT is new number of variable params, NPAR is old
206  IF (kint .EQ. npar) GO TO 280
207  IF (kint .GT. npar) THEN
208 C insert new variable parameter in list
209  DO 260 in= npar,lastin+1,-1
210  ix = nexofi(in)
211  niofex(ix) = in+1
212  nexofi(in+1)= ix
213  x(in+1) = x(in)
214  xt(in+1) = xt(in)
215  dirin(in+1) = dirin(in)
216  g2(in+1) = g2(in)
217  gstep(in+1) = gstep(in)
218  260 CONTINUE
219  ELSE
220 C remove variable parameter from list
221  DO 270 in= lastin+1,kint
222  ix = nexofi(in+1)
223  niofex(ix) = in
224  nexofi(in)= ix
225  x(in)= x(in+1)
226  xt(in)= xt(in+1)
227  dirin(in)= dirin(in+1)
228  g2(in)= g2(in+1)
229  gstep(in)= gstep(in+1)
230  270 CONTINUE
231  ENDIF
232  280 CONTINUE
233  ix = k
234  niofex(ix) = 0
235  npar = kint
236  CALL mnrset(1)
237 C lists are now arranged . . . .
238  IF (nvl .GT. 0) THEN
239  in = lastin+1
240  nexofi(in) = ix
241  niofex(ix) = in
242  sav = u(ix)
243  CALL mnpint(sav,ix,pinti)
244  x(in) = pinti
245  xt(in) = x(in)
246  werr(in) = wk
247  sav2 = sav + wk
248  CALL mnpint(sav2,ix,pinti)
249  vplu = pinti - x(in)
250  sav2 = sav - wk
251  CALL mnpint(sav2,ix,pinti)
252  vminu = pinti - x(in)
253  dirin(in) = 0.5 * (abs(vplu) +abs(vminu))
254  g2(in) = 2.0*up / dirin(in)**2
255  gsmin = 8.*epsma2*abs(x(in))
256  gstep(in) = max(gsmin, 0.1*dirin(in))
257  IF (amin .NE. undefi) THEN
258  small = sqrt(epsma2*(amin+up)/up)
259  gstep(in) = max(gsmin, small*dirin(in))
260  ENDIF
261  grd(in) = g2(in)*dirin(in)
262 C if parameter has limits
263  IF (nvarl(k) .GT. 1) THEN
264  IF (gstep(in).GT. 0.5) gstep(in)=0.5
265  gstep(in) = -gstep(in)
266  ENDIF
267  ENDIF
268  IF (ktofix .GT. 0) THEN
269  kinfix = niofex(ktofix)
270  IF (kinfix .GT. 0) CALL mnfixp(kinfix,ierr)
271  IF (ierr .GT. 0) GO TO 800
272  ENDIF
273  ierflg = 0
274  RETURN
275 C error on input, unable to implement request . . . .
276  800 CONTINUE
277  ierflg = 1
278  RETURN
279  END
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
z
nodes z
Definition: cstdatmos0.fNew2.h:93
up
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data h g *is for param c g data up(2, 1)/7.0d0/
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130
i
nodes i
Definition: cstdatmos0.fNew2.h:93
mnfree
subroutine mnfree(K)
Definition: mnfree.f:10
mnrset
subroutine mnrset(IOPT)
Definition: mnrset.f:10
true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
xs
! constants thru Cosmos real ! if multiplied to deg radian Torad ! light velocity m sec ! infinty ! kg m2 *Togpcm2 g cm2 ! g cm2 *Tokgpm2 kg m2 ! cm *Tom m ! m *Tocm cm ! g cm3 *Tokgpm3 kg m3 ! kg m3 *Togpcm3 g cm3 ! sec *Tonsec nsec ! Tesla m ! Avogadro *A2deninv ! mfp *n * xs
Definition: Zglobalc.h:18
p
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
mnpint
subroutine mnpint(PEXTI, I, PINTI)
Definition: mnpint.f:10
o
struct ob o[NpMax]
Definition: Zprivate.h:34
mnparm
subroutine mnparm(K, CNAMJ, UK, WK, A, B, IERFLG)
Definition: mnparm.f:25
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
a
nodes a
Definition: cstdatmos0.fNew2.h:93
b
dE dx *! Nuc Int sampling table b
Definition: cblkMuInt.h:130
m
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
Definition: ZlatfitBD.h:35
mnfixp
subroutine mnfixp(IINT, IERR)
Definition: mnfixp.f:10
false
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
h
dE dx *! Nuc Int sampling table h
Definition: cblkMuInt.h:130
g
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
maxp
!onst int maxp
Definition: Zprivate.h:3
mnwarn
subroutine mnwarn(COPT, CORG, CMES)
Definition: mnwarn.f:10
n
integer n
Definition: Zcinippxc.h:1
softenpik::half
integer, parameter half
Definition: csoftenPiK.f:108
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
f
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130