COSMOS v7.655  COSMOSv7655
(AirShowerMC)
mnmnot.f
Go to the documentation of this file.
1 *
2 * $Id: mnmnot.F,v 1.1.1.1 1996/03/07 14:31:30 mclareni Exp $
3 *
4 * $Log: mnmnot.F,v $
5 * Revision 1.1.1.1 1996/03/07 14:31:30 mclareni
6 * Minuit
7 *
8 *
9  SUBROUTINE mnmnot(FCN,ILAX,ILAX2,VAL2PL,VAL2MI,FUTIL)
10 *
11 * $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
12 *
13 * $Log: d506dp.inc,v $
14 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
15 * Minuit
16 *
17 *
18 *
19 *
20 * d506dp.inc
21 *
22 C ************ DOUBLE PRECISION VERSION *************
23  IMPLICIT DOUBLE PRECISION (a-h,o-z)
24 CC Performs a MINOS error analysis on one parameter.
25 CC The parameter ILAX is varied, and the minimum of the
26 CC function with respect to the other parameters is followed
27 CC until it crosses the value FMIN+UP.
28 CC
29 *
30 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
31 *
32 * $Log: d506cm.inc,v $
33 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
34 * Minuit
35 *
36 *
37 *
38 *
39 * d506cm.inc
40 *
41  parameter(mne=100 , mni=50)
42  parameter(mnihl=mni*(mni+1)/2)
43  CHARACTER*10 CPNAM
44  COMMON
45  1/mn7nam/ cpnam(mne)
46  2/mn7ext/ u(mne) ,alim(mne) ,blim(mne)
47  3/mn7err/ erp(mni) ,ern(mni) ,werr(mni) ,globcc(mni)
48  4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
49  5/mn7int/ x(mni) ,xt(mni) ,dirin(mni)
50  6/mn7fx2/ xs(mni) ,xts(mni) ,dirins(mni)
51  7/mn7der/ grd(mni) ,g2(mni) ,gstep(mni) ,gin(mne) ,dgrd(mni)
52  8/mn7fx3/ grds(mni) ,g2s(mni) ,gsteps(mni)
53  9/mn7fx1/ ipfix(mni) ,npfix
54  a/mn7var/ vhmat(mnihl)
55  b/mn7vat/ vthmat(mnihl)
56  c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
57 C
58  parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
59  parameter(zero=0.0, one=1.0, half=0.5)
60  COMMON
61  d/mn7npr/ maxint ,npar ,maxext ,nu
62  e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
63  e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
64  f/mn7tit/ cfrom ,cstatu ,ctitl ,cword ,cundef ,cvrsn ,covmes
65  g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
66  h/mn7min/ amin ,up ,edm ,fval3 ,epsi ,apsi ,dcovar
67  i/mn7cnv/ nfcn ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
68  j/mn7arg/ word7(maxp)
69  k/mn7log/ lwarn ,lrepor ,limset ,lnolim ,lnewmn ,lphead
70  l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
71  m/mn7rpt/ xpt(maxcpt) ,ypt(maxcpt)
72  n/mn7cpt/ chpt(maxcpt)
73  o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr ,ke2cr
74  CHARACTER CTITL*50, CWORD*(maxcwd), CUNDEF*10, CFROM*8,
75  + cvrsn*6, covmes(0:3)*22, cstatu*10, chpt*1
76  LOGICAL LWARN, LREPOR, LIMSET, LNOLIM, LNEWMN, LPHEAD
77  EXTERNAL fcn,futil
78  dimension xdev(mni),w(mni),gcc(mni)
79  CHARACTER*4 CPOS,CNEG,CSIG
80  parameter(cpos='POSI',cneg='NEGA')
81 C . . save and prepare start vals
82  isw2 = isw(2)
83  isw4 = isw(4)
84  sigsav = edm
85  istrav = istrat
86  dc = dcovar
87  lnewmn = .false.
88  apsi = epsi*0.5
89  abest=amin
90  mpar=npar
91  nfmxin = nfcnmx
92  DO 125 i= 1, mpar
93  125 xt(i) = x(i)
94  DO 130 j= 1, mpar*(mpar+1)/2
95  130 vthmat(j) = vhmat(j)
96  DO 135 i= 1, mpar
97  gcc(i) = globcc(i)
98  135 w(i) = werr(i)
99  it = niofex(ilax)
100  erp(it) = 0.
101  ern(it) = 0.
102  CALL mninex(xt)
103  ut = u(ilax)
104  IF (nvarl(ilax) .EQ. 1) THEN
105  alim(ilax) = ut -100.*w(it)
106  blim(ilax) = ut +100.*w(it)
107  ENDIF
108  ndex = it*(it+1)/2
109  xunit = sqrt(up/vthmat(ndex))
110  marc = 0
111  DO 162 i= 1, mpar
112  IF (i .EQ. it) GO TO 162
113  marc = marc + 1
114  imax = max(it,i)
115  indx = imax*(imax-1)/2 + min(it,i)
116  xdev(marc) = xunit*vthmat(indx)
117  162 CONTINUE
118 C fix the parameter in question
119  CALL mnfixp (it,ierr)
120  IF (ierr .GT. 0) THEN
121  WRITE (isyswr,'(A,I5,A,I5)')
122  + ' MINUIT ERROR. CANNOT FIX PARAMETER',ilax,' INTERNAL',it
123  GO TO 700
124  ENDIF
125 C . . . . . Nota Bene: from here on, NPAR=MPAR-1
126 C Remember: MNFIXP squeezes IT out of X, XT, WERR, and VHMAT,
127 C not W, VTHMAT
128  DO 500 isig= 1,2
129  IF (isig .EQ. 1) THEN
130  sig = 1.0
131  csig = cpos
132  ELSE
133  sig = -1.0
134  csig = cneg
135  ENDIF
136 C . sig=sign of error being calcd
137  IF (isw(5) .GT. 1) WRITE (isyswr,806) csig,ilax,cpnam(ilax)
138  806 FORMAT (/' DETERMINATION OF ',a4,'TIVE MINOS ERROR FOR PARAMETER',
139  + i3, 2x ,a)
140  IF (isw(2).LE.0) CALL mnwarn('D','MINOS','NO COVARIANCE MATRIX.')
141  nlimit = nfcn + nfmxin
142  istrat = max(istrav-1,0)
143  du1 = w(it)
144  u(ilax) = ut + sig*du1
145  u(ilax) = min(u(ilax),blim(ilax))
146  u(ilax) = max(u(ilax),alim(ilax))
147  delu = u(ilax) - ut
148 C stop if already at limit with negligible step size
149  IF (abs(delu)/(abs(ut)+abs(u(ilax))) .LT. epsmac) GO TO 440
150  fac = delu/w(it)
151  DO 185 i= 1, npar
152  185 x(i) = xt(i) + fac*xdev(i)
153  IF (isw(5) .GT. 1) WRITE (isyswr,801) ilax,ut,delu,u(ilax)
154  801 FORMAT (/' PARAMETER',i4,' SET TO',e11.3,' + ',e10.3,' = ',e12.3)
155 C loop to hit AMIN+UP
156  ke1cr = ilax
157  ke2cr = 0
158  xmidcr = u(ilax)
159  xdircr = delu
160 C
161  amin = abest
162  nfcnmx = nlimit - nfcn
163  CALL mncros(fcn,aopt,iercr,futil)
164  IF (abest-amin .GT. 0.01*up) GO TO 650
165  IF (iercr .EQ. 1) GO TO 440
166  IF (iercr .EQ. 2) GO TO 450
167  IF (iercr .EQ. 3) GO TO 460
168 C . error successfully calculated
169  eros = xmidcr-ut + aopt*xdircr
170  IF (isw(5) .GT. 1) WRITE (isyswr,808) csig,ilax,cpnam(ilax),eros
171  808 FORMAT (/9x,4hthe ,a4, 29htive minos error of PARAMETER,i3, 2h
172  +, ,a10, 4h, is ,e12.4)
173  GO TO 480
174 C . . . . . . . . failure returns
175  440 IF (isw(5) .GE. 1) WRITE(isyswr,807) csig,ilax,cpnam(ilax)
176  807 FORMAT (5x,'THE ',a4,'TIVE MINOS ERROR OF PARAMETER',i3,', ',a,
177  +', EXCEEDS ITS LIMIT.'/)
178  eros = undefi
179  GO TO 480
180  450 IF (isw(5) .GE. 1) WRITE (isyswr, 802) csig,ilax,nfmxin
181  802 FORMAT (9x,'THE ',a,'TIVE MINOS ERROR',i4,' REQUIRES MORE THAN',
182  + i5,' FUNCTION CALLS.'/)
183  eros = 0.
184  GO TO 480
185  460 IF (isw(5) .GE. 1) WRITE (isyswr, 805) csig,ilax
186  805 FORMAT (25x,a,'TIVE MINOS ERROR NOT CALCULATED FOR PARAMETER',i4/)
187  eros = 0.
188 C
189  480 IF (isw(5) .GT. 1) WRITE (isyswr,'(5X, 74(1H*))')
190  IF (sig .LT. zero) THEN
191  ern(it) = eros
192  IF (ilax2.GT.0 .AND. ilax2.LE.nu) val2mi = u(ilax2)
193  ELSE
194  erp(it) = eros
195  IF (ilax2.GT.0 .AND. ilax2.LE.nu) val2pl = u(ilax2)
196  ENDIF
197  500 CONTINUE
198 C . . parameter finished. reset v
199 C normal termination
200  itaur = 1
201  CALL mnfree(1)
202  DO 550 j= 1, mpar*(mpar+1)/2
203  550 vhmat(j) = vthmat(j)
204  DO 595 i= 1, mpar
205  werr(i) = w(i)
206  globcc(i) = gcc(i)
207  595 x(i) = xt(i)
208  CALL mninex (x)
209  edm = sigsav
210  amin = abest
211  isw(2) = isw2
212  isw(4) = isw4
213  dcovar = dc
214  GO TO 700
215 C new minimum
216  650 lnewmn = .true.
217  isw(2) = 0
218  dcovar = 1.
219  isw(4) = 0
220  sav = u(ilax)
221  itaur = 1
222  CALL mnfree(1)
223  u(ilax) = sav
224  CALL mnexin(x)
225  edm = bigedm
226 C in any case
227  700 CONTINUE
228  itaur = 0
229  nfcnmx = nfmxin
230  istrat = istrav
231  RETURN
232  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
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
mnmnot
subroutine mnmnot(FCN, ILAX, ILAX2, VAL2PL, VAL2MI, FUTIL)
Definition: mnmnot.f:10
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
o
struct ob o[NpMax]
Definition: Zprivate.h:34
mninex
subroutine mninex(PINT)
Definition: mninex.f:10
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
mncros
subroutine mncros(FCN, AOPT, IERCR, FUTIL)
Definition: mncros.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
e10
*************************block data cblkTracking *************************implicit none data *ExactThick *Freec *RatioToE0 *MagChgDist *TimeStructure *Truncn *Truncx data *IncMuonPolari *KEminObs *ThinSampling *EthinRatio *Generate *LpmEffect *MagPairEmin e10
Definition: cblkTracking.h:9
maxp
!onst int maxp
Definition: Zprivate.h:3
mnexin
subroutine mnexin(PINT)
Definition: mnexin.f:10
mnwarn
subroutine mnwarn(COPT, CORG, CMES)
Definition: mnwarn.f:10
n
integer n
Definition: Zcinippxc.h:1
is
block data cblkIncident data *Za1ry is
Definition: cblkIncident.h:5
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