1 gaskelld 1.1 subroutine limits_update(main,vertex,orig,recon,doing_deuterium,
2 > doing_pion,doing_kaon,doing_delta,doing_rho,contrib,slop)
3
4 include 'structures.inc'
5 include 'radc.inc'
6 record /event_main/ main
7 record /event/ vertex, orig, recon
8 record /contrib/ contrib
9 record /slop/ slop
10 integer i
11 logical doing_deuterium, doing_pion, doing_kaon, doing_delta, doing_rho
12
13 ! Update the "contribution limits" records
14
15 ! ... GENERATION values
16 call update_range(vertex.e.delta, contrib.gen.e.delta)
17 call update_range(vertex.e.yptar, contrib.gen.e.yptar)
18 call update_range(vertex.e.xptar, contrib.gen.e.xptar)
19 call update_range(vertex.p.delta, contrib.gen.p.delta)
20 call update_range(vertex.p.yptar, contrib.gen.p.yptar)
21 call update_range(vertex.p.xptar, contrib.gen.p.xptar)
22 gaskelld 1.1 call update_range(main.Trec, contrib.gen.Trec)
23
24 ! ........ another tricky shift
25 if (doing_deuterium .or. doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho) then
26 call update_range(vertex.e.E-main.Ein_shift,contrib.gen.sumEgen)
27 else
28 call update_range(vertex.e.E+vertex.p.E-main.Ein_shift,contrib.gen.sumEgen)
29 endif
30
31 ! ... TRUE values
32 ! ........ tricky shift here, remember this'll get compared with edge.e.E,
33 ! ........ compensate for main.target.Coulomb to
34 ! ........ copy the shift made to the edge in generate
35 call update_range(orig.e.E-main.Ee_shift, contrib.tru.e.E)
36 call update_range(orig.e.xptar, contrib.tru.e.xptar)
37 call update_range(orig.e.yptar, contrib.tru.e.yptar)
38 call update_range(orig.e.xptar, contrib.tru.e.xptar)
39 call update_range(orig.p.E, contrib.tru.p.E)
40 call update_range(orig.p.yptar, contrib.tru.p.yptar)
41 call update_range(orig.p.xptar, contrib.tru.p.xptar)
42
43 gaskelld 1.1 ! ........ another tricky shift
44 call update_range(orig.Em-main.Ein_shift+main.Ee_shift,contrib.tru.Em)
45 call update_range(orig.Pm, contrib.tru.Pm)
46 call update_range(orig.Trec, contrib.tru.Trec)
47
48 ! ... SPECTROMETER values
49 call update_range(main.SP.e.delta, contrib.SP.e.delta)
50 call update_range(main.SP.e.yptar, contrib.SP.e.yptar)
51 call update_range(main.SP.e.xptar, contrib.SP.e.xptar)
52 call update_range(main.SP.p.delta, contrib.SP.p.delta)
53 call update_range(main.SP.p.yptar, contrib.SP.p.yptar)
54 call update_range(main.SP.p.xptar, contrib.SP.p.xptar)
55
56 ! ... VERTEX values
57 call update_range(vertex.Trec, contrib.vertex.Trec)
58 call update_range(vertex.Em, contrib.vertex.Em)
59 call update_range(vertex.Pm, contrib.vertex.Pm)
60
61 ! ... RADIATION stuff
62 ! ??? should be looking at Egamma2+3 cause we do use limits on that, indirectly
63
64 gaskelld 1.1 do i = 1, 3
65 call update_range(Egamma_used(i), contrib.rad.Egamma(i))
66 enddo
67 call update_range(Egamma_used(1)+Egamma_used(2)+Egamma_used(3),
68 > contrib.rad.Egamma_total)
69
70 ! Update the "slop limits" records
71 ! ... MC slops
72 call update_range(main.RECON.e.delta-main.SP.e.delta,slop.MC.e.delta)
73 call update_range(main.RECON.e.yptar-main.SP.e.yptar,slop.MC.e.yptar)
74 call update_range(main.RECON.e.xptar-main.SP.e.xptar,slop.MC.e.xptar)
75 call update_range(main.RECON.p.delta-main.SP.p.delta,slop.MC.p.delta)
76 call update_range(main.RECON.p.yptar-main.SP.p.yptar,slop.MC.p.yptar)
77 call update_range(main.RECON.p.xptar-main.SP.p.xptar,slop.MC.p.xptar)
78
79 ! ... total slops
80 ! ........ that tricky shift again, slops accounted for by the shift not
81 ! ........ included in slop.total.Em.
82 call update_range(recon.Em-(orig.Em-main.Ein_shift+main.Ee_shift),
83 > slop.total.Em)
84 call update_range(abs(recon.Pm)-abs(orig.Pm), slop.total.Pm)
85 gaskelld 1.1
86 return
87 end
88
89 !-------------------------------------------------------------------
90
91 subroutine update_range(val,range)
92
93 include 'structures.inc'
94 record /range/ range
95 real*8 val
96
97 range.lo = min(range.lo, val)
98 range.hi = max(range.hi, val)
99
100 return
101 end
102
103 !----------------------------------------------------------------------
104 ! THE routine to GENERATE the (max of 7) random quantities we need to get
105 ! a full description of our event.
106 gaskelld 1.1
107 subroutine generate(main,vertex,orig,success)
108
109 implicit none
110 include 'simulate.inc'
111
112 integer i, ii
113 real*8 Emin, Emax
114 real*8 ranprob,ranth1,ranth,ranph,t3,t4,t5,t6
115 real*8 pferlo,pferhi
116 real*8 m_spec !spectator (A-1) mass based on missing energy
117 real*8 gauss1
118 logical success
119 real*8 grnd !random # generator.
120 record /event_main/ main
121 record /event/ vertex, orig
122
123 real*8 nsig_max
124 parameter(nsig_max=3.0d0) !max #/sigma for gaussian ran #s.
125
126
127 gaskelld 1.1 ! Randomize the position of the interaction inside the available region.
128 ! gen.xwid and gen.ywid are the intrinsic beam widths (one sigma value).
129 ! Use a gaussian beam distribution with +/- 3.0 sigma (add raster afterwards).
130 if (debug(2)) write(6,*)'gen: entering'
131 main.target.x = gauss1(nsig_max)*gen.xwid+targ.xoffset
132 main.target.y = gauss1(nsig_max)*gen.ywid+targ.yoffset
133
134 ! fr_pattern=1 - rectangle - targ.fr1/fr2 are the x/y raster half-widths.
135 ! fr_pattern=2 - circle - targ.fr1/fr2 are the inner and outer radii.
136
137 if (targ.fr_pattern .eq. 1) then !square raster
138 t3=grnd()*pi
139 t4=grnd()*pi
140 t5=cos(t3)*targ.fr1
141 t6=cos(t4)*targ.fr2
142 elseif (targ.fr_pattern .eq. 2) then !circular raster
143 t3=grnd()*2.*pi
144 t4=sqrt(grnd())*(targ.fr2-targ.fr1)+targ.fr1
145 t5=cos(t3)*t4
146 t6=sin(t3)*t4
147 else !no raster
148 gaskelld 1.1 t5=0.0
149 t6=0.0
150 endif
151
152 main.target.x = main.target.x+t5
153 main.target.y = main.target.y+t6
154 main.target.z = (0.5-grnd())*targ.length+targ.zoffset
155 main.target.rastery = t6 !'raster' contribution to vert. pos.
156
157 ! Take fluctuations of the beam energy into account, and remember to correct
158 ! for ionization loss in the target and Coulomb acceleration of incoming
159 ! electron. Remove targ.zoffset from the z position of the scattering
160 ! in order to get the position relative to the center of the target.
161
162 call trip_thru_target (1, main.target.z-targ.zoffset, Ebeam, 0.0d0,
163 > main.target.Eloss(1), main.target.teff(1),Me,1)
164 if (.not.using_Eloss) main.target.Eloss(1) = 0.0
165 if (using_Coulomb) then
166 main.target.Coulomb=targ.Coulomb_constant*(3.-(grnd())**(2./3.))
167 else
168 main.target.Coulomb=0.0
169 gaskelld 1.1 endif
170 vertex.Ein = Ebeam + (grnd()-0.5)*dEbeam +
171 > main.target.Coulomb - main.target.Eloss(1)
172
173 ! ... deterimine known variation in Ein from Ebeam_vertex_ave and in
174 ! ... Ee due to Coulomb energy to compare limits to generated event by event.
175 ! ... (USED TO SHIFT 'LIMIT' VALUES IN UPDATE_RANGE CALLS. ARE THEY NEEDED???)
176
177 main.Ein_shift = vertex.Ein - Ebeam_vertex_ave
178 main.Ee_shift = main.target.Coulomb - targ.Coulomb.ave
179
180 ! Initialize success code and fractional weight due to radiation and
181 ! generation tricks
182
183 5 success = .false.
184 main.gen_weight = 1.0
185
186 ! Generated quantities: (phase_space NOT YET IMPLEMENTED).
187 !
188 ! phase_space: Generate electron E,yptar,xptar and hadron yptar,xptar??
189 ! doing_hyd_elast: Generate electron angles. Solve for everything else.
190 gaskelld 1.1 ! doing_deuterium: Generate electron energy and angles, proton angles.
191 ! Solve for proton momentum, p_fermi.
192 ! doing_eep, A>2: generate electron and hadron energy, angles. Solve for Em,Pm.
193 ! doing_pion: generate electron energy and angles, hadron angles, p_fermi, Em.
194 ! Solve for hadron momentum.
195 ! doing_kaon: as doing_pion.
196 ! doing_delta: as doing_pion.
197 ! doing_rho: as doing_pion.
198 ! doing_semi: Generate electron E,yptar,xptar and hadron E, yptar,xptar
199 !
200 ! The above is summarized in the following table:
201 !
202 ! ELECTRON HADRON
203 ! ------------------ ------------------
204 ! E yptar xptar E yptar xptar p_fermi Em
205 !
206 !H(e,e'p) X X
207 !D(e,e'p) X X X X X
208 !A(e,e'p) X X X X X X
209 !----------------------------------------------------------------------
210 !H(e,e'pi) X X X X X
211 gaskelld 1.1 !A(e,e'pi) X X X X X X X
212 !H(e,e'K) X X X X X
213 !A(e,e'K) X X X X X X X
214 !H(e,e'p)pi X X X X X
215 !A(e,e'p)pi X X X X X X X
216 !----------------------------------------------------------------------
217 !phase_space X X X ? X X
218 !
219 ! So our procedure is the following:
220 ! 1) Always generate electron yptar and xptar
221 ! 2) generate hadron yptar and xptar for all cases except H(e,e'p).
222 ! 3) Generate hadron E for A(e,e'p)
223 ! 4) generate electron E for all but hydrogen elastic.
224 ! 5) generate p_fermi, Em for A(e,e'pi) and A(e,e'K).
225 !
226 ! After we generate xptar/yptar/energy, we calculate physics angles (theta/phi),
227 ! momenta, unit vectors, etc... here and/or in complete_ev.
228 !
229 ! Note that there are also jacobians associated with some and/or all of
230 ! the above.
231 ! 1: We generate uniformly in xptar/yptar, not theta/phi. We define the
232 gaskelld 1.1 ! phase space volume (genvol contribution) as the product of the xptar/yptar
233 ! range, and have a jacobian for each event taking into account the mapping
234 ! between the solid angle on the unit sphere, and the dxptar/dyptar volume
235 ! (the jacobian is 1/cos**3(dtheta), where dtheta is the angle between the
236 ! event and the central spectrometer vector
237 ! 2: For the D(e,e'p), we take Em as fixed in order to calculate the proton
238 ! energy. There is a jacobian ( |dEp'/dEm| ). It comes from integrating
239 ! over the energy conservation delta function: delta(E_D - E_p - E_n - Em).
240
241 ! Generate Electron Angles (all cases):
242 vertex.e.yptar=gen.e.yptar.min+grnd()*(gen.e.yptar.max-gen.e.yptar.min)
243 vertex.e.xptar=gen.e.xptar.min+grnd()*(gen.e.xptar.max-gen.e.xptar.min)
244
245 ! Generate Hadron Angles (all but H(e,e'p)):
246 if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon
247 1 .or.doing_delta.or.doing_semi.or.doing_rho) then
248 vertex.p.yptar=gen.p.yptar.min+grnd()*
249 1 (gen.p.yptar.max-gen.p.yptar.min)
250 vertex.p.xptar=gen.p.xptar.min+grnd()*
251 1 (gen.p.xptar.max-gen.p.xptar.min)
252 endif
253 gaskelld 1.1
254 ! Generate Hadron Momentum (A(e,e'p) or semi-inclusive production).
255 if (doing_heavy .or. doing_semi) then
256 Emin = max(gen.p.E.min, gen.sumEgen.min - gen.e.E.max)
257 Emax = min(gen.p.E.max, gen.sumEgen.max - gen.e.E.min)
258 if (Emin.gt.Emax) goto 100
259 main.gen_weight=main.gen_weight*(Emax-Emin)/(gen.p.E.max-gen.p.E.min)
260 vertex.p.E = Emin + grnd()*(Emax-Emin)
261 vertex.p.P = sqrt(vertex.p.E**2 - Mh2)
262 vertex.p.delta = 100.*(vertex.p.P-spec.p.P)/spec.p.P
263 endif
264
265 ! Generate Electron Energy (all but hydrogen elastic)
266 if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon
267 1 .or.doing_delta.or.doing_rho.or.doing_semi) then
268 Emin=gen.e.E.min
269 Emax=gen.e.E.max
270 if (doing_deuterium .or. doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho) then
271 Emin = max(Emin,gen.sumEgen.min)
272 Emax = min(Emax,gen.sumEgen.max)
273 else if (doing_heavy) then ! A(e,e'p)
274 gaskelld 1.1 Emin = max(Emin, gen.sumEgen.min - vertex.p.E)
275 Emax = min(Emax, gen.sumEgen.max - vertex.p.E)
276 endif
277 if (Emin.gt.Emax) goto 100
278 main.gen_weight=main.gen_weight*(Emax-Emin)/(gen.e.E.max-gen.e.E.min)
279 vertex.e.E = Emin + grnd()*(Emax-Emin)
280 vertex.e.P = vertex.e.E
281 vertex.e.delta = 100.*(vertex.e.P-spec.e.P)/spec.e.P
282 endif !not (doing_hyd_elast)
283
284
285 ! Calculate the electron and proton PHYSICS angles from the spectrometer angles.
286 ! Note that the proton angles are not yet know for hydrogen elastic.
287
288 call physics_angles(spec.e.theta,spec.e.phi,
289 & vertex.e.xptar,vertex.e.yptar,vertex.e.theta,vertex.e.phi)
290 call physics_angles(spec.p.theta,spec.p.phi,
291 & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi)
292
293 ! Generate Fermi Momentum and Em for A(e,e'pi) and A(e,e'K).
294 pfer=0.0
295 gaskelld 1.1 pferx=0.0
296 pfery=0.0
297 pferz=0.0
298 vertex.Em=0.0
299 efer=targ.Mtar_struck !used for pion/kaon xsec calcs.
300 if(doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon.or.
301 1 doing_deutdelta.or.doing_hedelta.or.doing_deutrho.or.doing_herho
302 2 .or.doing_deutsemi)then
303 ranprob=grnd()
304 ii=1
305 do while (ranprob.gt.mprob(ii) .and. ii.lt.nump)
306 ii=ii+1
307 enddo
308 if (ii.eq.1) then
309 pferlo=0
310 else
311 pferlo=(pval(ii-1)+pval(ii))/2
312 endif
313 if (ii.eq.nump) then
314 pferhi=pval(nump)
315 else
316 gaskelld 1.1 pferhi=(pval(ii)+pval(ii+1))/2
317 endif
318 pfer=pferlo+(pferhi-pferlo)*grnd()
319 ranth1=grnd()*2.-1.0
320 ranth=acos(ranth1)
321 ranph=grnd()*2.*pi
322 pferx=sin(ranth)*cos(ranph)
323 pfery=sin(ranth)*sin(ranph)
324 pferz=cos(ranth)
325
326 if (doing_deutpi.or.doing_deutkaon.or.doing_deutdelta
327 1 .or.doing_deutrho .or.doing_deutsemi) then !Em = binding energy
328 vertex.Em = Mp + Mn - targ.M
329 m_spec = targ.M - targ.Mtar_struck + vertex.Em != Mn(Mp) for pi+(-)
330 efer = targ.M - sqrt(m_spec**2+pfer**2)
331 endif
332 if (doing_hepi .or. doing_hekaon .or. doing_hedelta .or. doing_herho) then
333 call generate_em(pfer,vertex.Em) !Generate Em
334 m_spec = targ.M - targ.Mtar_struck + vertex.Em != M^*_{A-1}
335 efer = targ.M - sqrt(m_spec**2+pfer**2)
336 endif
337 gaskelld 1.1 endif
338
339 ! Compute all non-generated quantities
340
341 if (debug(5)) write(6,*)'gen: calling comp_ev with false, main, vertex'
342 if (debug(3)) write(6,*)'gen: calling comp_ev with false, main, vertex'
343 if (debug(3)) write(6,*)'gen: Ein, E =',vertex.Ein,vertex.e.E
344 call complete_ev(main,vertex,success)
345
346 main.sigcc = 1.0
347
348 if (debug(2)) write(6,*)'gen: initial success =',success
349 if (.not.success) goto 100
350
351 ! ........ temporary storage of Trec for generated event
352
353 main.Trec = vertex.Trec
354
355 ! Radiate the event, if requested. If we get an event weight of zero (i.e.
356 ! if we CAN'T radiate our event into the acceptance) then success is
357 ! false. generate_rad also set orig kinematics, and cuts on Em/Pm.
358 gaskelld 1.1 ! If not using_rad, must do these here.
359
360 if (using_rad) then
361 call generate_rad(main,vertex,orig,success)
362 if (debug(2)) write(6,*)'gen: after gen_rad, success =',success
363 else
364 success = .true.
365 if (doing_heavy) success =
366 > (vertex.Em .ge. VERTEXedge.Em.min .and.
367 > vertex.Em .le. VERTEXedge.Em.max .and.
368 > vertex.Pm .ge. VERTEXedge.Pm.min .and.
369 > vertex.Pm .le. VERTEXedge.Pm.max)
370 if (success) then
371 do i = 1, neventfields
372 orig.all(i) = vertex.all(i)
373 enddo
374 endif
375 endif
376
377
378 C DJG need to decay the rho here before we begin transporting through the
379 gaskelld 1.1 C DJG spectrometer
380 if(doing_rho) then
381 call rho_decay(orig.p.xptar,orig.p.yptar,orig.p.delta,spec.p.P,Mh2,
382 1 main.epsilon,orig.Q2)
383 endif
384
385 100 if (debug(2)) write(6,*)'gen: final success =',success
386 if (debug(2)) write(6,*)'gen: ending'
387 return
388 end
389
390 !---------------------------------------------------------------------
391
392 subroutine complete_ev(main,vertex,success)
393
394 implicit none
395 include 'simulate.inc'
396
397 real*8 a, b, c, r, t, QA, QB, QC, radical
398 real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z
399 real*8 new_y_x,new_y_y,new_y_z,dummy
400 gaskelld 1.1 real*8 px,py,pz,qx,qy,qz
401 real*8 oop_x,oop_y
402 real*8 krel,krelx,krely,krelz
403 real*8 MM
404
405 real*8 Ehad2,E_rec
406 real*8 W2
407 real*8 grnd !random # generator.
408
409 logical success
410 record /event_main/ main
411 record /event/ vertex
412
413 !-----------------------------------------------------------------------
414 ! Calculate everything left in the /event/ structure, given all necessary
415 ! GENERATION values (some set of xptar,yptar,delta for both arms and p_fermi,
416 ! and p_fermi, depending on the scattering process: see table in generate.f
417 !
418 ! The SINGLE element of /event/ NOT computed here is sigcc.
419 !
420 ! Another small anomaly is that main.jacobian IS computed here.
421 gaskelld 1.1 ! This is because all the necessary terms have to be
422 ! computed here anyway to calculate /event/ qties.
423 !
424 !-----------------------------------------------------------------------
425
426 ! Initialize
427
428 success = .false.
429 main.jacobian = 1.0
430
431 ! ... unit vector components of outgoing e,p
432 ! ... z is DOWNSTREAM, x is DOWN and y is LEFT looking downstream.
433
434 if (debug(4)) write(6,*)'comp_ev: at 1'
435 vertex.ue.x = sin(vertex.e.theta)*cos(vertex.e.phi)
436 vertex.ue.y = sin(vertex.e.theta)*sin(vertex.e.phi)
437 vertex.ue.z = cos(vertex.e.theta)
438 if (.not.doing_hyd_elast) then
439 vertex.up.x = sin(vertex.p.theta)*cos(vertex.p.phi)
440 vertex.up.y = sin(vertex.p.theta)*sin(vertex.p.phi)
441 vertex.up.z = cos(vertex.p.theta)
442 gaskelld 1.1 endif
443 if (debug(4)) write(6,*)'comp_ev: at 2'
444
445 ! First finish off the e side
446 ! Calculate scattered electron energy for hydrogen/deuterium (e,e'p)
447
448 if (doing_hyd_elast) then
449 vertex.e.E = vertex.Ein*Mh/(Mh+vertex.Ein*(1.-vertex.ue.z))
450
451 if (vertex.e.E.gt.vertex.Ein) return
452 vertex.e.P = vertex.e.E
453 vertex.e.delta = (vertex.e.P - spec.e.P)*100./spec.e.P
454 if (debug(4)) write(6,*)'comp_ev: at 3'
455 endif
456
457 ! The q vector
458
459 if (debug(5)) write(6,*)'comp_ev: Ein,E,uez=',vertex.Ein,vertex.e.E,vertex.ue.z
460
461 vertex.nu = vertex.Ein - vertex.e.E
462 vertex.Q2 = 2*vertex.Ein*vertex.e.E*(1.-vertex.ue.z)
463 gaskelld 1.1 vertex.q = sqrt(vertex.Q2 + vertex.nu**2)
464 vertex.xbj = vertex.Q2/2./Mp/vertex.nu
465
466
467 vertex.uq.x = - vertex.e.P*vertex.ue.x / vertex.q
468 vertex.uq.y = - vertex.e.P*vertex.ue.y / vertex.q
469 vertex.uq.z =(vertex.Ein - vertex.e.P*vertex.ue.z)/ vertex.q
470 if (abs(vertex.uq.x**2+vertex.uq.y**2+vertex.uq.z**2-1).gt.0.01)
471 > stop 'Error in q vector normalization'
472 if (debug(4)) write(6,*)'comp_ev: at 5'
473
474
475 ! Now complete the p side, along with vertex.Em, vertex.Pm, vertex.Mrec.
476 ! NOTE: Coherent pion/kaon production (bound final state) is treated as
477 ! hydrogen, but with targ.Mtar_struck=targ.M, targ.Mtar_rec=bound final state.
478
479 if (doing_hyd_elast) then !p = q
480 vertex.Em = 0.0
481 vertex.Pm = 0.0
482 vertex.Mrec = 0.0
483 vertex.up.x = vertex.uq.x
484 gaskelld 1.1 vertex.up.y = vertex.uq.y
485 vertex.up.z = vertex.uq.z
486 vertex.p.P = vertex.q
487 vertex.p.theta = acos(vertex.up.z)
488 if (abs(vertex.up.x/sin(vertex.p.theta)).gt.1)
489 > write(6,*) 'cos(phi)=',vertex.up.x/sin(vertex.p.theta)
490 vertex.p.phi = atan2(vertex.up.y,vertex.up.x)
491 if (vertex.p.phi.lt.0.) vertex.p.phi=vertex.p.phi+2.*pi
492 call spectrometer_angles(spec.p.theta,spec.p.phi,
493 & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi)
494 vertex.p.E = sqrt(vertex.p.P**2+Mh2)
495 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
496 if (debug(4)) write(6,*)'comp_ev: at 6'
497
498 elseif (doing_deuterium) then !need Ep, and a jacobian.
499
500 vertex.Em = targ.Mtar_struck + targ.Mrec - targ.M !=2.2249 MeV
501 vertex.Mrec = targ.M - targ.Mtar_struck + vertex.Em !=targ.Mrec
502
503 a = -1.*vertex.q*(vertex.uq.x*vertex.up.x+vertex.uq.y*vertex.up.y+vertex.uq.z*vertex.up.z)
504 b = vertex.q**2
505 gaskelld 1.1 c = vertex.nu + targ.M
506 t = c**2 - b + Mh2 - vertex.Mrec**2
507 QA = 4.*(a**2 - c**2)
508 QB = 4.*c*t
509 QC = -4.*a**2*Mh2 - t**2
510 radical = QB**2 - 4.*QA*QC
511 if (radical.lt.0) return
512 vertex.p.E = (-QB - sqrt(radical))/2./QA
513
514 ! Check for two solutions
515 ! if ( (-QB + sqrt(radical))/2./QA .gt. Mh ) then
516 ! write(6,*) 'There are two valid solutions for the hadron momentum'
517 ! write(6,*) 'We always pick one, so this may be a problem, and needs to be checked'
518 ! write(6,*) 'solns=',(-QB - sqrt(radical))/2./QA,(-QB + sqrt(radical))/2./QA
519 ! endif
520
521 ! Check for two solutions, but only print warning if both are within
522 ! event generation limits.
523 Ehad2 = (-QB + sqrt(radical))/2./QA
524 if (Ehad2.gt.edge.p.E.min .and. Ehad2.lt.edge.p.E.max .and. ntried.le.5000) then
525 write(6,*) 'The low-momentum solution to E_hadron is within the spectrometer generation'
526 gaskelld 1.1 write(6,*) 'limits. If it is in the acceptance, Its the end of the world as we know it!!!'
527 write(6,*) 'E_hadron solns=',vertex.p.E,Ehad2
528 endif
529
530
531 if (vertex.p.E.le.Mh) return
532 vertex.p.P = sqrt(vertex.p.E**2 - Mh2)
533 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
534
535 ! ........ the Jacobian here is |dEp'/dEm|
536 main.jacobian = (t*(c-vertex.p.E) + 2*c*vertex.p.E*(vertex.p.E-c)) /
537 1 (2*(a**2-c**2)*vertex.p.E + c*t)
538 main.jacobian = abs(main.jacobian)
539
540
541 elseif (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho) then
542
543 if (doing_rho) then
544 Mh = Mrho
545 ! DJG give the rho mass some width (non-relativistic Breit-Wigner)
546 Mh = Mh + 0.5*150.2*tan((2.*grnd()-1.)*atan(2.*500./150.2))
547 gaskelld 1.1 Mh2 = Mh*Mh
548 ntup.rhomass=Mh
549 c write(6,*) 'rho mass is', Mh
550 endif
551
552
553
554 vertex.Pm = pfer !vertex.Em generated at beginning.
555 vertex.Mrec = targ.M - targ.Mtar_struck + vertex.Em
556 a = -1.*vertex.q*(vertex.uq.x*vertex.up.x+vertex.uq.y*vertex.up.y+vertex.uq.z*vertex.up.z)
557 b = vertex.q**2
558 c = vertex.nu + targ.M
559
560 ! For nuclei, correct for fermi motion and missing energy. Also, check
561 ! second solution to quadratic equation - there are often two valid
562 ! solutions, and we always pick the larger one (which is the forward going
563 ! one in the center of mass) and HOPE that the smaller one is never in the
564 ! acceptance. If the low momentum solution IS within the acceptance, we
565 ! have big problems.
566 if (doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon.or.
567 1 doing_deutdelta.or.doing_hedelta.or.doing_deutrho.or.doing_herho) then
568 gaskelld 1.1 a = a - abs(pfer)*(pferx*vertex.up.x+pfery*vertex.up.y+pferz*vertex.up.z)
569 b = b + pfer**2 + 2*vertex.q*abs(pfer)*
570 1 (pferx*vertex.uq.x+pfery*vertex.uq.y+pferz*vertex.uq.z)
571 *** c = c - sqrt(vertex.Mrec**2+pfer**2)
572 c = vertex.nu + efer !same as above, but this way if we redefine
573 !'efer', it's the same everywhere.
574 endif
575 t = c**2 - b + Mh2 - targ.Mrec_struck**2
576 QA = 4.*(a**2 - c**2)
577 QB = 4.*c*t
578 QC = -4.*a**2*Mh2 - t**2
579
580 ! write(6,*) ' '
581 ! write(6,*) ' '
582 ! write(6,*) 'E0=',vertex.Ein
583 ! write(6,*) 'P_elec,P_prot=',vertex.e.P/1000.,vertex.p.P/1000.
584 ! write(6,*) 'thetae,phie=',vertex.e.theta*180./pi,vertex.e.phi*180./pi
585 ! write(6,*) 'thetap,phip=',vertex.p.theta*180./pi,vertex.p.phi*180./pi
586 ! write(6,*) 'q,nu,costhetapq=',vertex.q,vertex.nu,(vertex.uq.x*vertex.up.x+vertex.uq.y*vertex.up.y+vertex.uq.z*vertex.up.z)
587 ! write(6,*) 'a,b,c=',a/1000.,b/1000000.,c/1000.
588 ! write(6,*) 't=',t/1000000.
589 gaskelld 1.1 ! write(6,*) 'A,B,C=',QA/1.d6,QB/1.d9,QC/1.d12
590 ! write(6,*) 'rad=',QB**2 - 4.*QA*QC
591 ! write(6,*) 'e1,e2=',(-QB-sqrt(radical))/2000./QA,(-QB+sqrt(radical))/2000./QA
592 ! write(6,*) 'E_pi1,2=',vertex.nu+targ.M-(-QB-sqrt(radical))/2./QA,
593 ! > vertex.nu+targ.M-(-QB+sqrt(radical))/2./QA
594
595
596 radical = QB**2 - 4.*QA*QC
597 if (radical.lt.0) return
598 vertex.p.E = (-QB - sqrt(radical))/2./QA
599 Ehad2 = (-QB + sqrt(radical))/2./QA
600
601 if (doing_delta.or.doing_rho) then !choose one of the two solutions.
602 ! write(6,*) ' e1, e2=',vertex.p.E,Ehad2
603 if (grnd().gt.0.5) vertex.p.E = Ehad2
604 else !verify that 'backwards' soln. is no good.
605 if (Ehad2.gt.edge.p.E.min .and. Ehad2.lt.edge.p.E.max .and. ntried.le.5000) then
606 write(6,*) 'The low-momentum solution to E_hadron is within the spectrometer generation'
607 write(6,*) 'limits. If it is in the acceptance, Its the end of the world as we know it!!!'
608 write(6,*) 'E_hadron solns=',vertex.p.E,Ehad2
609 endif
610 gaskelld 1.1 endif
611
612 E_rec=c-vertex.p.E !energy of recoil system
613 if (E_rec.le.targ.Mrec_struck) return !non-physical solution
614
615 if (vertex.p.E.le.Mh) return
616 vertex.p.P = sqrt(vertex.p.E**2 - Mh2)
617 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
618 ! write(6,*) 'p,e=',vertex.p.P,vertex.p.E
619
620 elseif (doing_phsp) then
621
622 vertex.p.P = spec.p.P !????? single arm phsp??
623 vertex.p.E= sqrt(Mh2+vertex.p.P**2)
624 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
625 if (debug(4)) write(6,*)'comp_ev: at 7.5',Mh2,vertex.p.E
626
627 endif
628
629 if (debug(4)) write(6,*)'comp_ev: at 7'
630
631 gaskelld 1.1 ! Compute some pion and kaon stuff. Some of these should be OK for proton too.
632
633 if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then
634 W2 = targ.M**2 + 2.*targ.M*vertex.nu - vertex.Q2
635 main.W = sqrt(abs(W2)) * W2/abs(W2)
636 main.epsilon=1./(1. + 2.*(1+vertex.nu**2/vertex.Q2)*tan(vertex.e.theta/2.)**2)
637 main.theta_pq=acos(vertex.up.x*vertex.uq.x+vertex.up.y*vertex.uq.y+vertex.up.z*vertex.uq.z)
638 main.t = vertex.Q2 - Mh2 + 2*vertex.nu*vertex.p.E -
639 > 2*vertex.p.P*vertex.q*cos(main.theta_pq)
640 main.tmin = vertex.Q2 - Mh2 + 2*vertex.p.E*vertex.nu -
641 > 2*vertex.p.P*vertex.q
642 main.q2 = vertex.Q2
643
644 ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE.
645 ! Therefore, define a new system with the z axis parallel to q, and
646 ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q
647 ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane
648 ! particles closer to the downstream beamline than q.
649 ! phi=90 is above the horizontal plane when q points to the right, and
650 ! below the horizontal plane when q points to the left.
651 ! Also take into account the different definitions of x, y and z in
652 gaskelld 1.1 ! replay and SIMC:
653 ! As seen looking downstream: replay SIMC (old_simc)
654 ! x right down (left)
655 ! y down left (up)
656 ! z all have z pointing downstream
657 !
658 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
659
660 qx = -vertex.uq.y !convert to 'replay' coord. system
661 qy = vertex.uq.x
662 qz = vertex.uq.z
663 px = -vertex.up.y
664 py = vertex.up.x
665 pz = vertex.up.z
666
667 dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
668 new_x_x = -qx*qz/dummy
669 new_x_y = -qy*qz/dummy
670 new_x_z = (qx**2 + qy**2)/dummy
671
672 dummy = sqrt(qx**2 + qy**2)
673 gaskelld 1.1 new_y_x = qy/dummy
674 new_y_y = -qx/dummy
675 new_y_z = 0.0
676
677 p_new_x = px*new_x_x + py*new_x_y + pz*new_x_z
678 p_new_y = px*new_y_x + py*new_y_y + pz*new_y_z
679
680 main.phi_pq = atan2(p_new_y,p_new_x) !atan2(y,x)=atan(y/x)
681 if (main.phi_pq.lt.0.d0) main.phi_pq=main.phi_pq+2.*pi
682 ! if (p_new_y.lt.0.) then
683 ! main.phi_pq = 2*pi - main.phi_pq
684 ! endif
685
686 ! if ((p_new_x**2+p_new_y**2).eq.0.) then
687 ! main.phi_pq = 0.0
688 ! else
689 ! main.phi_pq = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2))
690 ! endif
691 ! if (p_new_y.lt.0.) then
692 ! main.phi_pq = 2*pi - main.phi_pq
693 ! endif
694 gaskelld 1.1
695 if (debug(2)) then
696 write(6,*)'comp_ev: nu =',vertex.nu
697 write(6,*)'comp_ev: Q2 =',vertex.Q2
698 write(6,*)'comp_ev: theta_e =',vertex.e.theta
699 write(6,*)'comp_ev: epsilon =',main.epsilon
700 write(6,*)'comp_ev: theta_pq =',main.theta_pq
701 write(6,*)'comp_ev: phi_pq =',main.phi_pq
702 write(6,*)'comp_ev: E_hadron =',vertex.p.E
703 endif
704 endif !end of pion/kaon specific stuff.
705
706 if (debug(4)) write(6,*)'comp_ev: at 8'
707
708 ! Compute the Pm vector in in SIMC LAB system, with x down, and y to the left.
709 ! Computer Parallel, Perpendicular, and Out of Plane componenants.
710 ! Parallel is component along q_hat. Out/plane is component along
711 ! (e_hat) x (q_hat) (oop_x,oop_y are components of out/plane unit vector)
712 ! Perp. component is what's left: along (q_hat) x (oop_hat).
713 ! So looking along q, out of plane is down, perp. is left.
714
715 gaskelld 1.1 vertex.Pmx = vertex.p.P*vertex.up.x - vertex.q*vertex.uq.x
716 vertex.Pmy = vertex.p.P*vertex.up.y - vertex.q*vertex.uq.y
717 vertex.Pmz = vertex.p.P*vertex.up.z - vertex.q*vertex.uq.z
718 vertex.Pmiss = sqrt(vertex.Pmx**2+vertex.Pmy**2+vertex.Pmz**2)
719 vertex.Emiss = vertex.nu + targ.M - vertex.p.E
720
721 !STILL NEED SIGN FOR PmPer!!!!!!
722
723 oop_x = -vertex.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat)
724 oop_y = vertex.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat)
725 vertex.PmPar = (vertex.Pmx*vertex.uq.x + vertex.Pmy*vertex.uq.y + vertex.Pmz*vertex.uq.z)
726 vertex.PmOop = (vertex.Pmx*oop_x + vertex.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2)
727 vertex.PmPer = sqrt( max(0.d0, vertex.Pm**2 - vertex.PmPar**2 - vertex.PmOop**2 ) )
728
729 if (debug(4)) write(6,*)'comp_ev: at 9',vertex.Pmx,vertex.Pmy,vertex.Pmz
730
731 ! Calculate Em, Pm, Mrec, Trec for all cases not already done.
732 ! For doing_heavy, get Mrec from nu+M=Ep+Erec, and Erec**2=Mrec**2+Pm**2
733 ! For (e,e'pi/K), could go back and determine momentum of recoil struck
734 ! particle, and get Mrec and Trec seperately for struck nucleon(hyperon)
735 ! and A-1 system. For now, just take Trec for the A-1 system, and ignore
736 gaskelld 1.1 ! the recoiling struck nucleon (hyperon), so Trec=0 for hydrogen target.
737
738 if (doing_hyd_elast) then
739 vertex.Trec = 0.0
740 else if (doing_deuterium) then
741 vertex.Pm = vertex.Pmiss
742 vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec
743 else if (doing_heavy) then
744 vertex.Pm = vertex.Pmiss
745 vertex.Mrec = sqrt(vertex.Emiss**2-vertex.Pmiss**2)
746 vertex.Em = targ.Mtar_struck + vertex.Mrec - targ.M
747 vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec
748 else if (doing_hydpi .or. doing_hydkaon .or. doing_hyddelta .or. doing_hydrho) then
749 vertex.Trec = 0.0
750 else if (doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon
751 1 .or.doing_deutdelta.or.doing_hedelta.or.doing_deutrho.or.doing_herho) then
752 vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec
753 else if (doing_semi) then
754 vertex.Pm = vertex.Pmiss
755 vertex.Em = vertex.Emiss
756 endif
757 gaskelld 1.1
758 if (debug(5)) write(6,*) 'vertex.Pm,vertex.Trec,vertex.Em',vertex.Pm,vertex.Trec,vertex.Em
759 if (debug(4)) write(6,*)'comp_ev: at 10'
760
761
762 ! calculate krel for deuteron/heavy pion(kaon). Deuteron is straightforward.
763 ! A>2 case is some approximation for 3He (DJG).
764
765 if (doing_deutpi .or. doing_deutkaon .or. doing_deutdelta .or. doing_deutrho) then
766 if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) write(6,*) 'BAD MM!!!!! Emiss,Pmiss=',vertex.Emiss, vertex.Pmiss
767 MM = sqrt(max(0.d0,vertex.Emiss**2-vertex.Pmiss**2))
768 krel = sqrt( max(0.d0,MM**2-4.*targ.Mrec_struck**2) )
769 else if (doing_hepi .or. doing_hekaon .or. doing_hedelta .or. doing_herho) then
770 if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) write(6,*) 'BAD MM!!!!! Emiss,Pmiss=',vertex.Emiss, vertex.Pmiss
771 MM = sqrt(max(0.d0,vertex.Emiss**2-vertex.Pmiss**2))
772 krelx = vertex.Pmx + 1.5*pferx*pfer
773 krely = vertex.Pmy + 1.5*pfery*pfer
774 krelz = vertex.Pmz + 1.5*pferz*pfer
775 krel = sqrt(krelx**2+krely**2+krelz**2)
776 if(vertex.Em.lt.6.0) krel = -krel !bound state test???
777 endif
778 gaskelld 1.1 ntup.krel = krel
779
780 if(doing_semi) then
781 CDJG if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) then
782 CDJG I should be testing that the missing mass is above two pion
783 CDJG threshold! Otherwise, it's just exclusive
784 if ((vertex.Emiss**2-vertex.Pmiss**2).lt.(Mp+Mpi0)**2) then
785 success=.false.
786 return
787 endif
788 endif
789
790 if(doing_semi) then
791 vertex.zhad = vertex.p.E/vertex.nu
792 vertex.pt2 = vertex.p.P**2*(1.0-cos(main.theta_pq)**2)
793 if(vertex.zhad.gt.1.0) then
794 success=.false.
795 return
796 endif
797 endif
798
799 gaskelld 1.1
800 ! Determine PHYSICS scattering angles theta/phi for the two spectrometer
801 ! vectors, and the Jacobian which we must include in our xsec computation
802 ! to take into account the fact that we're not generating in the physics
803 ! solid angle d(omega) but in 'spectrometer angles' yptar,xptar.
804 ! NOTE: the conversion from spectrometer to physics angles was done when
805 ! the angles were first generated (except for the proton angle for hydrogen
806 ! elastic, where it was done earlier in complete_ev).
807 !
808 ! call physics_angles(spec.e.theta,spec.e.phi,
809 ! & vertex.e.xptar,vertex.e.yptar,vertex.e.theta,vertex.e.phi)
810 ! call physics_angles(spec.p.theta,spec.p.phi,
811 ! & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi)
812
813 r = sqrt(1.+vertex.e.yptar**2+vertex.e.xptar**2) !r=cos(theta-theta0)
814 main.jacobian = main.jacobian / r**3 !1/cos**3(theta-theta0)
815
816 if (doing_heavy .or. doing_pion .or. doing_kaon .or.
817 1 doing_rho .or. doing_delta .or. doing_semi) then
818 r = sqrt(1.+vertex.p.yptar**2+vertex.p.xptar**2)
819 main.jacobian = main.jacobian / r**3 !1/cos**3(theta-theta0)
820 gaskelld 1.1 endif
821
822 if (debug(4)) write(6,*)'comp_ev: at 11'
823
824 ! The effective target thickness that the scattered particles see, and the
825 ! resulting energy losses
826
827 call trip_thru_target (2, main.target.z-targ.zoffset, vertex.e.E,
828 > vertex.e.theta, main.target.Eloss(2), main.target.teff(2),Me,1)
829 call trip_thru_target (3, main.target.z-targ.zoffset, vertex.p.E,
830 > vertex.p.theta, main.target.Eloss(3), main.target.teff(3),Mh,1)
831 if (.not.using_Eloss) then
832 main.target.Eloss(2) = 0.0
833 main.target.Eloss(3) = 0.0
834 endif
835 if (debug(4)) write(6,*)'comp_ev: at 12'
836
837 ! Initialize parameters necessary for radiative corrections
838
839 call radc_init_ev(main,vertex)
840
841 gaskelld 1.1 ! Success code
842
843 success = .true.
844 if(debug(2)) write(6,*)'comp_ev: at end...'
845 return
846 end
847
848 !---------------------------------------------------------------------
849
850 subroutine complete_recon_ev(recon,success)
851
852 implicit none
853 include 'simulate.inc'
854
855 real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z
856 real*8 new_y_x,new_y_y,new_y_z,dummy
857 real*8 px,py,pz,qx,qy,qz
858 real*8 W2
859 real*8 oop_x,oop_y
860 real*8 mm,mmA,mm2,mmA2,t
861
862 gaskelld 1.1 logical success
863 record /event/ recon
864
865 !-----------------------------------------------------------------------
866 ! Calculate everything left in the /event/ structure, given
867 ! recon.Ein, and the following for both recon.e AND recon.p:
868 ! delta,xptar,yptar, z,P,E,theta,phi (with E,P corrected for eloss, if desired).
869 !
870 ! The SINGLE element of /event/ NOT computed here is sigcc (for (e,e'p)).
871 !
872 !-----------------------------------------------------------------------
873
874 ! Initialize
875
876 success = .false.
877
878 recon.Ein = Ebeam_vertex_ave !lowered by most probable eloss (init.f)
879
880 ! ... unit vector components of outgoing e,p
881 ! ... z is DOWNSTREAM, x is DOWN and y is LEFT looking downstream.
882
883 gaskelld 1.1 if (debug(4)) write(6,*)'comp_rec_ev: at 1'
884 recon.ue.x = sin(recon.e.theta)*cos(recon.e.phi)
885 recon.ue.y = sin(recon.e.theta)*sin(recon.e.phi)
886 recon.ue.z = cos(recon.e.theta)
887 recon.up.x = sin(recon.p.theta)*cos(recon.p.phi)
888 recon.up.y = sin(recon.p.theta)*sin(recon.p.phi)
889 recon.up.z = cos(recon.p.theta)
890 if (debug(4)) write(6,*)'comp_rec_ev: at 2'
891
892 ! The q vector
893
894 if (debug(5)) write(6,*)'comp_rec_ev: Ein,E,uez=',recon.Ein,recon.e.E,recon.ue.z
895 recon.nu = recon.Ein - recon.e.E
896 recon.Q2 = 2*recon.Ein*recon.e.E*(1-recon.ue.z)
897 recon.q = sqrt(recon.Q2 + recon.nu**2)
898 recon.uq.x = - recon.e.P*recon.ue.x / recon.q
899 recon.uq.y = - recon.e.P*recon.ue.y / recon.q
900 recon.uq.z =(recon.Ein - recon.e.P*recon.ue.z)/ recon.q
901 W2 = targ.M**2 + 2.*targ.M*recon.nu - recon.Q2
902 recon.W = sqrt(abs(W2)) * W2/abs(W2)
903 recon.xbj = recon.Q2/2./Mp/recon.nu
904 gaskelld 1.1 if (debug(4)) write(6,*)'comp_rec_ev: at 5'
905
906 ! Now complete the p side
907
908 if(doing_phsp)then
909 recon.p.P=spec.p.P
910 recon.p.E=sqrt(Mh2+recon.p.P**2)
911 if (debug(4)) write(6,*)'comp_rec_ev: at 7.5',Mh2,recon.p.E
912 endif
913
914 ! Compute some pion/kaon stuff.
915
916 recon.epsilon=1./(1. + 2.*(1+recon.nu**2/recon.Q2)*tan(recon.e.theta/2.)**2)
917 recon.theta_pq=acos(min(1.0,recon.up.x*recon.uq.x+recon.up.y*recon.uq.y+recon.up.z*recon.uq.z))
918
919 ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE.
920 ! Therefore, define a new system with the z axis parallel to q, and
921 ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q
922 ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane
923 ! particles closer to the downstream beamline than q.
924 ! phi=90 is above the horizontal plane when q points to the right, and
925 gaskelld 1.1 ! below the horizontal plane when q points to the left.
926 ! Also take into account the different definitions of x, y and z in
927 ! replay and SIMC:
928 ! As seen looking downstream: replay SIMC (old_simc)
929 ! x right down (left)
930 ! y down left (up)
931 ! z all have z pointing downstream
932 !
933 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
934
935 qx = -recon.uq.y !convert to 'replay' coord. system
936 qy = recon.uq.x
937 qz = recon.uq.z
938 px = -recon.up.y
939 py = recon.up.x
940 pz = recon.up.z
941
942 dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
943 new_x_x = -qx*qz/dummy
944 new_x_y = -qy*qz/dummy
945 new_x_z = (qx**2 + qy**2)/dummy
946 gaskelld 1.1
947 dummy = sqrt(qx**2 + qy**2)
948 new_y_x = qy/dummy
949 new_y_y = -qx/dummy
950 new_y_z = 0.0
951
952 p_new_x = px*new_x_x + py*new_x_y + pz*new_x_z
953 p_new_y = px*new_y_x + py*new_y_y + pz*new_y_z
954
955 if ((p_new_x**2+p_new_y**2).eq.0.) then
956 recon.phi_pq = 0.0
957 else
958 recon.phi_pq = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2))
959 endif
960 if (p_new_y.lt.0.) then
961 recon.phi_pq = 2*pi - recon.phi_pq
962 endif
963
964 if (debug(2)) then
965 write(6,*)'comp_rec_ev: nu =',recon.nu
966 write(6,*)'comp_rec_ev: Q2 =',recon.Q2
967 gaskelld 1.1 write(6,*)'comp_rec_ev: theta_e =',recon.e.theta
968 write(6,*)'comp_rec_ev: epsilon =',recon.epsilon
969 write(6,*)'comp_rec_ev: theta_pq =',recon.theta_pq
970 write(6,*)'comp_rec_ev: phi_pq =',recon.phi_pq
971 endif
972
973 if (debug(4)) write(6,*)'comp_rec_ev: at 8'
974
975 ! Compute the Pm vector in in SIMC LAB system, with x down, and y to the left.
976 ! Computer Parallel, Perpendicular, and Out of Plane componenants.
977 ! Parallel is component along q_hat. Out/plane is component along
978 ! (e_hat) x (q_hat) (oop_x,oop_y are components of out/plane unit vector)
979 ! Perp. component is what's left: along (q_hat) x (oop_hat).
980 ! So looking along q, out of plane is down, perp. is left.
981
982 recon.Pmx = recon.p.P*recon.up.x - recon.q*recon.uq.x
983 recon.Pmy = recon.p.P*recon.up.y - recon.q*recon.uq.y
984 recon.Pmz = recon.p.P*recon.up.z - recon.q*recon.uq.z
985 recon.Pm = sqrt(recon.Pmx**2+recon.Pmy**2+recon.Pmz**2)
986
987 !STILL NEED SIGN FOR PmPer!!!!!!
988 gaskelld 1.1
989 oop_x = -recon.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat)
990 oop_y = recon.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat)
991 recon.PmPar = (recon.Pmx*recon.uq.x + recon.Pmy*recon.uq.y + recon.Pmz*recon.uq.z)
992 recon.PmOop = (recon.Pmx*oop_x + recon.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2)
993 recon.PmPer = sqrt( max(0.d0, recon.Pm**2 - recon.PmPar**2 - recon.PmOop**2 ) )
994
995 if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then
996 recon.Em = recon.nu + targ.Mtar_struck - recon.p.E
997 mm2 = recon.Em**2 - recon.Pm**2
998 mm = sqrt(abs(mm2)) * abs(mm2)/mm2
999 mmA2= (recon.nu + targ.M - recon.p.E)**2 - recon.Pm**2
1000 mmA = sqrt(abs(mmA2)) * abs(mmA2)/mmA2
1001 t = recon.Q2 - Mh2
1002 > + 2*(recon.nu*recon.p.E - recon.p.P*recon.q*cos(recon.theta_pq))
1003 ntup.mm = mm
1004 ntup.mmA = mmA
1005 ntup.t = t
1006 endif
1007
1008 if(doing_semi.or.doing_rho) then
1009 gaskelld 1.1 recon.zhad = recon.p.E/recon.nu
1010 recon.pt2 = recon.p.P**2*(1.0-cos(recon.theta_pq)**2)
1011 endif
1012
1013 ! Calculate Trec, Em. Trec for (A-1) system (eep), or for struck nucleon (pi/K)
1014 ! Note that there are other ways to calculate 'Em' for the pion/kaon case.
1015 ! This Em for pi/kaon gives roughly E_m + E_rec + T_{A-1} (I think)
1016 if (doing_hyd_elast) then
1017 recon.Trec = 0.0
1018 recon.Em = recon.nu + targ.M - recon.p.E - recon.Trec
1019 else if (doing_deuterium .or. doing_heavy) then
1020 recon.Trec = sqrt(recon.Pm**2+targ.Mrec**2) - targ.Mrec
1021 recon.Em = recon.nu + targ.Mtar_struck - recon.p.E - recon.Trec
1022 else if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho) then
1023 recon.Em = recon.nu + targ.Mtar_struck - recon.p.E
1024 endif
1025
1026 if (debug(5)) write(6,*) 'recon.Pm,recon.Trec,recon.Em',recon.Pm,recon.Trec,recon.Em
1027 if (debug(4)) write(6,*)'comp_rec_ev: at 10'
1028
1029 success=.true.
1030 gaskelld 1.1 return
1031 end
1032
1033 !------------------------------------------------------------------------
1034
1035 subroutine complete_main(force_sigcc,main,vertex,vertex0,recon,success)
1036
1037 implicit none
1038 include 'simulate.inc'
1039
1040 integer i, iPm1
1041 real*8 a, b, r, frac, peepi, peeK, peedelta, peerho, peepiX
1042 real*8 survivalprob
1043 real*8 weight, width, sigep, deForest, tgtweight
1044 logical force_sigcc, success
1045 record /event_main/ main
1046 record /event/ vertex, vertex0, recon
1047
1048 !-----------------------------------------------------------------------
1049 ! Calculate everything left in the /main/ structure that hasn't been
1050 ! required up til now. This routine is called ONLY after we
1051 gaskelld 1.1 ! know an event IS going to contribute, don't want to be doing these
1052 ! calculations if they're not needed!
1053 !
1054 ! A small anomaly is that sigcc from the /event/ structure is computed
1055 ! here since it's not needed til now, and was not calculated by
1056 ! COMPLETE_ev.
1057 !-----------------------------------------------------------------------
1058
1059
1060 ! Initialize success code
1061
1062 if (debug(2)) write(6,*)'comp_main: entering...'
1063 success = .false.
1064
1065 ! The spectral function weighting
1066
1067 if (doing_hyd_elast.or.doing_pion.or.doing_kaon.or.doing_delta.or.doing_phsp.or.doing_rho.or.doing_semi) then !no SF.
1068 main.SF_weight=1.0
1069 else if (doing_deuterium .or. doing_heavy) then
1070 main.SF_weight = 0.0
1071 do i=1,nrhoPm
1072 gaskelld 1.1 weight = 0.0
1073
1074 ! ... use linear interpolation to determine rho(pm)
1075
1076 r = (vertex.Pm-Pm_theory(i).min)/Pm_theory(i).bin
1077 if (r.ge.0 .and. r.le.Pm_theory(i).n) then
1078 iPm1 = nint(r)
1079 if (iPm1.eq.0) iPm1 = 1
1080 if (iPm1.eq.Pm_theory(i).n) iPm1 = Pm_theory(i).n-1
1081 frac = r+0.5-float(iPm1)
1082 b = theory(i,iPm1)
1083 a = theory(i,iPm1+1)-b
1084 weight = a*frac + b
1085 endif
1086
1087 if (doing_heavy) then
1088
1089 width=Emsig_theory(i)/2.0
1090 if (vertex.Em.lt.E_Fermi) weight = 0.0
1091 weight=weight/pi/Em_int_theory(i) * width/
1092 > ((vertex.Em-Em_theory(i))**2+width**2)
1093 gaskelld 1.1 endif
1094 main.SF_weight = main.SF_weight+weight*nprot_theory(i)
1095 enddo ! <i=1,nrhoPm>
1096 endif
1097
1098 ! ... if we have come up with weight<=, quit now (avoid weight<0 in ntuple)
1099 ! ... unless FORCE_SIGCC is on (to get central sigcc).
1100
1101 if (debug(5))write(6,*)'comp_main: calculating cross section'
1102 if (main.SF_weight.le.0 .and. .not.force_sigcc) return
1103
1104 ! ... Cross section, for BOTH vertex and recon (where possible)
1105 ! ... Note that all of these are cross section per nucleon. For A(e,e'p)
1106 ! ... this becomes cross section per nucleus because of the weighting of
1107 ! ... the spectral function. For pion/kaon production, we explicitly
1108 ! ... apply a weight for the number of nucleons involved (protons for pi+ or Y0,
1109 ! ... neutrons for pi- or Y-) (Y=hyperon)
1110
1111 tgtweight = 1.0
1112
1113 if (doing_phsp) then
1114 gaskelld 1.1 main.sigcc = 1.0
1115 main.sigcc_recon = 1.0
1116
1117 elseif (doing_hyd_elast) then
1118 main.sigcc = sigep(vertex)
1119 main.sigcc_recon = sigep(recon)
1120
1121 elseif (doing_deuterium.or.doing_heavy) then
1122 main.sigcc = deForest(vertex)
1123 main.sigcc_recon = deForest(recon)
1124
1125 elseif (doing_pion) then
1126 main.sigcc = peepi(vertex,main)
1127 main.sigcc_recon = 1.0
1128 if (which_pion.eq.1 .or. which_pion.eq.11) then !OK for coherent???
1129 tgtweight = targ.N
1130 else
1131 tgtweight = targ.Z
1132 endif
1133
1134 elseif (doing_kaon) then
1135 gaskelld 1.1 main.sigcc = peeK(vertex,main,survivalprob)
1136 main.sigcc_recon = 1.0
1137 if (which_kaon.eq.2 .or. which_kaon.eq.12) then !OK for coherent???
1138 tgtweight = targ.N
1139 else
1140 tgtweight = targ.Z
1141 endif
1142
1143 elseif (doing_delta) then
1144 main.sigcc = peedelta(vertex,main) !Need new xsec model.
1145 main.sigcc_recon = 1.0
1146
1147 elseif (doing_rho) then
1148 main.sigcc = peerho(vertex,main)
1149 main.sigcc_recon = 1.0
1150
1151 elseif (doing_semi) then
1152 main.sigcc = peepiX(vertex,vertex0,main,survivalprob,.FALSE.)
1153 main.sigcc_recon = 1.0
1154 main.sigcent = peepiX(vertex,vertex0,main,survivalprob,.TRUE.)
1155
1156 gaskelld 1.1 else
1157 main.sigcc = 1.0
1158 main.sigcc_recon = 1.0
1159 endif
1160
1161 if (debug(3)) then
1162 write(6,*)'======================================'
1163 write(6,*)'complete_main:'
1164 write(6,*)'theta_e =',vertex.e.theta
1165 write(6,*)'q mag =',vertex.q
1166 write(6,*)'nu =',vertex.nu
1167 write(6,*)'Q2 =',vertex.Q2/1000000.
1168 write(6,*)'Ein =',vertex.Ein
1169 write(6,*)'hadron mom =',vertex.p.P
1170 write(6,*)'e mom =',vertex.e.P
1171 write(6,*)'mass =',Mh
1172 write(6,*)'epsilon =',main.epsilon
1173 write(6,*)'phi_pq =',main.phi_pq
1174 write(6,*)'theta_pq =',main.theta_pq
1175 write(6,*)'======================================'
1176 endif
1177 gaskelld 1.1
1178 ! The total contributing weight from this event -- this weight is
1179 ! proportional to # experimental counts represented by the event.
1180 ! Apply survival probability to kaons if we're not modeling decay.
1181
1182 main.weight = main.SF_weight*main.jacobian*main.gen_weight*main.sigcc
1183 main.weight = main.weight * tgtweight !correct for #/nucleons involved
1184 if ((doing_kaon.or.doing_semika) .and. .not.doing_decay)
1185 > main.weight = main.weight*survivalprob
1186 if (debug(5))write(6,*) 'gen_weight = ',main.gen_weight,
1187 > main.jacobian,main.sigcc
1188
1189 success = .true.
1190
1191 if (debug(2)) write(6,*)'comp_main: ending, success =',success
1192 return
1193 end
1194
1195
1196 subroutine physics_angles(theta0,phi0,dx,dy,theta,phi)
1197
1198 gaskelld 1.1 !Generate physics angles in lab frame. Theta is angle from beamline.
1199 !phi is projection on x-y plane (so phi=0 is down, sos=pi/2, hms=3*pi/2.
1200 !
1201 !theta=acos( (cos(theta0)-dy*sin(theta0)*sin(phi0))/ sqrt(1+dx**2+dy**2) )
1202 !phi=atan( (dy*cos(theta0)+sin(theta0)*sin(phi0)) / (sin(theta0)*cos(phi0)+dx) )
1203 !
1204 ! Note that these formulae assume phi0=pi/2 or 3*pi/2 (thus the sin(phi0)
1205 ! gives a -/+ sign for the HMS/SOS). Thus, we set the cos(phi0) term to zero.
1206
1207 real*8 dx,dy !dx/dy (xptar/yptar) for event.
1208 real*8 theta0,phi0 !central physics angles of spectrometer.
1209 real*8 theta,phi !physics angles for event.
1210 real*8 r,sinth,costh,sinph,cosph !intermediate variables.
1211 real*8 tmp
1212
1213 include 'constants.inc'
1214
1215 costh = cos(theta0)
1216 sinth = sin(theta0)
1217 sinph = sin(phi0)
1218 cosph = cos(phi0)
1219 gaskelld 1.1 r = sqrt( 1. + dx**2 + dy**2 )
1220
1221 if (abs(cosph).gt.0.0001) then !phi not at +/- pi/2
1222 write(6,*) 'theta,phi will be incorrect if phi0 <> pi/2 or 3*pi/2'
1223 write(6,*) 'phi0=',phi0,'=',phi0*180/pi,'degrees'
1224 endif
1225
1226 tmp=(costh - dy*sinth*sinph) / r
1227 if (abs(tmp).gt.1) write(6,*) 'tmp=',tmp
1228 theta = acos( (costh - dy*sinth*sinph) / r )
1229 if (dx.ne.0.0) then
1230 phi = atan( (dy*costh + sinth*sinph) / dx ) !gives -90 to 90 deg.
1231 if (phi.le.0) phi=phi+pi !make 0 to 180 deg.
1232 if (sinph.lt.0.) phi=phi+pi !add pi to phi for HMS
1233 else
1234 phi = phi0
1235 endif
1236
1237 return
1238 end
1239
1240 gaskelld 1.1
1241
1242 subroutine spectrometer_angles(theta0,phi0,dx,dy,theta,phi)
1243
1244 !Generate spectrometer angles from physics angles in lab frame.
1245 !Theta is angle from beamline.
1246 !phi is projection on x-y plane (so phi=0 is down, sos=pi/2, hms=3*pi/2.
1247
1248 real*8 dx,dy !dx/dy (xptar/yptar) for event.
1249 real*8 theta0,phi0 !central physics angles of spectrometer.
1250 real*8 theta,phi !physics angles for event.
1251 real*8 x,y,z !intermediate variables.
1252 real*8 x0,y0,z0 !intermediate variables.
1253 real*8 cos_dtheta,y_event
1254
1255 include 'constants.inc'
1256
1257 x = sin(theta)*cos(phi)
1258 y = sin(theta)*sin(phi)
1259 z = cos(theta)
1260 x0 = sin(theta0)*cos(phi0)
1261 gaskelld 1.1 y0 = sin(theta0)*sin(phi0)
1262 z0 = cos(theta0)
1263
1264 cos_dtheta = x*x0 + y*y0 + z*z0
1265 dx = x / cos_dtheta
1266 dy = sqrt(1/cos_dtheta**2-1.-dx**2)
1267
1268 y_event = y/cos_dtheta !projected to plane perp. to spectrometer.
1269 if (y_event .lt. y0) dy = -dy
1270
1271 return
1272 end
|
Return to
event.f
CVS log
Up to [HallC] / simc_semi