simc_semi/simc.f - annotate

Return to simc.f CVS log

Up to [HallC] / simc_semi

1 gaskelld 1.1 program simc 2 3 ! Last modified: 4 ! 5 ! This program represents variations on themes by 6 ! N. C. R. Makins (Argonne National Lab), 7 ! T. G. O'Neill (Argonne National Lab), and 8 ! Seemingly Countless Others (Virtually Everywhere). 9 ! 10 implicit none 11 include 'simulate_init.inc' 12 include 'histograms_init.inc' 13 include 'radc.inc' 14 include 'hbook.inc' 15 include 'sos/struct_sos.inc' 16 include 'hms/struct_hms.inc' 17 include 'hrsr/struct_hrsr.inc' 18 include 'hrsl/struct_hrsl.inc' 19 include 'shms/struct_shms.inc' 20 21 integer4 i, ninform 22 gaskelld 1.1 real8 r, sum_sigcc, tmpnum 23 real8 domega_e, domega_p !populated e/hadron solid angles. 24 logical success 25 logical pass_cuts 26 character filename80, genfile80, histfile80, timestring123 27 character timestring223,genifile80 28 record /event/ vertex, vertex0, orig, recon 29 record /event_main/ main 30 record /contrib/ contrib 31 record /histograms/ H 32 record /event_central/ central 33 record /sums_twoarm/ sumerr, sumerr2, aveerr, resol 34 35 real8 one 36 parameter (one=1.0d0) !double precision 1 for subroutine calls 37 38 real*8 grnd
39 gaskelld 1.4 real*8 ang_targ_earm,ang_targ_parm
40 gaskelld 1.1 41 ! INITIALIZE 42 !----------- 43 44 ! ... initialize histogram area for HBOOK 45 46 call hlimit(PawSize) 47 48 ! ... read in the data base 49 50 call dbase_read(H) 51 52 if (debug(3)) write(6,) 'Main after dbrd: p,e th', 53 > spec.p.theta,spec.e.theta,using_P_arm_montecarlo,using_E_arm_montecarlo 54 55 ! ... open diagnostic ntuple file, if requested 56 57 i = index(base,' ') 58 if (Nntu.gt.0) then 59 filename = 'worksim/'//base(1:i-1)//'.rzdat' 60 call NtupleInit(filename) 61 gaskelld 1.1 endif 62 63 ! ... set up some quantities that the radiative corrections routines will need 64 ! ... N.B. Call this even if we're not using radiative corrections -- the 65 ! ... target thicknesses seen by the incoming and 66 ! ... scattered particles (in radiation lengths) are determined here, and 67 ! ... are needed for the multiple scattering calculations 68 ! ... done in the Monte Carlos. 69 70 if (debug(2)) write(6,)'sim: calling radc_init' 71 call radc_init 72 if (debug(2)) write(6,)'sim: done with radc_init' 73 74 ! ... compute some quantities for a central event 75 76 call calculate_central(central,vertex0) 77 if (debug(2)) write(6,)'sim: done with calc_central' 78 if (debug(2)) write(6,)'central.sigcc=',central.sigcc 79 if (debug(4)) write(6,)'sim: at 1' 80 81 targetfac=targ.mass_amu/3.75914d+6/(targ.abundancy/100.) 82 gaskelld 1.1 > abs(cos(targ.angle))/(targ.thick1000.) 83 if (debug(4)) write(6,)'sim: at 2' 84 85 ! ... Note: EXPER.charge is in mC and the luminosity comes out in fm^-2 86 ! ... now luminosity is in ub^-1 87 88 luminosity = EXPER.charge/targetfac 89 if (debug(4)) write(6,)'sim: at 3' 90 91 ! ... initialize the random number generator and number of attempts 92 93 r = grnd() 94 ntried = 0 95
96 gaskelld 1.4 ! GAW - insert calls to initialize target field for both arms 97 ! mkj 10-20-2003 add if statements to determine the angle magnitude 98 ! and sign between target direction and e-arm or p-arm. 99 !
100 gaskelld 1.6 if(using_tgt_field) then 101 if ( degradabs(targ_Bphi-spec.e.phi) < .01) then 102 if (targ_Bangle .ge. spec.e.theta) then 103 ang_targ_earm = -1sin(spec.e.phi)(targ_Bangle-spec.e.theta) 104 else 105 ang_targ_earm = +1sin(spec.e.phi)(spec.e.theta-targ_Bangle) 106 endif 107 elseif ( degradabs(targ_Bphi-spec.e.phi)-180.0 < .01) then 108 ang_targ_earm = +1sin(spec.e.phi)(targ_Bangle+spec.e.theta) 109 else 110 write(,) ' Error determining angle between target and e-arm' 111 write(,) ' targ_Bangle =',targ_Bangledegrad,' targ_Bphi =',targ_Bphidegrad 112 write(,) ' Central e-arm angle =',spec.e.thetadegrad,' e-arm phi =',spec.e.phidegrad 113 endif 114 c 115 if ( degradabs(targ_Bphi-spec.p.phi) < .01) then 116 if (targ_Bangle .ge. spec.p.theta) then 117 ang_targ_parm = -1sin(spec.p.phi)(targ_Bangle-spec.p.theta) 118 else 119 ang_targ_parm = +1sin(spec.p.phi)(targ_Bangle-spec.p.theta) 120 endif 121 gaskelld 1.6 elseif ( degradabs(targ_Bphi-spec.p.phi)-180.0 < .01) then 122 ang_targ_parm = +1sin(spec.p.phi)(targ_Bangle+spec.p.theta) 123 else 124 write(,) ' Error determining angle between target and p-arm' 125 write(,) ' targ_Bangle =',targ_Bangledegrad,' targ_Bphi =',targ_Bphidegrad 126 write(,) ' Central p-arm angle =',spec.p.thetadegrad,' p-arm phi =',spec.p.thetadegrad 127 endif
128 gaskelld 1.4 c 129 write(,) ' targ_Bangle =',targ_Bangledegrad,' targ_Bphi =',targ_Bphidegrad
130 gaskelld 1.6 write(,) ' Angle between target and e-arm',ang_targ_earmdegrad 131 write(,) ' Angle between target and p-arm',ang_targ_parmdegrad
132 gaskelld 1.4 c 133
134 gaskelld 1.6 call trgInit(tgt_field_file,ang_targ_earmdegrad,0., 135 1 ang_targ_parmdegrad,0.) 136 endif
137 gaskelld 1.4
138 gaskelld 1.1 ! LOOP OVER SIMULATED EVENTS 139 !--------------------------- 140 write(6,) ' ' 141 call date (timestring1) 142 call time (timestring1(11:23)) 143 nevent = 0 144 ninform = 1000 145 sum_sigcc = 0.0 !sum of main.sigcc (for generating ave sigcc) 146 147 do while (nevent.lt.abs(ngen)) 148 149 ! reset decay distance, initial hadron mass (modified if particle decays), 150 ! and some ntuple variables. 151 152 Mh2_final = Mh2 153 ntup.radphot = 0.0 154 ntup.radarm = 0.0 155 decdist = 0.0 !decay distance. 156 ntup.resfac = 0.0 157 ntup.sigcm = 0.0 158 ntup.krel = 0.0 159 gaskelld 1.1 ntup.mm = 0.0 160 ntup.mmA = 0.0 161 ntup.t = 0.0 162 163 ! Setup for this event 164 165 ! ... keep the human interested 166 167 ntried = ntried + 1 168 169 if(debug(2)) then 170 write(6,)'*******************************************' 171 write(6,)'SIM: ntried =',ntried 172 write(6,)' ncontribute =',ncontribute 173 endif 174 175 if (mod(ntried,ninform).eq.1) then 176 write (6,'(1x,a,i9,a,i8,a,e11.4)') 'Generating Event', 177 > ntried, ' ... ', nevent,' successes so far - Monitor:', 178 > wtcontributeluminosity/ntried 179 if (ntried.ge.5000) ninform = 20000 180 gaskelld 1.1 endif 181 182 ! ... generate an event 183 184 call generate(main,vertex,orig,success) 185 if(debug(2)) write(6,)'sim: after gen, success =',success 186 187 ! Run the event through various manipulations, checking to see whether 188 ! it made it at each stage 189 190 ! ... run through spectrometers and check SP cuts 191 192 if(debug(3)) write(6,)'sim: before mc: orig.p.E =',orig.p.E 193 if(success)call montecarlo(orig,main,recon,success) 194 if(debug(2)) write(6,*)'sim: after mc, success =',success
195 gaskelld 1.2
196 gaskelld 1.1 ! ... calculate everything else about the event 197 198 if (success) then 199 call complete_recon_ev(recon,success) 200 endif
201 gaskelld 1.2
202 gaskelld 1.1 if(debug(2)) write(6,)'sim: after comp_ev, success =',success 203 if(debug(5)) write(6,) 'recon.Em,recon.Pm',recon.Em,recon.Pm 204 ! ... calculate remaining pieces of the main structure 205 206 if (success) call complete_main(.false.,main,vertex,vertex0,recon,success) 207 ! ... Apply SPedge cuts to success if hard_cuts is set. 208 209 pass_cuts = .not. ( 210 > recon.e.delta .le. (SPedge.e.delta.min+slop.MC.e.delta.used) .or. 211 > recon.e.delta .ge. (SPedge.e.delta.max-slop.MC.e.delta.used) .or. 212 > recon.e.yptar .le. (SPedge.e.yptar.min+slop.MC.e.yptar.used) .or. 213 > recon.e.yptar .ge. (SPedge.e.yptar.max-slop.MC.e.yptar.used) .or. 214 > recon.e.xptar .le. (SPedge.e.xptar.min+slop.MC.e.xptar.used) .or. 215 > recon.e.xptar .ge. (SPedge.e.xptar.max-slop.MC.e.xptar.used) .or. 216 > recon.p.delta .le. (SPedge.p.delta.min+slop.MC.p.delta.used) .or. 217 > recon.p.delta .ge. (SPedge.e.delta.max-slop.MC.p.delta.used) .or. 218 > recon.p.yptar .le. (SPedge.p.yptar.min+slop.MC.p.yptar.used) .or. 219 > recon.p.yptar .ge. (SPedge.p.yptar.max-slop.MC.p.yptar.used) .or. 220 > recon.p.xptar .le. (SPedge.p.xptar.min+slop.MC.p.xptar.used) .or. 221 > recon.p.xptar .ge. (SPedge.p.xptar.max-slop.MC.p.xptar.used) ) 222 223 gaskelld 1.1 if (hard_cuts) then !apply spec. limit and Em cuts. 224 if (.not.pass_cuts) success = .false. 225 if (doing_eep .and. (recon.Em .gt. cuts.Em.max)) success = .false. 226 endif 227 228 if (success) sum_sigcc = sum_sigcc + main.sigcc 229 230 ! Em and Pm histos are NEW. Check that they don't suffer from energy 231 ! shifts (e.g. coulomb distortions - see update_range call in limits_update). 232 233 if(ntried.gt.0)then 234 call inc(H.geni.e.delta, vertex.e.delta, one) 235 call inc(H.geni.e.yptar, vertex.e.yptar, one) 236 call inc(H.geni.e.xptar, -vertex.e.xptar, one) 237 call inc(H.geni.p.delta, vertex.p.delta, one) 238 call inc(H.geni.p.yptar, vertex.p.yptar, one) 239 call inc(H.geni.p.xptar, -vertex.p.xptar, one) 240 call inc(H.geni.Em, vertex.Em, one) 241 call inc(H.geni.Pm, vertex.Pm, one) 242 endif 243 244 gaskelld 1.1 if (success) then 245 if (debug(2)) write(6,)'sim: We have a success!' 246 247 ! ... Output ntuple and histograms if passed all cuts, or if not using 248 ! ... SPedge limits as hard cuts. 249 250 ! ... increment the "experimental" spectrometer arm and "contribution" 251 ! ... histograms 252 call inc(H.RECON.e.delta, main.RECON.e.delta, main.weight) 253 call inc(H.RECON.e.yptar, main.RECON.e.yptar, main.weight) 254 call inc(H.RECON.e.xptar, main.RECON.e.xptar, main.weight) 255 call inc(H.RECON.p.delta, main.RECON.p.delta, main.weight) 256 call inc(H.RECON.p.yptar, main.RECON.p.yptar, main.weight) 257 call inc(H.RECON.p.xptar, main.RECON.p.xptar, main.weight) 258 call inc(H.RECON.Em, recon.Em, one) 259 call inc(H.RECON.Pm, recon.Pm, one) 260 call inc(H.gen.e.delta, vertex.e.delta, one) 261 call inc(H.gen.e.yptar, vertex.e.yptar, one) 262 call inc(H.gen.e.xptar, -vertex.e.xptar, one) 263 call inc(H.gen.p.delta, vertex.p.delta, one) 264 call inc(H.gen.p.yptar, vertex.p.yptar, one) 265 gaskelld 1.1 call inc(H.gen.p.xptar, -vertex.p.xptar, one) 266 call inc(H.gen.Em, vertex.Em, one) 267 268 ! ... update counters and integrated weights AND keep track of resolution 269 ! ... FOR EVENTS INSIDE OF spectrometer population limits (events are only 270 ! ... generated to fill region inside of the given delta limits) and below 271 ! ... cuts.Em.max (used for elastic/quasielastic only). 272 273 if (debug(4)) write(6,)'sim: cut' 274 275 ncontribute = ncontribute + 1 276 if (debug(4)) write(6,)'sim: ncontribute =',ncontribute 277 if (.not.rad_proton_this_ev) ncontribute_no_rad_proton = 278 > ncontribute_no_rad_proton + 1 279 280 281 ! Using main.SP.e.delta means that we correct for eloss/coulomb, which we 282 ! shouldn't do. Using vertex.e.delta means radiation counts as error in recon. 283 ! --> USE MAIN.SP AND UNDERESTIMATE RESOLUTION (MISSING TARGET ELOSS) 284 285 ! For xptar/yptar, main.SP.e.xptar/yptar corrects for mult. scatt. --> USE VERTEX. 286 gaskelld 1.1 287 ! For z, we HAVE to use main.SP.e.z, but there are not corrections, so USE MAIN.SP 288 ! vertex.e.z is not defined (main.target.x/y/z --> y_E_arm 289 ! (offsets and rotation) --> spectrometer M.C. --> RECON.e.z 290 291 if (pass_cuts) then 292 npasscuts = npasscuts + 1 293 wtcontribute = wtcontribute + main.weight 294 sumerr.e.delta = sumerr.e.delta + (recon.e.delta - main.SP.e.delta) 295 sumerr.e.xptar = sumerr.e.xptar + (recon.e.xptar - vertex.e.xptar) 296 sumerr.e.yptar = sumerr.e.yptar + (recon.e.yptar - vertex.e.yptar) 297 sumerr.e.ytar = sumerr.e.ytar + (recon.e.z - main.SP.e.z) 298 sumerr.p.delta = sumerr.p.delta + (recon.p.delta - main.SP.p.delta) 299 sumerr.p.xptar = sumerr.p.xptar + (recon.p.xptar - vertex.p.xptar) 300 sumerr.p.yptar = sumerr.p.yptar + (recon.p.yptar - vertex.p.yptar) 301 sumerr.p.ytar = sumerr.p.ytar + (recon.p.z - main.SP.p.z) 302 sumerr2.e.delta= sumerr2.e.delta + (recon.e.delta - main.SP.e.delta)2 303 sumerr2.e.xptar= sumerr2.e.xptar + (recon.e.xptar - vertex.e.xptar)2 304 sumerr2.e.yptar= sumerr2.e.yptar + (recon.e.yptar - vertex.e.yptar)2 305 sumerr2.e.ytar = sumerr2.e.ytar + (recon.e.z - main.SP.e.z)2 306 sumerr2.p.delta= sumerr2.p.delta + (recon.p.delta - main.SP.p.delta)2 307 gaskelld 1.1 sumerr2.p.xptar= sumerr2.p.xptar + (recon.p.xptar - vertex.p.xptar)2 308 sumerr2.p.yptar= sumerr2.p.yptar + (recon.p.yptar - vertex.p.yptar)2 309 sumerr2.p.ytar = sumerr2.p.ytar + (recon.p.z - main.SP.p.z)2 310 endif 311 312 ! ... update the "contribution" and "slop" limits 313 314 call limits_update(main,vertex,orig,recon,doing_deuterium, 315 > doing_pion,doing_kaon,doing_delta,doing_rho,contrib,slop) 316 317 318 endif ! <success> 319 320 ! ... write out line to diagnostic ntuple file, if requested 321 322 if (Nntu.gt.0) call results_ntu_write(main,vertex,orig,recon,success) 323 324 ! END of LOOP over events 325 ! Cute thing here, if ngen < 0, then just make \|ngen\| ATTEMPTS rather than 326 ! insisting on ngen good events ... this is suitable for tests with new 327 ! and uncertain parameters, you don't want the program fishing around all 328 gaskelld 1.1 ! day for an event! 329 330 if (ngen.lt.0) then 331 nevent = nevent+1 332 else 333 if (success) nevent = nevent+1 334 endif 335 336 enddo ! <loop over ntried> 337 338 339 ! END GAME: NORMALIZE, GENERATE OUTPUT 340 !------------------------------------- 341 342 ! Our last event announcement ... and how long did it all take? 343 write (6,'(1x,''---> Last Event '',i9,'' ...'',i9,'' successes'')') ntried, nevent 344 call date (timestring2) 345 call time (timestring2(11:23)) 346 347 ! NORMALIZE! 348 349 gaskelld 1.1 ! ... put in the luminosity and efficiency factors 350 351 normfac = luminosity/ntriednevent 352 353 ! ... multiply in the relevant phase spaces (see event.f for description 354 ! ... of generated variables. Electron angles for all cases. 355 ! ... add hadron angles and electron energy for all but H(e,e'p). 356 ! ... add hadron energy for A(e,e'p) (what about phase_space?) 357 ! ... Cross sections are all in microbarn/MeVi/srj (i,j depend on reaction) 358 359 domega_e=(gen.e.yptar.max-gen.e.yptar.min)*(gen.e.xptar.max-gen.e.xptar.min)
360 gaskelld 1.2 if(doing_rho) then 361 domega_p = 4.pi 362 else 363 domega_p=(gen.p.yptar.max-gen.p.yptar.min)(gen.p.xptar.max-gen.p.xptar.min) 364 endif
365 gaskelld 1.1 366 genvol = domega_e 367 ! genvol_inclusive = genvol !may want dOmega, or dEdOmega 368 369 ! ... 2-fold to 5-fold. 370 if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon 371 1 .or.doing_delta.or.doing_rho .or. doing_semi) then 372 genvol = genvol domega_p * (gen.e.E.max-gen.e.E.min) 373 endif 374 375 if (doing_heavy.or.doing_semi) then !6-fold 376 genvol = genvol * (gen.p.E.max-gen.p.E.min) 377 endif 378 379 normfac = normfac * genvol 380 if (doing_phsp) normfac = 1.0 381 wtcontribute = wtcontributenormfac 382 383 ! Close diagnostic ntuple file, if used 384 385 if (Nntu.gt.0) call NtupleClose(filename) 386 gaskelld 1.1 387 ! Calculate resolutions 388 389 if (npasscuts.gt.1) then 390 tmpnum = dble(npasscuts) 391 aveerr.e.delta = sumerr.e.delta/tmpnum 392 aveerr.e.xptar = sumerr.e.xptar/tmpnum 393 aveerr.e.yptar = sumerr.e.yptar/tmpnum 394 aveerr.e.ytar = sumerr.e.ytar /tmpnum 395 aveerr.p.delta = sumerr.p.delta/tmpnum 396 aveerr.p.xptar = sumerr.p.xptar/tmpnum 397 aveerr.p.yptar = sumerr.p.yptar/tmpnum 398 aveerr.p.ytar = sumerr.p.ytar /tmpnum 399 resol.e.delta = sqrt(max(0.,(sumerr2.e.delta/tmpnum) - 400 > (sumerr.e.delta/tmpnum)2 )) 401 resol.e.xptar = sqrt(max(0.,(sumerr2.e.xptar/tmpnum) - 402 > (sumerr.e.xptar/tmpnum)2 )) 403 resol.e.yptar = sqrt(max(0.,(sumerr2.e.yptar/tmpnum) - 404 > (sumerr.e.yptar/tmpnum)2 )) 405 resol.e.ytar = sqrt(max(0.,(sumerr2.e.ytar/tmpnum) - 406 > (sumerr.e.ytar/tmpnum)2 )) 407 gaskelld 1.1 resol.p.delta = sqrt(max(0.,(sumerr2.p.delta/tmpnum) - 408 > (sumerr.p.delta/tmpnum)2 )) 409 resol.p.xptar = sqrt(max(0.,(sumerr2.p.xptar/tmpnum) - 410 > (sumerr.p.xptar/tmpnum)2 )) 411 resol.p.yptar = sqrt(max(0.,(sumerr2.p.yptar/tmpnum) - 412 > (sumerr.p.yptar/tmpnum)2 )) 413 resol.p.ytar = sqrt(max(0.,(sumerr2.p.ytar/tmpnum) - 414 > (sumerr.p.ytar/tmpnum)2 )) 415 endif 416 417 ! write(6,9910) 'e- delta=',10.aveerr.e.delta,10.resol.e.delta,'x10^-3' 418 ! write(6,9910) 'e- xptar=',1000.aveerr.e.xptar,1000.resol.e.xptar,'mr' 419 ! write(6,9910) 'e- yptar=',1000.aveerr.e.yptar,1000.resol.e.yptar,'mr' 420 ! write(6,9910) 'e- ytar =',10.aveerr.e.ytar,10.resol.e.ytar ,'mm' 421 ! write(6,9910) 'p delta=',10.aveerr.p.delta,10.resol.p.delta,'x10-3' 422 ! write(6,9910) 'p xptar=',1000.aveerr.p.xptar,1000.resol.p.xptar,'mr' 423 ! write(6,9910) 'p yptar=',1000.aveerr.p.yptar,1000.resol.p.yptar,'mr' 424 ! write(6,9910) 'p ytar =',10.aveerr.p.ytar,10.resol.p.ytar,'mm' 425 !9910 format(1x,a10,2f10.3,3x,a) 426 427 ! Produce output files 428 gaskelld 1.1 429 900 if (ngen.eq.0) goto 1000 430 431 ! ... the diagnostic histograms 432 433 i = index(base,' ') 434 genfile = ' ' 435 write(genfile,'(a,''.gen'')') 'outfiles/'//base(1:i-1) 436 genifile = ' ' 437 write(genifile,'(a,''.geni'')') 'outfiles/'//base(1:i-1) 438 histfile = ' ' 439 write(histfile,'(a,''.hist'')') 'outfiles/'//base(1:i-1) 440 write(6,'(1x,''... writing '',a)') genfile 441 442 write(6,) 'Come back soon!' 443 444 open(unit=7, file=genifile, status='unknown') 445 if (electron_arm.eq.1 .or. hadron_arm.eq.1) then 446 write(7,) 'HMS Trials: ',hSTOP_trials 447 write(7,) 'Slit hor/vert/corners ',hSTOP_slit_hor,hSTOP_slit_vert,hSTOP_slit_oct 448 write(7,) 'Q1 entrance/mid/exit ',hSTOP_Q1_in,hSTOP_Q1_mid,hSTOP_Q1_out 449 gaskelld 1.1 write(7,) 'Q2 entrance/mid/exit ',hSTOP_Q2_in,hSTOP_Q2_mid,hSTOP_Q2_out 450 write(7,) 'Q3 entrance/mid/exit ',hSTOP_Q3_in,hSTOP_Q3_mid,hSTOP_Q3_out 451 write(7,) 'Dipole entrance/exit ',hSTOP_D1_in,hSTOP_D1_out 452 write(7,) 'Events reaching hut ',hSTOP_hut 453 write(7,) 'DC1, DC2, Scin, Cal ',hSTOP_dc1,hSTOP_dc2,hSTOP_scin,hSTOP_cal 454 write(7,) 'Successes ',hSTOP_successes 455 write(7,) 456 endif 457 if (electron_arm.eq.2 .or. hadron_arm.eq.2) then 458 write(7,) 'SOS Trials: ',sSTOP_trials 459 write(7,) 'Slit hor/vert/corners ',sSTOP_slit_vert,sSTOP_slit_hor,sSTOP_slit_oct 460 write(7,) 'Quad entrance/mid/exit',sSTOP_quad_in,sSTOP_quad_mid,sSTOP_quad_out 461 write(7,) 'D1 entrance/exit ',sSTOP_bm01_in,sSTOP_bm01_out 462 write(7,) 'D2 entrance/exit ',sSTOP_bm02_in,sSTOP_bm02_out 463 write(7,) 'Vacuum exit ',sSTOP_exit 464 write(7,) 'Events reaching hut ',sSTOP_hut 465 write(7,) 'DC1, DC2, Scin, Cal ',sSTOP_dc1,sSTOP_dc2,sSTOP_scin,sSTOP_cal 466 write(7,) 'Successes ',sSTOP_successes 467 write(7,) 468 endif 469 if (electron_arm.eq.3 .or. hadron_arm.eq.3) then 470 gaskelld 1.1 write(7,) 'HRSr Trials: ',rSTOP_trials 471 write(7,) 'Slit hor/vert ',rSTOP_slit_vert,rSTOP_slit_hor 472 write(7,) 'Q1 entrance/mid/exit ',rSTOP_Q1_in,rSTOP_Q1_mid,rSTOP_Q1_out 473 write(7,) 'Q2 entrance/mid/exit ',rSTOP_Q2_in,rSTOP_Q2_mid,rSTOP_Q2_out 474 write(7,) 'Dipole entrance/exit ',rSTOP_D1_in,rSTOP_D1_out 475 write(7,) 'Q3 entrance/mid/exit ',rSTOP_Q3_in,rSTOP_Q3_mid,rSTOP_Q3_out 476 write(7,) 'Events reaching hut ',rSTOP_hut 477 write(7,) 'VDC1, VDC2 ',rSTOP_dc1,rSTOP_dc2 478 write(7,) 'S1, S2, Cal ',rSTOP_s1,rSTOP_s2,rSTOP_cal 479 write(7,) 480 endif 481 if (electron_arm.eq.4 .or. hadron_arm.eq.4) then 482 write(7,) 'HRSl Trials: ',lSTOP_trials 483 write(7,) 'Slit hor/vert ',lSTOP_slit_vert,lSTOP_slit_hor 484 write(7,) 'Q1 entrance/mid/exit ',lSTOP_Q1_in,lSTOP_Q1_mid,lSTOP_Q1_out 485 write(7,) 'Q2 entrance/mid/exit ',lSTOP_Q2_in,lSTOP_Q2_mid,lSTOP_Q2_out 486 write(7,) 'Dipole entrance/exit ',lSTOP_D1_in,lSTOP_D1_out 487 write(7,) 'Q3 entrance/mid/exit ',lSTOP_Q3_in,lSTOP_Q3_mid,lSTOP_Q3_out 488 write(7,) 'Events reaching hut ',lSTOP_hut 489 write(7,) 'VDC1, VDC2 ',lSTOP_dc1,lSTOP_dc2 490 write(7,) 'S1, S2, Cal ',lSTOP_s1,lSTOP_s2,lSTOP_cal 491 gaskelld 1.1 endif 492 if (electron_arm.eq.5 .or. hadron_arm.eq.5 .or. 493 > electron_arm.eq.6 .or. hadron_arm.eq.6) then 494 write(7,) 'SHMS Trials: ',shmsSTOP_trials 495 write(7,) 'Slit hor/vert/corners ',shmsSTOP_slit_hor,shmsSTOP_slit_vert,shmsSTOP_slit_oct 496 write(7,) 'Q1 entrance/mid/exit ',shmsSTOP_Q1_in,shmsSTOP_Q1_mid,shmsSTOP_Q1_out 497 write(7,) 'Q2 entrance/mid/exit ',shmsSTOP_Q2_in,shmsSTOP_Q2_mid,shmsSTOP_Q2_out 498 write(7,) 'D1 entrance/mid/exit ',shmsSTOP_D1_in,shmsSTOP_D1_mid,shmsSTOP_D1_out 499 write(7,) 'BP thingie in/out ',shmsSTOP_BP_in,shmsSTOP_BP_out 500 write(7,) 'Events reaching hut ',shmsSTOP_hut 501 write(7,) 'DC1, DC2, Scin, Cal ',shmsSTOP_dc1,shmsSTOP_dc2 502 write(7,) 'S1, S2, S3, Cal ',shmsSTOP_s1,shmsSTOP_s2,shmsSTOP_s3,shmsSTOP_cal 503 write(7,) 'Successes ',shmsSTOP_successes 504 write(7,) 505 endif 506 507 close(7) 508 509 open(unit=7, file=genfile, status='unknown') 510 511 write(7,) 'E arm Experimental Target Distributions:' 512 gaskelld 1.1 write(7,'(6a12)') 'delta','EXPERIM','yptar','EXPERIM','xptar','EXPERIM' 513 do i=1,nHbins 514 write(7,'(3(1x,2(e11.4,1x)))') 515 > H.RECON.e.delta.min+(i-0.5)H.RECON.e.delta.bin, 516 > H.RECON.e.delta.buf(i), H.RECON.e.yptar.min+(i-0.5)* 517 > H.RECON.e.yptar.bin, H.RECON.e.yptar.buf(i), 518 > H.RECON.e.xptar.min+(i-0.5)H.RECON.e.xptar.bin, 519 > H.RECON.e.xptar.buf(i) 520 enddo 521 522 write(7,) 'P arm Experimental Target Distributions:' 523 write(7,'(6a12)') 'delta','EXPERIM','yptar','EXPERIM','xptar','EXPERIM' 524 do i=1,nHbins 525 write(7,'(3(1x,2(e11.4,1x)))') 526 > H.RECON.p.delta.min+(i-0.5)H.RECON.p.delta.bin, 527 > H.RECON.p.delta.buf(i), H.RECON.p.yptar.min+(i-0.5) 528 > H.RECON.p.yptar.bin, H.RECON.p.yptar.buf(i), 529 > H.RECON.p.xptar.min+(i-0.5)H.RECON.p.xptar.bin, 530 > H.RECON.p.xptar.buf(i) 531 enddo 532 533 gaskelld 1.1 write(7,) 'Distributions of Contributing E arm Events:' 534 write(7,'(6a12)') 'delta','CONTRIB','yuptar','CONTRIB','xptar','CONTRIB' 535 do i=1,nHbins 536 write(7,'(3(1x,2(e11.4,1x)))') 537 > H.gen.e.delta.min+(i-0.5)H.gen.e.delta.bin, 538 > H.gen.e.delta.buf(i), H.gen.e.yptar.min+(i-0.5) 539 > H.gen.e.yptar.bin, H.gen.e.yptar.buf(i), 540 > H.gen.e.xptar.min+(i-0.5)H.gen.e.xptar.bin, H.gen.e.xptar.buf(i) 541 enddo 542 543 write(7,) 'Distributions of Contributing P arm Events:' 544 write(7,'(6a12)') 'delta','CONTRIB','yptar','CONTRIB','xptar','CONTRIB' 545 do i=1,nHbins 546 write(7,'(3(1x,2(e11.4,1x)))') 547 > H.gen.p.delta.min+(i-0.5)H.gen.p.delta.bin, 548 > H.gen.p.delta.buf(i), H.gen.p.yptar.min+(i-0.5) 549 > H.gen.p.yptar.bin, H.gen.p.yptar.buf(i), 550 > H.gen.p.xptar.min+(i-0.5)H.gen.p.xptar.bin, H.gen.p.xptar.buf(i) 551 enddo 552 553 write(7,) 'Original E arm Events:' 554 gaskelld 1.1 write(7,'(6a12)') 'delta','ORIGIN','yptar','ORIGIN','xptar','ORIGIN' 555 do i=1,nHbins 556 write(7,'(3(1x,2(e11.4,1x)))') 557 > H.gen.e.delta.min+(i-0.5)H.gen.e.delta.bin, 558 > H.gen.e.delta.buf(i), H.gen.e.yptar.min+(i-0.5) 559 > H.gen.e.yptar.bin, H.geni.e.yptar.buf(i), 560 > H.gen.e.xptar.min+(i-0.5)H.gen.e.xptar.bin, H.geni.e.xptar.buf(i) 561 enddo 562 563 write(7,) 'Original P arm Events:' 564 write(7,'(6a12)') 'delta','ORIGIN','yptar','ORIGIN','xptar','ORIGIN' 565 do i=1,nHbins 566 write(7,'(3(1x,2(e11.4,1x)))') 567 > H.geni.p.delta.min+(i-0.5)H.geni.p.delta.bin, 568 > H.geni.p.delta.buf(i), H.geni.p.yptar.min+(i-0.5) 569 > H.geni.p.yptar.bin, H.geni.p.yptar.buf(i), 570 > H.geni.p.xptar.min+(i-0.5)H.geni.p.xptar.bin, H.geni.p.xptar.buf(i) 571 enddo 572 573 write(7,) 'Original Em/Pm distributions:' 574 write(7,'(4a12)') 'Em','ORIGIN','Pm','ORIGIN' 575 gaskelld 1.1 do i=1,nHbins 576 write(7,'(3(1x,4(e11.4,1x)))') 577 > H.geni.Em.min+(i-0.5)H.geni.Em.bin, 578 > H.geni.Em.buf(i), H.geni.Pm.min+(i-0.5) 579 > H.geni.Pm.bin, H.geni.Pm.buf(i) 580 enddo 581 582 close(7) 583 584 ! Counters, etc report to the screen and to the diagnostic histogram file 585 ! call report(6,timestring1,timestring2,central,contrib,sum_sigcc,aveerr,resol) 586 open(unit=7,file=histfile,status='unknown') 587 call report(7,timestring1,timestring2,central,contrib,sum_sigcc,aveerr,resol) 588 close(unit=7) 589 590 1000 end 591 592 !------------------------------------------------------------------- 593 594 subroutine inc(hist,val,weight) 595 596 gaskelld 1.1 implicit none 597 include 'histograms.inc' 598 599 integer4 ibin 600 real8 val, weight 601 record /hist_entry/ hist 602 603 ibin= nint(0.5+(val-hist.min)/hist.bin) 604 if(ibin.ge.1.and.ibin.le.nHbins)then 605 hist.buf(ibin) = hist.buf(ibin) + weight 606 endif 607 608 return 609 end 610 611 !-------------------------------------------------------------------- 612 613 subroutine report(iun,timestring1,timestring2,central,contrib, 614 > sum_sigcc,aveerr,resol) 615 616 implicit none 617 gaskelld 1.1 include 'simulate.inc' 618 include 'radc.inc' 619 include 'brem.inc' 620 621 integer4 iun 622 real8 sum_sigcc 623 character23 timestring1, timestring2 624 record /contrib/ contrib 625 record /event_central/ central 626 record /sums_twoarm/ aveerr, resol 627 628 ! Output of diagnostics 629 630 ! Run time 631 write(iun,'(/1x,''BEGIN Time: '',a23)') timestring1 632 write(iun,'(1x,''END Time: '',a23)') timestring2 633 634 ! Kinematics 635 write(iun,) 'KINEMATICS:' 636 if (doing_eep) then 637 if (doing_hyd_elast) then 638 gaskelld 1.1 write(iun,) ' *-------- H(e,e''p) --------*' 639 else if (doing_deuterium) then 640 write(iun,) ' **-------- D(e,e''p) --------*' 641 else if (doing_heavy) then 642 write(iun,) ' **-------- A(e,e''p) --------*' 643 else 644 stop 'I don''t have ANY idea what (e,e''p) we''re doing!!!' 645 endif 646 else if (doing_semi) then 647 if (doing_semipi) then 648 if (targ.A .eq. 1) then 649 if(doing_hplus) then 650 write(iun,) ' **-------- H(e,e''pi+)X --------*' 651 else 652 write(iun,) ' **-------- H(e,e''pi-)X --------*' 653 endif 654 elseif (targ.A .eq. 2) then 655 if(doing_hplus) then 656 write(iun,) ' **-------- D(e,e''pi+)X --------*' 657 else 658 write(iun,) ' **-------- D(e,e''pi-)X --------*' 659 gaskelld 1.1 endif 660 else 661 stop 'I don''t have ANY idea what A(e,e''pi)X we''re doing!!!' 662 endif 663 else if (doing_semika) then 664 if (targ.A .eq. 1) then 665 if(doing_hplus) then 666 write(iun,) ' **-------- H(e,e''k+)X --------*' 667 else 668 write(iun,) ' **-------- H(e,e''k-)X --------*' 669 endif 670 elseif (targ.A .eq. 2) then 671 if(doing_hplus) then 672 write(iun,) ' **-------- D(e,e''k+)X --------*' 673 else 674 write(iun,) ' **-------- D(e,e''k-)X --------*' 675 endif 676 else 677 stop 'I don''t have ANY idea what A(e,e''k)X we''re doing!!!' 678 endif 679 else 680 gaskelld 1.1 stop 'I don''t have ANY idea what A(e,e''x)X we''re doing!!!' 681 endif 682 else if (doing_rho) then 683 if (targ.A .eq. 1) then 684 write(iun,) ' **-------- H(e,e''rho) --------*' 685 else 686 write(iun,) 'I am not set up for anything else yet!' 687 endif 688 else if (doing_delta) then 689 if (doing_hyddelta) then 690 write(6,) ' *-------- H(e,e''p)pi --------*' 691 else if (doing_deutpi) then 692 write(6,) ' **-------- D(e,e''p)pi --------*' 693 else if (doing_hepi) then 694 write(6,) ' **-------- A(e,e''p)pi --------*' 695 else 696 stop 'I don''t have ANY idea what (e,e''p)pi we''re doing!!!' 697 endif 698 else if (doing_pion) then 699 if (doing_hydpi) then 700 if (targ.A .eq. 1) then 701 gaskelld 1.1 write(iun,) ' **-------- H(e,e''pi) --------*' 702 else if (targ.A .ge.3) then 703 write(iun,) ' **-------- A(e,e''pi) --------*' 704 endif 705 else if (doing_deutpi) then 706 write(iun,) ' **-------- D(e,e''pi) --------*' 707 else if (doing_hepi) then 708 write(iun,) ' **-------- A(e,e''pi) --------*' 709 else 710 stop 'I don''t have ANY idea what (e,e''pi) we''re doing!!!' 711 endif 712 if (which_pion .eq. 0) then 713 write(iun,) ' **---- Default Final State ----*' 714 else if (which_pion .eq. 10) then 715 write(iun,) ' **---- Final State is A+pi ----*' 716 endif 717 else if (doing_kaon) then 718 if (doing_hydkaon) then 719 if (targ.A .eq. 1) then 720 write(iun,) ' **-------- H(e,e''K) --------*' 721 else if (targ.A .ge.3) then 722 gaskelld 1.1 write(iun,) ' **-------- A(e,e''K) --------*' 723 endif 724 else if (doing_deutkaon) then 725 write(iun,) ' **-------- D(e,e''K) --------*' 726 else if (doing_hekaon) then 727 write(iun,) ' **-------- A(e,e''K) --------*' 728 else 729 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!' 730 endif 731 if (which_kaon.eq.0) then 732 write(iun,) ' **---- producing a LAMBDA ----*' 733 else if (which_kaon.eq.1) then 734 write(iun,) ' **---- producing a SIGMA0 ----*' 735 else if (which_kaon.eq.2) then 736 write(iun,) ' **---- producing a SIGMA- ----*' 737 else if (which_kaon.eq.10) then 738 write(iun,) ' **---- WITH BOUND LAMBDA ----*' 739 else if (which_kaon.eq.11) then 740 write(iun,) ' **---- WITH BOUND SIGMA0 ----*' 741 else if (which_kaon.eq.12) then 742 write(iun,) ' **---- WITH BOUND SIGMA- ----*' 743 gaskelld 1.1 else 744 stop 'I don''t have ANY idea what (e,e''K) we''re doing!!!' 745 endif 746 else if (doing_phsp) then 747 write(iun,) ' **--- PHASE SPACE - NO physics, NO radiation (may not work)---*' 748 else 749 stop 'I don''t have ANY idea what we''re doing!!!' 750 endif 751 752 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'Ebeam', Ebeam, 'MeV' 753 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 754 > '(dE/E)beam', dEbeam/Ebeam, '(full wid)' 755 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'x-width', gen.xwid, 'cm' 756 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'y-width', gen.ywid, 'cm' 757 write(iun,'(9x,a12,'' = '',i15,2x,a16)') 758 > 'fr_pattern',targ.fr_pattern, '1=square,2=circ' 759 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'fr1', targ.fr1, 'cm' 760 write(iun,'(9x,a12,'' = '',f15.4,2x,a10)') 'fr2', targ.fr2, 'cm' 761 762 write(iun,) ' ' 763 write(iun,'(9x,18x,2(a15,2x))') '____E arm____','____P arm____' 764 gaskelld 1.1 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'angle', 765 > spec.e.thetadegrad, spec.p.thetadegrad, 'deg' 766 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'momentum', 767 > spec.e.p, spec.p.p, 'MeV/c' 768 769 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'x offset', 770 > spec.e.offset.x, spec.p.offset.x, 'cm' 771 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'y offset', 772 > spec.e.offset.y, spec.p.offset.y, 'cm' 773 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'z offset', 774 > spec.e.offset.z, spec.p.offset.z, 'cm' 775 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'xptar offset', 776 > spec.e.offset.xptar, spec.p.offset.xptar, 'mr' 777 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'yptar offset', 778 > spec.e.offset.yptar, spec.p.offset.yptar, 'mr' 779 780 write(iun,) ' VALUES FOR "CENTRAL" EVENT:' 781 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'delta', 782 > central.e.delta, central.p.delta, '%' 783 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'xptar', 784 > central.e.xptar, central.p.xptar, 'mr' 785 gaskelld 1.1 write(iun,'(9x,a12,'' = '',2(f15.4,2x),2x,a5)') 'yptar', 786 > central.e.yptar, central.p.yptar, 'mr' 787 788 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'Q2',central.Q2/1.d6,'(GeV/c)^2' 789 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'q', central.q, 'MeV/c' 790 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'nu', 791 > central.nu, 'MeV' 792 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'recon Em', central.Em, 'MeV' 793 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'recon Pm', central.Pm, 'MeV/c' 794 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'recon W', central.W, 'MeV/c' 795 write(iun,'(17x,a10,'' = '',f15.4,2x,a9)') 'recon MM', central.MM, 'MeV/c' 796 797 ! Target 798 write(iun,) 'TARGET specs:' 799 9911 format(2x,2(5x,a10,' = ',e12.6,1x,a5)) 800 write(iun,9911) 'A', targ.A, ' ', 'Z', targ.Z, ' ' 801 write(iun,9911) 'mass', targ.mass_amu, 'amu', 'mass', targ.M,'MeV' 802 write(iun,9911) 'Mrec', targ.mrec_amu, 'amu', 'Mrec', targ.Mrec,'MeV' 803 write(iun,9911) 'Mtar_struc',targ.Mtar_struck,'MeV', 804 > 'Mrec_struc',targ.Mrec_struck,'MeV' 805 write(iun,9911) 'rho', targ.rho, 'g/cm3', 'thick', targ.thick,'g/cm2' 806 gaskelld 1.1 write(iun,9911) 'angle', targ.angledegrad, 'deg', 'abundancy', 807 > targ.abundancy, '%' 808 write(iun,9911) 'X0', targ.X0, 'g/cm2', 'X0_cm', targ.X0_cm,'cm' 809 write(iun,9911) 'length',targ.length,'cm','zoffset',targ.zoffset,'cm' 810 write(iun,9911) 'xoffset',targ.xoffset,'cm','yoffset',targ.yoffset,'cm' 811 write(iun,'(t12,3a15)') '__ave__', '__lo__', '__hi__' 812 9912 format(1x,a15,3f15.5,2x,a6) 813 write(iun,9912) 'Coulomb', targ.Coulomb.ave, targ.Coulomb.min, 814 > targ.Coulomb.max, 'MeV' 815 write(iun,9912) 'Eloss_beam', targ.Eloss(1).ave, 816 > targ.Eloss(1).min, targ.Eloss(1).max, 'MeV' 817 write(iun,9912) 'Eloss_e', targ.Eloss(2).ave, targ.Eloss(2).min, 818 > targ.Eloss(2).max, 'MeV' 819 write(iun,9912) 'Eloss_p', targ.Eloss(3).ave, targ.Eloss(3).min, 820 > targ.Eloss(3).max, 'MeV' 821 write(iun,9912) 'teff_beam', targ.teff(1).ave, targ.teff(1).min, 822 > targ.teff(1).max, 'radlen' 823 write(iun,9912) 'teff_e', targ.teff(2).ave, targ.teff(2).min, 824 > targ.teff(2).max, 'radlen' 825 write(iun,9912) 'teff_p', targ.teff(3).ave, targ.teff(3).min, 826 > targ.teff(3).max, 'radlen' 827 gaskelld 1.1 9913 format(1x,a15,t25,f15.5,2x,a6) 828 write(iun,9913) 'musc_nsig_max', targ.musc_nsig_max, ' ' 829 write(iun,9913) 'musc_max_beam', targ.musc_max(1)1000., 'mr' 830 write(iun,9913) 'musc_max_e', targ.musc_max(2)1000., 'mr' 831 write(iun,9913) 'musc_max_p', targ.musc_max(3)1000., 'mr' 832 833 ! Flags 834 write(iun,) 'FLAGS:' 835 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_eep', doing_eep, 836 > 'doing_kaon', doing_kaon, 'doing_pion', doing_pion 837 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_semi', doing_semi, 838 > 'doing_rho', doing_rho, 'doing_hplus', doing_hplus 839 write(iun,'(5x,2(2x,a19,''='',l2))') 'doing_semipi',doing_semipi, 840 > 'doing_semika', doing_semika 841 write(iun,'(5x,2(2x,a19,''='',l2))') 'doing_delta',doing_delta, 842 > 'doing_phsp', doing_phsp 843 write(iun,'(5x,2(2x,a19,''='',l2))') 'which_pion', which_pion, 844 > 'which_kaon', which_kaon 845 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hyd_elast', doing_hyd_elast, 846 > 'doing_deuterium', doing_deuterium, 'doing_heavy', doing_heavy 847 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydpi', doing_hydpi, 848 gaskelld 1.1 > 'doing_deutpi', doing_deutpi, 'doing_hepi', doing_hepi 849 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydkaon', doing_hydkaon, 850 > 'doing_deutkaon', doing_deutkaon, 'doing_hekaon', doing_hekaon 851 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydsemi', doing_hydsemi, 852 > 'doing_deutsemi', doing_deutsemi, 'do_fermi', do_fermi 853 write(iun,'(5x,3(2x,a19,''='',l2))') 'doing_hydrho', doing_hydrho, 854 > 'doing_deutrho', doing_deutrho, 'doing_herho', doing_herho 855 write(iun,'(5x,(2x,a19,''='',l2),2(2x,a19,''='',i2))') 'mc_smear', 856 > mc_smear,'electron_arm',electron_arm,'hadron_arm',hadron_arm 857 write(iun,'(5x,3(2x,a19,''='',l2)))') 'using_Eloss', using_Eloss, 858 > 'using_Coulomb',using_Coulomb,'deForest_flag',deForest_flag 859 write(iun,'(5x,3(2x,a19,''='',l2)))') 'correct_Eloss', correct_Eloss, 860 > 'correct_raster',correct_raster, 'doing_decay', doing_decay 861 write(iun,'(5x,2(2x,a19,''='',l2))') 862 > 'using_E_arm_montecarlo', using_E_arm_montecarlo, 863 > 'using_P_arm_montecarlo', using_P_arm_montecarlo 864 if (electron_arm.eq.5 .or. hadron_arm.eq.5 .or. 865 > electron_arm.eq.6 .or. hadron_arm.eq.6) 866 > write(iun,'(7x,a19,''='',l2)') 'use_first_cer',use_first_cer 867 write(iun,'(7x,a11,''='',f10.3,a4)') 'ctau',ctau,'cm' 868 869 gaskelld 1.1 ! Counters 870 write(iun,) 'COUNTERS:' 871 write(iun,'(12x,''Ngen (request) = '',i10)') ngen 872 write(iun,'(12x,''Ntried = '',i10)') ntried 873 write(iun,'(12x,''Ncontribute = '',i10)') ncontribute 874 write(iun,'(12x,''Nco_no_rad_prot= '',i10)') ncontribute_no_rad_proton 875 write(iun,'(12x,''-> %no_rad_prot= '',f10.3)') 876 > (100.ncontribute_no_rad_proton/max(dble(ncontribute),0.1d0)) 877 write(iun,'(/1x,''INTEGRATED WEIGHTS (number of counts in delta/Em cuts!):'')') 878 write(iun,'(1x,'' MeV: wtcontr= '',e16.8)') wtcontribute/nevent 879 880 ! Radiative corrections 881 write(iun,) 'RADIATIVE CORRECTIONS:' 882 if (.not.using_rad) then 883 write(iun,'(x,a14,''='',l3)') 'using_rad', using_rad 884 else 885 write(iun,'(4(x,a14,''='',l3))') 'use_expon',use_expon, 886 > 'include_hard',include_hard,'calc_spence',calculate_spence 887 write(iun,'(2(x,a14,''='',l3))') 'using_rad', using_rad, 888 > 'use_offshell_rad', use_offshell_rad 889 write(iun,'(4(x,a14,''='',i3))') 'rad_flag',rad_flag, 890 gaskelld 1.1 > 'extrad_flag', extrad_flag, 'one_tail', one_tail, 891 > 'intcor_mode', intcor_mode 892 write(iun,'(x,a14,''='',f11.3)') 'dE_edge_test',dE_edge_test 893 write(iun,'(x,a14,''='',f11.3)') 'Egamma_max', Egamma_tot_max 894 895 9914 format(1x,a18,' = ',4f11.3) 896 write(iun,) 'Central Values:' 897 write(iun,9914) 'hardcorfac', central.rad.hardcorfac 898 write(iun,9914) 'etatzai', central.rad.etatzai 899 write(iun,9914) 'frac(1:3)', central.rad.frac 900 write(iun,9914) 'lambda(1:3)', central.rad.lambda 901 write(iun,9914) 'bt(1:2)', central.rad.bt 902 write(iun,9914) 'c_int(0:3)', central.rad.c_int 903 write(iun,9914) 'c_ext(0:3)', central.rad.c_ext 904 write(iun,9914) 'c(0:3)', central.rad.c 905 write(iun,9914) 'g_int', central.rad.g_int 906 write(iun,9914) 'g_ext', central.rad.g_ext 907 write(iun,9914) 'g(0:3)', central.rad.g 908 endif 909 910 ! Miscellaneous 911 gaskelld 1.1 write(iun,) 'MISCELLANEOUS:' 912 9915 format(12x,a14,' = ',e16.6,1x,a6) 913 write(iun,) 'Note that central.sigcc is for central delta,theta,phi in both spectrometers' 914 write(iun,) ' and may give non-physical kinematics, esp. for Hydrogen' 915 write(iun,) 'Note also that AVE.sigcc is really AVER.weight (the two arenot exactly equal)' 916 write(iun,9915) 'CENTRAL.sigcc', central.sigcc, ' ' 917 write(iun,9915) 'AVERAGE.sigcc', sum_sigcc/nevent, ' ' 918 write(iun,9915) 'charge', EXPER.charge, 'mC' 919 write(iun,9915) 'targetfac', targetfac, ' ' 920 write(iun,9915) 'luminosity', luminosity, 'ub^-1' 921 write(iun,9915) 'luminosity', luminosity(hbarc/100000.)*2, 'GeV^2' 922 write(iun,9915) 'genvol', genvol, ' ' 923 write(iun,9915) 'normfac', normfac, ' ' 924 9916 format(12x,a14,' = ',f6.1,' to ',f6.1,' in ',i4,' bins of ',f7.2,1x,a5) 925 if (doing_heavy) then 926 write(iun,'(12x,''Theory file: '',a)') 927 > theory_file(1:index(theory_file,' ')-1) 928 endif 929 930 ! Resolution summary 931 932 gaskelld 1.1 write(iun,) 'RECON SUMMARY: Ave.Error Resolution' 933 write(iun,99165) 'Electron arm: delta =',10.aveerr.e.delta,10.resol.e.delta,'x10^-3' 934 write(iun,99165) 'xptar =',1000.aveerr.e.xptar,1000.resol.e.xptar,'mr' 935 write(iun,99165) 'yptar =',1000.aveerr.e.yptar,1000.resol.e.yptar,'mr' 936 write(iun,99165) 'ytar =',10.aveerr.e.ytar,10.resol.e.ytar ,'mm' 937 write(iun,99165) 'Hadron arm: delta =',10.aveerr.p.delta,10.resol.p.delta,'x10-3' 938 write(iun,99165) 'xptar =',1000.aveerr.p.xptar,1000.resol.p.xptar,'mr' 939 write(iun,99165) 'yptar =',1000.aveerr.p.yptar,1000.resol.p.yptar,'mr' 940 write(iun,99165) 'ytar =',10.aveerr.p.ytar,10.resol.p.ytar,'mm' 941 99165 format(2x,a22,2f12.5,2x,a) 942 943 944 ! Input spectrometer limits. 945 946 write(iun,) 'Input Spectrometer Limits:' 947 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.delta.min/max', 948 > SPedge.e.delta.min,SPedge.e.delta.max,'%' 949 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.yptar.min/max', 950 > SPedge.e.yptar.min,SPedge.e.yptar.max,'rad' 951 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.e.xptar.min/max', 952 > SPedge.e.xptar.min,SPedge.e.xptar.max,'rad' 953 gaskelld 1.1 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.delta.min/max', 954 > SPedge.p.delta.min,SPedge.p.delta.max,'%' 955 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.yptar.min/max', 956 > SPedge.p.yptar.min,SPedge.p.yptar.max,'rad' 957 write(iun,'(9x,a25,'' = '',2(2x,f15.4),a5)') 'SPedge.p.xptar.min/max', 958 > SPedge.p.xptar.min,SPedge.p.xptar.max,'rad' 959 960 961 ! Edges used in generation and checking, as well as range of events found 962 ! to contribute within those edges 963 964 ! ... on GENERATED qties 965 write(iun,) 'Limiting VERTEX values (vertex.e/p.,Em,Pm,Trec)' 966 write(iun,) ' USED limits are gen.e/p., and VERTEXedge.Em,Pm,Trec' 967 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') '______used______', 968 > '_____found______' 969 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' 970 9917 format(1x,a18,t21,2f12.3,t50,2f10.3,2x,a5) 971 write(iun,9917) 'E arm delta', gen.e.delta.min, gen.e.delta.max, 972 > contrib.gen.e.delta.lo, contrib.gen.e.delta.hi, '%' 973 write(iun,9917) 'E arm yptar', gen.e.yptar.min1000., 974 gaskelld 1.1 > gen.e.yptar.max1000., 975 > contrib.gen.e.yptar.lo1000., contrib.gen.e.yptar.hi1000.,'mr' 976 write(iun,9917) 'E arm xptar', gen.e.xptar.min1000., 977 > gen.e.xptar.max1000., 978 > contrib.gen.e.xptar.lo1000., contrib.gen.e.xptar.hi1000.,'mr' 979 write(iun,9917) 'P arm delta', gen.p.delta.min, gen.p.delta.max, 980 > contrib.gen.p.delta.lo, contrib.gen.p.delta.hi, '%' 981 write(iun,9917) 'P arm yptar', gen.p.yptar.min1000., 982 > gen.p.yptar.max1000., 983 > contrib.gen.p.yptar.lo1000., contrib.gen.p.yptar.hi1000.,'mr' 984 write(iun,9917) 'P arm xptar', gen.p.xptar.min1000., 985 > gen.p.xptar.max1000., 986 > contrib.gen.p.xptar.lo1000., contrib.gen.p.xptar.hi1000.,'mr' 987 write(iun,9917) 'sumEgen', gen.sumEgen.min, gen.sumEgen.max, 988 > contrib.gen.sumEgen.lo, contrib.gen.sumEgen.hi, 'MeV' 989 990 ! ... on VERTEX qties 991 ! write(iun,) 'Limiting VERTEX values:' 992 ! write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') 993 ! > '______used______', '_____found______' 994 ! write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' 995 gaskelld 1.1 write(iun,9917) 'Trec', VERTEXedge.Trec.min, VERTEXedge.Trec.max, 996 > contrib.vertex.Trec.lo, contrib.vertex.Trec.hi, 'MeV' 997 write(iun,9917) 'Em', VERTEXedge.Em.min, VERTEXedge.Em.max, 998 > contrib.vertex.Em.lo, contrib.vertex.Em.hi, 'MeV' 999 write(iun,9917) 'Pm', VERTEXedge.Pm.min, VERTEXedge.Pm.max, 1000 > contrib.vertex.Pm.lo, contrib.vertex.Pm.hi, 'MeV/c' 1001 if ((doing_deuterium .or. doing_pion .or. doing_kaon .or. doing_delta) .and. using_rad) 1002 > write(iun,) ' ** NOTE: sumEgen.min only used in GENERATE_RAD' 1003 1004 ! ... on TRUE qties 1005 write(iun,) 'Limiting ORIGINAL values: orig.e/p.,Em,Pm,Trec (no edge.* limits for Pm,Trec)' 1006 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') 1007 > '______used______', '_____found______' 1008 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' 1009 write(iun,9917) 'E arm E', edge.e.E.min, edge.e.E.max, 1010 > contrib.tru.e.E.lo, contrib.tru.e.E.hi, 'MeV' 1011 write(iun,9917) 'E arm yptar', edge.e.yptar.min1000., 1012 > edge.e.yptar.max1000., 1013 > contrib.tru.e.yptar.lo1000., contrib.tru.e.yptar.hi1000.,'mr' 1014 write(iun,9917) 'E arm xptar', edge.e.xptar.min1000., edge.e.xptar.max1000., 1015 > contrib.tru.e.xptar.lo1000., contrib.tru.e.xptar.hi1000., 'mr' 1016 gaskelld 1.1 write(iun,9917) 'P arm E', edge.p.E.min, edge.p.E.max, 1017 > contrib.tru.p.E.lo, contrib.tru.p.E.hi, 'MeV' 1018 write(iun,9917) 'P arm yptar', edge.p.yptar.min1000., 1019 > edge.p.yptar.max1000., 1020 > contrib.tru.p.yptar.lo1000., contrib.tru.p.yptar.hi1000.,'mr' 1021 write(iun,9917) 'P arm xptar', edge.p.xptar.min1000., edge.p.xptar.max1000., 1022 > contrib.tru.p.xptar.lo1000., contrib.tru.p.xptar.hi1000., 'mr' 1023 write(iun,9917) 'Em', max(-999999.999d0,edge.Em.min), min(999999.999d0,edge.Em.max), 1024 > contrib.tru.Em.lo, contrib.tru.Em.hi, 'MeV' 1025 write(iun,9917) 'Pm', 0., 0., contrib.tru.Pm.lo, 1026 > contrib.tru.Pm.hi, 'MeV' 1027 write(iun,9917) 'Trec', 0., 0., contrib.tru.Trec.lo, 1028 > contrib.tru.Trec.hi, 'MeV' 1029 1030 ! ... on SPECTROMETER qties 1031 write(iun,) 'Limiting SPECTROMETER values' 1032 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') 1033 > '______used______', '_____found______' 1034 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' 1035 write(iun,9917) 'E arm delta', SPedge.e.delta.min, 1036 > SPedge.e.delta.max, contrib.SP.e.delta.lo, 1037 gaskelld 1.1 > contrib.SP.e.delta.hi, '%' 1038 write(iun,9917) 'E arm yptar', SPedge.e.yptar.min1000., 1039 > SPedge.e.yptar.max1000., 1040 > contrib.SP.e.yptar.lo1000., contrib.SP.e.yptar.hi1000., 'mr' 1041 write(iun,9917) 'E arm xptar', SPedge.e.xptar.min1000., 1042 > SPedge.e.xptar.max1000., 1043 > contrib.SP.e.xptar.lo1000., contrib.SP.e.xptar.hi1000., 'mr' 1044 write(iun,9917) 'P arm delta', SPedge.p.delta.min, 1045 > SPedge.p.delta.max, contrib.SP.p.delta.lo, 1046 > contrib.SP.p.delta.hi, '%' 1047 write(iun,9917) 'P arm yptar', SPedge.p.yptar.min1000., 1048 > SPedge.p.yptar.max1000., 1049 > contrib.SP.p.yptar.lo1000., contrib.SP.p.yptar.hi1000., 'mr' 1050 write(iun,9917) 'P arm xptar', SPedge.p.xptar.min1000., 1051 > SPedge.p.xptar.max1000., 1052 > contrib.SP.p.xptar.lo1000., contrib.SP.p.xptar.hi1000., 'mr' 1053 1054 ! ... on RADIATION qties 1055 if (using_rad) then 1056 write(iun,) 'Limiting RADIATION values CONTRIBUTING to the (Em,Pm) distributions:' 1057 write(iun,'(t25,2x,a16,2x,t50,2x,a16,2x)') 1058 gaskelld 1.1 > '______used______', '_____found______' 1059 write(iun,'(t25,2a10,t50,2a10)') 'min','max', 'lo', 'hi' 1060 write(iun,9917) 'Egamma(1)', 0., Egamma1_max, 1061 > contrib.rad.Egamma(1).lo, contrib.rad.Egamma(1).hi, 'MeV' 1062 write(iun,9917) 'Egamma(2)', 0., Egamma2_max, contrib.rad.Egamma(2).lo, 1063 > contrib.rad.Egamma(2).hi, 'MeV' 1064 write(iun,9917) 'Egamma(3)', 0., Egamma3_max, contrib.rad.Egamma(3).lo, 1065 > contrib.rad.Egamma(3).hi, 'MeV' 1066 write(iun,9917) 'Egamma_total', 0., Egamma_tot_max, 1067 > contrib.rad.Egamma_total.lo, contrib.rad.Egamma_total.hi,'MeV' 1068 endif 1069 1070 ! ... on slops 1071 write(iun,) 'ACTUAL and LIMITING SLOP values used/obtained:' 1072 write(iun,'(t25,3a10)') '__used__', '__min__', '__max__' 1073 9918 format(1x,a10,a12,t25,3f10.3,2x,a5) 1074 write(iun,9918) 'slop.MC ', 'E arm delta', slop.MC.e.delta.used, 1075 > slop.MC.e.delta.lo, slop.MC.e.delta.hi, '%' 1076 write(iun,9918) ' ', 'E arm yptar', slop.MC.e.yptar.used1000., 1077 > slop.MC.e.yptar.lo1000., slop.MC.e.yptar.hi1000., 'mr' 1078 write(iun,9918) ' ', 'E arm xptar', slop.MC.e.xptar.used1000., 1079 gaskelld 1.1 > slop.MC.e.xptar.lo1000., slop.MC.e.xptar.hi1000., 'mr' 1080 write(iun,9918) ' ', 'P arm delta', slop.MC.p.delta.used, 1081 > slop.MC.p.delta.lo, slop.MC.p.delta.hi, '%' 1082 write(iun,9918) ' ', 'P arm yptar', slop.MC.p.yptar.used1000., 1083 > slop.MC.p.yptar.lo1000., slop.MC.p.yptar.hi1000., 'mr' 1084 write(iun,9918) ' ', 'P arm xptar', slop.MC.p.xptar.used1000., 1085 > slop.MC.p.xptar.lo1000., slop.MC.p.xptar.hi1000., 'mr' 1086 write(iun,9918) 'slop.total', 'Em', slop.total.Em.used, 1087 > slop.total.Em.lo, slop.total.Em.hi, 'MeV' 1088 write(iun,9918) ' ', 'Pm', 0., slop.total.Pm.lo, 1089 > slop.total.Pm.hi, 'MeV/c' 1090 1091 write(iun,'(/)') 1092 return 1093 end 1094 1095 !----------------------------------------------------------------------- 1096 1097 subroutine calculate_central(central,vertex0) 1098 1099 implicit none 1100 gaskelld 1.1 include 'simulate.inc' 1101 include 'radc.inc' 1102 1103 integer i 1104 real8 m_spec 1105 logical success 1106 record /event_main/ main0 1107 record /event/ vertex0, recon0 1108 record /event_central/ central 1109 1110 ! JRA 2002: Redo this so that central values are chosen as in generate, 1111 ! and then call complete_ev and complete_recon_ev, just like a normal event. 1112 1113 if (debug(2)) write(6,)'calc_cent: entering...' 1114 main0.target.x = 0.0 1115 main0.target.y = 0.0 1116 main0.target.z = targ.zoffset 1117 main0.target.rastery = 0.0 1118 main0.target.Eloss(1) = targ.Eloss(1).ave 1119 main0.target.Coulomb = targ.Coulomb.ave 1120 main0.target.teff(1) = targ.teff(1).ave 1121 gaskelld 1.1 vertex0.Ein = Ebeam_vertex_ave 1122 main0.Ein_shift = 0.0 1123 main0.Ee_shift = 0.0 1124 main0.gen_weight = 1.0 1125 1126 ! Set all of these to central values, but then complete_ev will force 1127 ! the variables that are not normally generated to their appropriate values. 1128 vertex0.e.delta = 0.0 1129 vertex0.e.yptar = 0.0 1130 vertex0.e.xptar = 0.0 1131 vertex0.e.theta = spec.e.theta 1132 vertex0.e.phi = spec.e.phi 1133 vertex0.e.P = spec.e.P(1.+vertex0.e.delta/100.) 1134 vertex0.e.E = vertex0.e.P 1135 1136 vertex0.p.delta = 0.0 1137 vertex0.p.yptar = 0.0 1138 vertex0.p.xptar = 0.0 1139 vertex0.p.theta = spec.p.theta 1140 vertex0.p.phi = spec.p.phi 1141 vertex0.p.P = spec.p.P(1.+vertex0.p.delta/100.) 1142 gaskelld 1.1 vertex0.p.E = sqrt(vertex0.p.P2 + Mh2) 1143 1144 pfer = 0.0 1145 pferx = 0.0 1146 pfery = 0.0 1147 pferz = 0.0 1148 vertex0.Em = targ.Mtar_struck + targ.Mrec - targ.M 1149 m_spec = targ.M - targ.Mtar_struck + vertex0.Em !=targ.Mrec 1150 efer = targ.M - sqrt(m_spec2+pfer2) !=M-Mrec 1151 1152 ! Old version - not appropriate for all event types. 1153 ! Pop in generation values for central event 1154 ! 1155 ! if (debug(2)) write(6,)'calc_cent: entering...' 1156 ! main0.target.x = 0.0 1157 ! main0.target.y = 0.0 1158 ! main0.target.z = 0.0 1159 ! vertex0.Ein = Ebeam_vertex_ave 1160 ! main0.target.Coulomb = targ.Coulomb.ave 1161 ! main0.target.Eloss(1) = targ.Eloss(1).ave 1162 ! main0.target.teff(1) = targ.teff(1).ave 1163 gaskelld 1.1 ! vertex0.e.delta = 0.0 1164 ! vertex0.e.yptar = 0.0 1165 ! vertex0.e.xptar = 0.0 1166 ! vertex0.e.theta = spec.e.theta 1167 ! vertex0.e.phi = spec.e.phi 1168 ! vertex0.e.P = spec.e.P(1.+vertex0.e.delta/100.) 1169 ! vertex0.e.E = vertex0.e.P 1170 ! vertex0.p.delta = 0.0 1171 ! vertex0.p.yptar = 0.0 1172 ! vertex0.p.xptar = 0.0 1173 ! vertex0.p.theta = spec.p.theta 1174 ! vertex0.p.phi = spec.p.phi 1175 ! vertex0.p.P = spec.p.P(1.+vertex0.p.delta/100.) 1176 ! vertex0.p.E = sqrt(vertex0.p.P*2 + Mh2) 1177 1178 ! Complete_recon_ev for vertex0. Note: complete_recon_ev doesn't 1179 ! call radc_init_ev or calculate teff(2,3). 1180 1181 ! JRA: Do we want complete_ev or complete_recon_ev? Do we want to calculate 1182 ! and/or dump other central values (for pion/kaon production, for example). 1183 1184 gaskelld 1.1 c call complete_ev(main0,vertex0,success) 1185 c if (debug(2)) write(6,)'calc_cent: done with complete_ev' 1186 c if (.not.success) stop 'COMPLETE_EV failed trying to complete a CENTRAL event!' 1187 1188 call complete_recon_ev(vertex0,success) 1189 if (debug(2)) write(6,)'calc_cent: done with complete_recon_ev' 1190 if (.not.success) stop 'COMPLETE_EV failed trying to complete a CENTRAL event!' 1191 central.Q2 = vertex0.Q2 1192 central.q = vertex0.q 1193 central.nu = vertex0.nu 1194 central.Em = vertex0.Em 1195 central.Pm = vertex0.Pm 1196 central.W = vertex0.W 1197 central.e.xptar = vertex0.e.xptar 1198 central.e.yptar = vertex0.e.yptar 1199 central.e.delta = vertex0.e.delta 1200 central.p.xptar = vertex0.p.xptar 1201 central.p.yptar = vertex0.p.yptar 1202 central.p.delta = vertex0.p.delta 1203 1204 if (central.Em2-central.Pm2 .lt. 0) then 1205 gaskelld 1.1 central.MM = -sqrt(abs(central.Em2-central.Pm2)) 1206 else 1207 central.MM = sqrt(central.Em2-central.Pm2) 1208 endif 1209 1210 main0.target.teff(2) = targ.teff(2).ave 1211 main0.target.teff(3) = targ.teff(3).ave 1212 if (debug(2)) write(6,)'calc_cent: calling radc_init_ev' 1213 if (debug(4)) write(6,)'calc_cent: Ein =',vertex0.Ein 1214 call radc_init_ev(main0,vertex0) 1215 if (debug(2)) write(6,)'calc_cent: done with radc_init_ev' 1216 if (using_rad) then 1217 central.rad.hardcorfac = hardcorfac 1218 central.rad.etatzai= etatzai 1219 central.rad.g_int = g_int 1220 central.rad.g_ext = g_ext 1221 do i = 0, 3 1222 if (i.gt.0) central.rad.frac(i) = frac(i) 1223 if (i.gt.0) central.rad.lambda(i) = lambda(i) 1224 if (i.gt.0.and.i.lt.3) central.rad.bt(i) = bt(i) 1225 central.rad.c_int(i) = c_int(i) 1226 gaskelld 1.1 central.rad.c_ext(i) = c_ext(i) 1227 central.rad.c(i) = c(i) 1228 central.rad.g(i) = g(i) 1229 enddo 1230 endif 1231 if (debug(4)) write(6,*)'calc_cent: at 1' 1232 1233 do i = 1, neventfields 1234 recon0.all(i) = vertex0.all(i) 1235 enddo 1236 call complete_main(.true.,main0,vertex0,vertex0,recon0,success) 1237 central.sigcc = main0.sigcc 1238
1239 gaskelld 1.3 write(6,) 'central event' 1240 write(6,) 'Pt',sqrt(vertex0.pt2)/1.d3 1241 write(6,) 'z', vertex0.zhad 1242 write(6,) 'lab cross section (nb/Gev2/sr2)',central.sigcc1000.01000.0*1000.0 1243
1244 gaskelld 1.1 if (debug(2)) write(6,)'calc_cent: ending...' 1245 return 1246 end 1247 1248 !------------------------------------------------------------------------- 1249 1250 subroutine montecarlo(orig,main,recon,success) 1251 1252 implicit none 1253 include 'simulate.inc' 1254 1255 Track-coordinate and spectrometer common blocks 1256 1257 real8 x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm,delta_E_arm 1258 real8 x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm,delta_P_arm 1259 real8 xfp, yfp, dxfp, dyfp 1260 real8 eloss_E_arm, eloss_P_arm, r, beta, dangles(2), dang_in(2) 1261 logical success 1262 logical ok_E_arm, ok_P_arm 1263 record /event/ orig, recon 1264 record /event_main/ main 1265 gaskelld 1.1 1266 real8 beta_electron 1267 parameter (beta_electron = 1.) 1268 real8 tmpfact 1269 real*8 fry !fast raster y position.
1270 gaskelld 1.4 real*8 frx !fast raster x position - left as looking downstream
1271 gaskelld 1.1 real8 m2 !mass for call to mc_hms(sos). Changes if decay 1272 real8 pathlen 1273
1274 gaskelld 1.3 real*8 dx_tmp,dy_tmp 1275
1276 gaskelld 1.4 real*8 ctheta,stheta,phad,pelec 1277
1278 gaskelld 1.1 real*8 zero 1279 parameter (zero=0.0d0) !double precision zero for subroutine calls 1280 1281 ! Prepare the event for the Monte Carlo's and/or spectrometer cuts 1282 1283 success = .false. 1284 ntup.resfac = 0.0 !resfac (see simulate.inc) 1285 if (correct_raster) then 1286 fry = -main.target.rastery
1287 gaskelld 1.4 frx = -main.target.rasterx !This should make frx positive for beam left
1288 gaskelld 1.1 else 1289 fry = 0.0
1290 gaskelld 1.4 frx = 0.0
1291 gaskelld 1.1 endif 1292 1293 !BEAM MULTIPLE SCATTERING: NOT YET IMPLEMENTED, JUST GETTING IT READY! 1294 1295 ! ... multiple scattering of the beam. Generate angles for beam deflection, 1296 ! ... and then apply those angles to the electron and proton. 1297 ! ... The yptar offset goes directly to yptar of both arms (small angle approx). 1298 ! ... xptar is multiplied by cos(theta) to get the xptar offsets. 1299 1300 if (mc_smear) then 1301 call target_musc(orig.Ein,beta_electron,main.target.teff(1),dang_in) 1302 else 1303 dang_in(1)=0.0 !yp offset, goes directly to yp of both arms 1304 dang_in(2)=0.0 !xp offset, mult. by cos_th for xp of both arms 1305 endif 1306 1307 if (debug(3)) write(6,) 'mc: using p,e arm mc =', 1308 > using_P_arm_montecarlo,using_E_arm_montecarlo 1309 1310 !____ P arm ________________ 1311 1312 gaskelld 1.1 ! Go from TRUE to SPECTROMETER quantities by computing target distortions 1313 ! ... ionization loss correction (if requested) 1314 1315 if (using_Eloss) then 1316 if (debug(3)) write(6,)'mc: p arm stuff0 =', 1317 > orig.p.E,main.target.Eloss(3),spec.p.P 1318 main.SP.p.delta = (sqrt(abs((orig.p.E-main.target.Eloss(3))*2 1319 > -Mh2))-spec.p.P) / spec.p.P100. 1320 else 1321 if (debug(3)) write(6,)'mc: p arm stuff1 =',orig.p.delta 1322 main.SP.p.delta = orig.p.delta 1323 endif 1324 1325 ! ... multiple scattering 1326 1327 if (mc_smear) then 1328 beta = orig.p.p/orig.p.E 1329 call target_musc(orig.p.p, beta, main.target.teff(3), dangles) 1330 else 1331 dangles(1)=0.0 1332 dangles(2)=0.0 1333 gaskelld 1.1 endif 1334 main.SP.p.yptar = orig.p.yptar + dangles(1) + dang_in(1) 1335 main.SP.p.xptar = orig.p.xptar + dangles(2) + dang_in(2)spec.p.cos_th 1336 1337 ! CASE 1: Using the spectrometer Monte Carlo 1338 1339 if (using_P_arm_montecarlo) then 1340 1341 ! ... change to P arm spectrometer coordinates (TRANSPORT system), 1342 1343 if (abs(cos(spec.p.phi)).gt.0.0001) then !phi not at +/- pi/2 1344 write(6,) 'y_P_arm, z_P_arm will be incorrect if spec.p.phi <> pi/2 or 3pi/2' 1345 write(6,) 'spec.p.phi=',spec.p.phi,'=',spec.p.phi180/pi,'degrees' 1346 endif 1347 delta_P_arm = main.SP.p.delta 1348 x_P_arm = -main.target.y 1349 y_P_arm = -main.target.xspec.p.cos_th - main.target.zspec.p.sin_thsin(spec.p.phi) 1350 z_P_arm = main.target.zspec.p.cos_th + main.target.xspec.p.sin_thsin(spec.p.phi) 1351 1352 ! ... Apply spectrometer offset (using spectrometer coordinate system). 1353 x_P_arm = x_P_arm - spec.p.offset.x 1354 gaskelld 1.1 y_P_arm = y_P_arm - spec.p.offset.y 1355 z_P_arm = z_P_arm - spec.p.offset.z 1356 dx_P_arm = main.SP.p.xptar - spec.p.offset.xptar 1357 dy_P_arm = main.SP.p.yptar - spec.p.offset.yptar 1358
1359 gaskelld 1.4 1360 c write(,) 'sign_hms_part =' ,sign_hms_part 1361 if (using_tgt_field) then ! do target field tracking - GAW 1362 if (debug(6)) then 1363 write(,) '------------------------------' 1364 write(,'("frx,fry,x_E_arm = ",3f12.5)') frx,fry,x_P_arm 1365 endif 1366 call track_from_tgt(x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm, 1367 > sign_hadronspec.p.P*(1+delta_P_arm/100.),Mh,1,ok_P_arm) 1368 endif 1369 ! GAW - end 99/11/3 1370
1371 gaskelld 1.1 C DJG need to decay the rho here before we begin transporting through the 1372 C DJG spectrometer 1373 c m2 = Mh2 1374 c if(doing_rho) then 1375 c call rho_decay(dx_P_arm,dy_P_arm,delta_P_arm,spec.p.P,m2, 1376 c 1 main.epsilon,orig.Q2) 1377 c endif 1378 C DJG moved this to the last part of generate!!! 1379 1380 ! ........ drift this position back to z=0, the plane through the target center 1381 1382 x_P_arm = x_P_arm - z_P_armdx_P_arm 1383 y_P_arm = y_P_arm - z_P_armdy_P_arm 1384 z_P_arm = 0.0 1385 1386 main.SP.p.z=y_P_arm 1387 1388 ! ... Monte Carlo through the P arm! if we succeed, continue ... 1389 ! ... Here's what's passed: 1390 ! spectrometer central momentum 1391 ! spectrometer central theta angle 1392 gaskelld 1.1 ! momentum delta 1393 ! x, y, z, theta, and phi in TRANSPORT coordinates 1394 ! x, y, theta, and phi at the focal plane 1395 ! radiation length of target (NOT USED!) 1396 ! particle mass (modified if the hadron decays) 1397 ! multiple scattering flag 1398 ! wcs smearing flag 1399 ! decay flag 1400 ! resmult=resolution factor for DCs (see simulate.inc) 1401 ! vertical position at target (for reconstruction w/raster MEs). 1402 ! ok flag 1403
1404 gaskelld 1.2 m2 = Mh2
1405 gaskelld 1.1 pathlen = 0.0 1406 if (hadron_arm.eq.1) then 1407 call mc_hms(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, 1408 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, 1409 > m2, mc_smear, mc_smear, doing_decay,
1410 gaskelld 1.5 > ntup.resfac, fry, ok_P_arm, pathlen)
1411 gaskelld 1.1 else if (hadron_arm.eq.2) then 1412 call mc_sos(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, 1413 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, 1414 > m2, mc_smear, mc_smear, doing_decay, 1415 > ntup.resfac, fry, ok_P_arm, pathlen) 1416 else if (hadron_arm.eq.3) then 1417 call mc_hrsr(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, 1418 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, 1419 > m2, mc_smear, mc_smear, doing_decay, 1420 > ntup.resfac, fry, ok_P_arm, pathlen) 1421 else if (hadron_arm.eq.4) then 1422 call mc_hrsl(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, 1423 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, 1424 > m2, mc_smear, mc_smear, doing_decay, 1425 > ntup.resfac, fry, ok_P_arm, pathlen) 1426 else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then 1427 call mc_shms(spec.p.P, spec.p.theta, delta_P_arm, x_P_arm, 1428 > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, 1429 > m2, mc_smear, mc_smear, doing_decay, 1430 > ntup.resfac, fry, ok_P_arm, pathlen, hadron_arm, use_first_cer) 1431 endif 1432 gaskelld 1.1
1433 gaskelld 1.4 1434 C DJG Do polarized target field stuff if needed 1435 C DJG Note that the call to track_to_tgt is with -fry: for some reason that routine then 1436 C DJG sets xx=-fry, and calls mc_hms_recon with xx, so it all works out. Whatever. 1437 C DJG To summarize: fry points "down", frx points beam left (facing downstream) 1438 C DJG For spectrometers to the right of the beamline, need to pass ctheta,stheta to track_to_tgt 1439 C DJG For spectrometers to the left of the beamline, need to pass ctheta,-stheta to track_to_tgt 1440 1441 if (using_tgt_field) then 1442 phad = spec.p.P(1.+delta_P_arm/100.0) 1443 ctheta = cos(spec.p.theta) 1444 stheta = sin(spec.p.theta) 1445 if((hadron_arm.eq.1).or.(hadron_arm.eq.3)) then !spectrometers of the right (HMS, HRSR) 1446 call track_to_tgt(delta_P_arm,y_P_arm,dx_P_arm,dy_P_arm,frx,-fry, 1447 1 sign_hadronphad,sqrt(m2),ctheta,stheta,1,ok_P_arm,hadron_arm) 1448 else if((hadron_arm.eq.2).or.(hadron_arm.eq.4).or.(hadron_arm.eq.5)) then !left spects (SOS,HRSL,SHMS) 1449 call track_to_tgt(delta_P_arm,y_P_arm,dx_P_arm,dy_P_arm,frx,-fry, 1450 1 sign_hadronphad,sqrt(m2),ctheta,-stheta,1,ok_P_arm,hadron_arm) 1451 else 1452 write(6,) 'Target field reconstruction not set up for your spectrometer' 1453 stop 1454 gaskelld 1.4 endif 1455 endif 1456 1457
1458 gaskelld 1.1 if (.not.ok_P_arm) then 1459 if (debug(3)) write(6,)'mc: ok_P_arm =',ok_P_arm 1460 return 1461 endif 1462 1463 main.RECON.p.delta = delta_P_arm 1464 main.RECON.p.yptar = dy_P_arm 1465 main.RECON.p.xptar = dx_P_arm 1466 main.RECON.p.z = y_P_arm 1467 main.FP.p.x = xfp 1468 main.FP.p.dx = dxfp 1469 main.FP.p.y = yfp 1470 main.FP.p.dy = dyfp 1471 main.FP.p.path = pathlen 1472 1473 ! CASE 2: Not using the detailed Monte Carlo, just copy the IN event to the 1474 ! OUT record 1475 1476 else 1477 main.RECON.p.delta = main.SP.p.delta 1478 main.RECON.p.yptar = main.SP.p.yptar 1479 gaskelld 1.1 main.RECON.p.xptar = main.SP.p.xptar 1480 endif 1481 1482 ! Fill recon quantities. 1483 1484 recon.p.delta = main.RECON.p.delta 1485 recon.p.yptar = main.RECON.p.yptar 1486 recon.p.xptar = main.RECON.p.xptar 1487 recon.p.z = main.RECON.p.z 1488 recon.p.P = spec.p.P(1.+recon.p.delta/100.) 1489 recon.p.E = sqrt(recon.p.P**2 + Mh2)
1490 gaskelld 1.3 dx_tmp = recon.p.xptar + spec.p.offset.xptar 1491 dy_tmp = recon.p.yptar + spec.p.offset.yptar 1492 1493 call physics_angles(spec.p.theta,spec.p.phi,dx_tmp, 1494 > dy_tmp,recon.p.theta,recon.p.phi)
1495 gaskelld 1.1 1496 ! ... correct for energy loss - use most probable (last flag = 4) 1497 1498 if (correct_Eloss) then 1499 call trip_thru_target (3, zero, recon.p.E, 1500 > recon.p.theta, eloss_P_arm, r,Mh,4) 1501 recon.p.E = recon.p.E + eloss_P_arm 1502 recon.p.E = max(recon.p.E,sqrt(Mh2+0.000001)) !can get P~0 when calculating hadron momentum-->P<0 after eloss 1503 recon.p.P = sqrt(recon.p.E**2-Mh2)
1504 gaskelld 1.3 C DJG Should not correct delta for Eloss - delta is a SPECTROMETER variable 1505 C recon.p.delta = (recon.p.P-spec.p.P)/spec.p.P*100.
1506 gaskelld 1.1 endif 1507 1508 !___ E arm ______________ 1509 1510 ! Go from TRUE to SPECTROMETER quantities by computing target distortions 1511 1512 ! ... ionization loss correction (if requested) and Coulomb deceleration 1513 1514 if (debug(3)) write(6,)'mc: e arm stuff0 =', 1515 > orig.e.E,main.target.Eloss(2),spec.e.P 1516 main.SP.e.delta=100. (orig.e.E - main.target.Eloss(2) 1517 > - main.target.Coulomb - spec.e.P) / spec.e.P 1518 1519 ! ... multiple scattering 1520 1521 if (mc_smear) then 1522 call target_musc(orig.e.p, beta_electron, main.target.teff(2), dangles) 1523 else 1524 dangles(1)=0.0 1525 dangles(2)=0.0 1526 endif 1527 gaskelld 1.1 1528 main.SP.e.yptar = orig.e.yptar + dangles(1) + dang_in(1) 1529 main.SP.e.xptar = orig.e.xptar + dangles(2) + dang_in(2)spec.p.cos_th 1530 1531 ! CASE 1: Using the spectrometer Monte Carlo 1532 1533 if (using_E_arm_montecarlo) then 1534 1535 ! ... change to E arm spectrometer coordinates (TRANSPORT system), 1536 1537 if (abs(cos(spec.e.phi)).gt.0.0001) then !phi not at +/- pi/2 1538 write(6,) 'y_E_arm, z_E_arm will be incorrect if spec.e.phi <> pi/2 or 3pi/2' 1539 write(6,) 'spec.e.phi=',spec.e.phi,'=',spec.e.phi180/pi,'degrees' 1540 endif 1541 delta_E_arm = main.SP.e.delta 1542 x_E_arm = -main.target.y 1543 y_E_arm = -main.target.xspec.e.cos_th - main.target.zspec.e.sin_thsin(spec.e.phi) 1544 z_E_arm = main.target.zspec.e.cos_th + main.target.xspec.e.sin_th*sin(spec.e.phi) 1545 1546 ! ... Apply spectrometer offset (using spectrometer coordinate system). 1547 x_E_arm = x_E_arm - spec.e.offset.x 1548 gaskelld 1.1 y_E_arm = y_E_arm - spec.e.offset.y 1549 z_E_arm = z_E_arm - spec.e.offset.z 1550 dx_E_arm = main.SP.e.xptar - spec.e.offset.xptar 1551 dy_E_arm = main.SP.e.yptar - spec.e.offset.yptar 1552
1553 gaskelld 1.4 ! GAW -project to z=0 to compare with reconstructed target positions 1554 if (using_tgt_field) then ! do target field tracking - GAW 1555 1556 if (debug(6)) then 1557 write(,) '------------------------------' 1558 write(,'("frx,fry,x_E_arm = ",3f12.5)') frx,fry,x_E_arm 1559 endif 1560 call track_from_tgt(x_E_arm,y_E_arm,z_E_arm,dx_E_arm,dy_E_arm, 1561 > -spec.e.P(1+delta_E_arm/100.),Me,-1,ok_E_arm) 1562 endif 1563
1564 gaskelld 1.1 ! ........ drift this position back to z=0, the plane through the target center 1565 1566 x_E_arm = x_E_arm - z_E_armdx_E_arm 1567 y_E_arm = y_E_arm - z_E_armdy_E_arm 1568 z_E_arm = 0.0 1569 1570 main.SP.e.z=y_E_arm 1571 1572 ! ... Monte Carlo through the E arm! if we succeed, continue ... 1573 ! ... Here's what's passed: 1574 ! spectrometer central momentum 1575 ! spectrometer central theta angle 1576 ! momentum delta 1577 ! x, y, z, theta, and phi in TRANSPORT coordinates 1578 ! x, y, theta, and phi at the focal plane 1579 ! radiation length of target (NOT USED!) 1580 ! particle mass 1581 ! multiple scattering flag 1582 ! wcs smearing flag 1583 ! resmult=resolution factor for DCs (see simulate.inc) 1584 ! decay flag (hardwired to .false. for electron arm). 1585 gaskelld 1.1 ! ok flag 1586 1587 m2 = me2 1588 pathlen = 0.0 1589 if (electron_arm.eq.1) then 1590 call mc_hms(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, 1591 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, 1592 > me2, mc_smear, mc_smear, .false., 1593 > tmpfact, fry, ok_E_arm, pathlen) 1594 else if (electron_arm.eq.2) then 1595 call mc_sos(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, 1596 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, 1597 > me2, mc_smear, mc_smear, .false., 1598 > tmpfact, fry, ok_E_arm, pathlen) 1599 else if (electron_arm.eq.3) then 1600 call mc_hrsr(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, 1601 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, 1602 > me2, mc_smear, mc_smear, .false., 1603 > tmpfact, fry, ok_E_arm, pathlen) 1604 else if (electron_arm.eq.4) then 1605 call mc_hrsl(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, 1606 gaskelld 1.1 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, 1607 > me2, mc_smear, mc_smear, .false., 1608 > tmpfact, fry, ok_E_arm, pathlen) 1609 else if (electron_arm.eq.5 .or. electron_arm.eq.6) then 1610 call mc_shms(spec.e.P, spec.e.theta, delta_E_arm, x_E_arm, 1611 > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, 1612 > me2, mc_smear, mc_smear, .false., 1613 > tmpfact, fry, ok_E_arm, pathlen, electron_arm, use_first_cer) 1614 endif 1615 ntup.resfac=ntup.resfac+tmpfact
1616 gaskelld 1.4 1617 1618 C DJG Do polarized target field stuff if needed 1619 C DJG Note that the call to track_to_tgt is with -fry: for some reason that routine then 1620 C DJG sets xx=-fry, and calls mc_hms_recon with xx, so it all works out. Whatever. 1621 C DJG To summarize: fry points "down", frx points beam left (facing downstream) 1622 C DJG For spectrometers to the right of the beamline, need to pass ctheta,stheta to track_to_tgt 1623 C DJG For spectrometers to the left of the beamline, need to pass ctheta,-stheta to track_to_tgt 1624 1625 if (using_tgt_field) then 1626 pelec = spec.e.P(1.+delta_E_arm/100.0) 1627 ctheta = cos(spec.e.theta) 1628 stheta = sin(spec.e.theta) 1629 if((electron_arm.eq.1).or.(electron_arm.eq.3)) then !spectrometers on the right (HMS, HRSR) 1630 call track_to_tgt(delta_E_arm,y_E_arm,dx_E_arm,dy_E_arm,frx,-fry, 1631 1 -1.0pelec,sqrt(me2),ctheta,stheta,-1,ok_E_arm,electron_arm) 1632 else if((electron_arm.eq.2).or.(electron_arm.eq.4).or.(electron_arm.eq.5)) then !left spects(SOS,HRSL,SHMS) 1633 call track_to_tgt(delta_E_arm,y_E_arm,dx_E_arm,dy_E_arm,frx,-fry, 1634 1 -1.0pelec,sqrt(me2),ctheta,-stheta,-1,ok_E_arm,electron_arm) 1635 else 1636 write(6,) 'Target field reconstruction not set up for your spectrometer' 1637 gaskelld 1.4 stop 1638 endif 1639 endif
1640 gaskelld 1.1 1641 if (.not.ok_E_arm) then 1642 if (debug(3)) write(6,)'mc: ok_E_arm =',ok_E_arm 1643 return 1644 endif 1645 1646 main.RECON.e.delta = delta_E_arm 1647 main.RECON.e.yptar = dy_E_arm 1648 main.RECON.e.xptar = dx_E_arm 1649 main.RECON.e.z = y_E_arm 1650 main.FP.e.x = xfp 1651 main.FP.e.dx = dxfp 1652 main.FP.e.y = yfp 1653 main.FP.e.dy = dyfp 1654 main.FP.e.path = pathlen 1655 1656 ! CASE 2: Not using the detailed Monte Carlo, just copy the IN event to the 1657 ! OUT record 1658 1659 else 1660 main.RECON.e.delta = main.SP.e.delta 1661 gaskelld 1.1 main.RECON.e.yptar = main.SP.e.yptar 1662 main.RECON.e.xptar = main.SP.e.xptar 1663 endif 1664 1665 ! Fill recon quantities. 1666 1667 recon.e.delta = main.RECON.e.delta 1668 recon.e.yptar = main.RECON.e.yptar 1669 recon.e.xptar = main.RECON.e.xptar 1670 recon.e.z = main.RECON.e.z 1671 recon.e.P = spec.e.P(1.+recon.e.delta/100.) 1672 recon.e.E = recon.e.P
1673 gaskelld 1.3 1674 dx_tmp = recon.e.xptar + spec.e.offset.xptar 1675 dy_tmp = recon.e.yptar + spec.e.offset.yptar 1676 1677 call physics_angles(spec.e.theta,spec.e.phi,dx_tmp, 1678 > dy_tmp,recon.e.theta,recon.e.phi)
1679 gaskelld 1.1 1680 1681 ! ... correct for energy loss and correct for Coulomb deceleration 1682 1683 if (correct_Eloss) then 1684 call trip_thru_target (2, zero, recon.e.E, recon.e.theta, 1685 > eloss_E_arm, r, Me, 4) 1686 recon.e.E = recon.e.E + eloss_E_arm 1687 endif 1688 recon.e.E = recon.e.E + targ.Coulomb.ave 1689 recon.e.P = recon.e.E
1690 gaskelld 1.3 C DJG Should not correct delta for energy loss - delta is a SPECTROMETER 1691 C variable!! 1692 c recon.e.delta = (recon.e.P-spec.e.P)/spec.e.P*100.
1693 gaskelld 1.1 1694 ! Made it! 1695 success = .true. 1696 1697 return 1698 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood