1 jones 1.1 subroutine limits_update(main,vertex,orig,recon,doing_deuterium,
2 > doing_pion,doing_kaon,doing_dvcs,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_dvcs
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 jones 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_dvcs) 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 jones 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 jones 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 jones 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 jones 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 gauss1
117 logical success
118 real*8 grnd !random # generator.
119 record /event_main/ main
120 record /event/ vertex, orig
121
122 real*8 nsig_max
123 parameter(nsig_max=3.0d0) !max #/sigma for gaussian ran #s.
124
125
126 ! Randomize the position of the interaction inside the available region.
127 jones 1.1 ! gen.xwid and gen.ywid are the intrinsic beam widths (one sigma value).
128 ! Use a gaussian beam distribution with +/- 3.0 sigma (add raster afterwards).
129 if (debug(2)) write(6,*)'gen: entering'
130 main.target.x = gauss1(nsig_max)*gen.xwid+targ.xoffset
131 main.target.y = gauss1(nsig_max)*gen.ywid+targ.yoffset
132
133 ! fr_pattern=1 - rectangle - targ.fr1/fr2 are the x/y raster half-widths.
134 ! fr_pattern=2 - circle - targ.fr1/fr2 are the inner and outer radii.
135
136 if (targ.fr_pattern .eq. 1) then !square raster
137 t3=grnd()*pi
138 t4=grnd()*pi
139 t5=cos(t3)*targ.fr1
140 t6=cos(t4)*targ.fr2
141 elseif (targ.fr_pattern .eq. 2) then !circular raster
142 t3=grnd()*2.*pi
143 t4=sqrt(grnd())*(targ.fr2-targ.fr1)+targ.fr1
144 t5=cos(t3)*t4
145 t6=sin(t3)*t4
146 else !no raster
147 t5=0.0
148 jones 1.1 t6=0.0
149 endif
150
151 main.target.x = main.target.x+t5
152 main.target.y = main.target.y+t6
153 main.target.z = (0.5-grnd())*targ.length+targ.zoffset
154 main.target.rastery = t6 !'raster' contribution to vert. pos.
155 main.target.rasterx = t5 !points left as look downstream. GAW
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 jones 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: fixed Em, generate electron angles
190 jones 1.1 ! doing_deuterium: fixed Em, generate electron fully and proton angles, calc Ep
191 ! doing_eep, A>2: generate electron and hadron energy and angles (calc Em/Pm).
192 ! doing_pion: generate electron energy/angles, hadron angles, p_fermi, Em.
193 ! doing_kaon: as doing_pion.
194 !
195 ! The above is summarized in the following table:
196 !
197 ! ELECTRON HADRON
198 ! ------------------ ------------------
199 ! E yptar xptar E yptar xptar p_fermi Em
200 !
201 !H(e,e'p) X X
202 !D(e,e'p) X X X X X
203 !A(e,e'p) X X X X X X
204 !----------------------------------------------------------------------
205 !H(e,e'pi) X X X X X
206 !A(e,e'pi) X X X X X X X
207 !H(e,e'K) X X X X X
208 !A(e,e'K) X X X X X X X
209 !----------------------------------------------------------------------
210 !phase_space X X X ? X X
211 jones 1.1 !
212 ! So our procedure is the following:
213 ! 1) Always generate electron yptar and xptar
214 ! 2) generate hadron yptar and xptar for all cases except H(e,e'p), D(e,e'p)
215 ! 3) Generate hadron E for A(e,e'p)
216 ! 4) generate p_fermi, Em for A(e,e'pi) and A(e,e'K) - calculate vertex.Mrec
217 ! (invariant mass M* of the A-1 recoil system)
218 ! 5) generate electron E for all but hydrogen elastic.
219 !
220 ! After we generate xptar/yptar/energy, we calculate physics angles (theta/phi),
221 ! momenta, unit vectors, etc... here and/or in complete_ev.
222 !
223 ! Note that there are also jacobians associated with some and/or all of
224 ! the above.
225 ! 1: We generate uniformly in xptar/yptar, not theta/phi. We define the
226 ! phase space volume (genvol contribution) as the product of the xptar/yptar
227 ! range, and have a jacobian for each event taking into account the mapping
228 ! between the solid angle on the unit sphere, and the dxptar/dyptar volume
229 ! (the jacobian is 1/cos**3(dtheta), where dtheta is the angle between the
230 ! event and the central spectrometer vector
231 ! 2: For the D(e,e'p), we take Em as fixed in order to calculate the proton
232 jones 1.1 ! energy. There is a jacobian ( |dEp'/dEm| ). It comes from integrating
233 ! over the energy conservation delta function: delta(E_D - E_p - E_n - Em).
234
235 ! Generate Electron Angles (all cases):
236 vertex.e.yptar=gen.e.yptar.min+grnd()*(gen.e.yptar.max-gen.e.yptar.min)
237 vertex.e.xptar=gen.e.xptar.min+grnd()*(gen.e.xptar.max-gen.e.xptar.min)
238
239 ! Generate Hadron Angles (all but H(e,e'p)):
240 if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon.or.doing_dvcs) then
241 vertex.p.yptar=gen.p.yptar.min+grnd()*(gen.p.yptar.max-gen.p.yptar.min)
242 vertex.p.xptar=gen.p.xptar.min+grnd()*(gen.p.xptar.max-gen.p.xptar.min)
243 endif
244
245 ! Generate Hadron Momentum (A(e,e'p) only).
246 if (doing_heavy) then
247 Emin = max(gen.p.E.min, gen.sumEgen.min - gen.e.E.max)
248 Emax = min(gen.p.E.max, gen.sumEgen.max - gen.e.E.min)
249 if (Emin.gt.Emax) goto 100
250 main.gen_weight=main.gen_weight*(Emax-Emin)/(gen.p.E.max-gen.p.E.min)
251 vertex.p.E = Emin + grnd()*(Emax-Emin)
252 vertex.p.P = sqrt(vertex.p.E**2 - Mh2)
253 jones 1.1 vertex.p.delta = 100.*(vertex.p.P-spec.p.P)/spec.p.P
254 endif
255
256 ! Generate Fermi Momentum and Em for A(e,e'pi) and A(e,e'K).
257 pfer=0.0
258 pferx=0.0
259 pfery=0.0
260 pferz=0.0
261 vertex.Em=0.0
262 if(doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon.or.doing_hedvcs)then
263 ranprob=grnd()
264 ii=1
265 do while (ranprob.gt.mprob(ii) .and. ii.lt.nump)
266 ii=ii+1
267 enddo
268 if (ii.eq.1) then
269 pferlo=0
270 else
271 pferlo=(pval(ii-1)+pval(ii))/2
272 endif
273 if (ii.eq.nump) then
274 jones 1.1 pferhi=pval(nump)
275 else
276 pferhi=(pval(ii)+pval(ii+1))/2
277 endif
278 pfer=pferlo+(pferhi-pferlo)*grnd()
279 ranth1=grnd()*2.-1.0
280 ranth=acos(ranth1)
281 ranph=grnd()*2.*pi
282 pferx=sin(ranth)*cos(ranph)
283 pfery=sin(ranth)*sin(ranph)
284 pferz=cos(ranth)
285
286 if(doing_deutpi.or.doing_deutkaon.or.doing_hedvcs)then !Em = binding energy
287 vertex.Em = Mp + Mn - targ.M
288 endif
289 if (doing_hepi .or. doing_hekaon.or.doing_hedvcs) then!Generate Em
290 call generate_em(pfer,vertex.Em) !Generate Em
291 endif
292 endif
293
294 ! Generate Electron Energy (all but hydrogen elastic)
295 jones 1.1 if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon.or.
296 1 doing_dvcs) then
297 Emin=gen.e.E.min
298 Emax=gen.e.E.max
299 if (doing_deuterium .or. doing_pion .or. doing_kaon .or.
300 1 doing_dvcs) then
301 Emin = max(Emin,gen.sumEgen.min)
302 Emax = min(Emax,gen.sumEgen.max)
303 else if (doing_heavy) then ! A(e,e'p)
304 Emin = max(Emin, gen.sumEgen.min - vertex.p.E)
305 Emax = min(Emax, gen.sumEgen.max - vertex.p.E)
306 endif
307
308 if (Emin.gt.Emax) goto 100
309 main.gen_weight=main.gen_weight*(Emax-Emin)/(gen.e.E.max-gen.e.E.min)
310 vertex.e.E = Emin + grnd()*(Emax-Emin)
311 vertex.e.P = vertex.e.E
312 vertex.e.delta = 100.*(vertex.e.P-spec.e.P)/spec.e.P
313 endif !not (doing_hyd_elast)
314
315
316 jones 1.1 ! Calculate the electron and proton PHYSICS angles from the spectrometer angles.
317
318 call physics_angles(spec.e.theta,spec.e.phi,
319 & vertex.e.xptar,vertex.e.yptar,vertex.e.theta,vertex.e.phi)
320 call physics_angles(spec.p.theta,spec.p.phi,
321 & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi)
322
323 ! Compute all non-generated quantities
324
325 if (debug(5)) write(6,*)'gen: calling comp_ev with false, main, vertex'
326 if (debug(3)) write(6,*)'gen: calling comp_ev with false, main, vertex'
327 if (debug(3)) write(6,*)'gen: Ein, E =',vertex.Ein,vertex.e.E
328 call complete_ev(main,vertex,success)
329
330 main.sigcc = 1.0
331
332 if (debug(2)) write(6,*)'gen: initial success =',success
333 if (.not.success) goto 100
334
335 ! ........ temporary storage of Trec for generated event
336
337 jones 1.1 main.Trec = vertex.Trec
338
339 ! Radiate the event, if requested. If we get an event weight of zero (i.e.
340 ! if we CAN'T radiate our event into the acceptance) then success is
341 ! false. generate_rad also set orig kinematics, and cuts on Em/Pm.
342 ! If not using_rad, must do these here.
343
344 if (using_rad) then
345 call generate_rad(main,vertex,orig,success)
346 if (debug(2)) write(6,*)'gen: after gen_rad, success =',success
347 else
348 success = .true.
349 if (doing_heavy) success =
350 > (vertex.Em .ge. VERTEXedge.Em.min .and.
351 > vertex.Em .le. VERTEXedge.Em.max .and.
352 > vertex.Pm .ge. VERTEXedge.Pm.min .and.
353 > vertex.Pm .le. VERTEXedge.Pm.max)
354 if (success) then
355 do i = 1, neventfields
356 orig.all(i) = vertex.all(i)
357 enddo
358 jones 1.1 endif
359 endif
360
361
362 100 if (debug(2)) write(6,*)'gen: final success =',success
363 if (debug(2)) write(6,*)'gen: ending'
364 return
365 end
366
367 !---------------------------------------------------------------------
368
369 subroutine complete_ev(main,vertex,success)
370
371 implicit none
372 include 'simulate.inc'
373
374 real*8 a, b, c, r, t, QA, QB, QC, radical
375 real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z
376 real*8 targ_new_x,targ_new_y
377 real*8 new_y_x,new_y_y,new_y_z,dummy
378 real*8 px,py,pz,qx,qy,qz
379 jones 1.1 real*8 targx,targy,targz,pol
380 real*8 cos_dtheta,y_spec,y_event
381 real*8 oop_x,oop_y
382 real*8 krel,krelx,krely,krelz
383 real*8 MM
384
385 logical success
386 record /event_main/ main
387 record /event/ vertex
388
389 !-----------------------------------------------------------------------
390 ! Calculate everything left in the /event/ structure, given all necessary
391 ! GENERATION values (some set of xptar,yptar,delta for both arms and p_fermi,
392 ! and p_fermi, depending on the scattering process: see table in generate.f
393 !
394 ! The SINGLE element of /event/ NOT computed here is sigcc.
395 !
396 ! Another small anomaly is that main.jacobian IS computed here (if not
397 ! doing_recon) ... this is because all the necessary terms have to be
398 ! computed here anyway to calculate /event/ qties.
399 !
400 jones 1.1 !-----------------------------------------------------------------------
401
402 ! Initialize
403
404 success = .false.
405 main.jacobian = 1.0
406 ! ... unit vector components of outgoing e,p
407 ! ... z is DOWNSTREAM, x is DOWN and y is LEFT looking downstream.
408
409 if (debug(4)) write(6,*)'comp_ev: at 1'
410 vertex.ue.x = sin(vertex.e.theta)*cos(vertex.e.phi)
411 vertex.ue.y = sin(vertex.e.theta)*sin(vertex.e.phi)
412 vertex.ue.z = cos(vertex.e.theta)
413 if (.not.doing_hyd_elast) then
414 vertex.up.x = sin(vertex.p.theta)*cos(vertex.p.phi)
415 vertex.up.y = sin(vertex.p.theta)*sin(vertex.p.phi)
416 vertex.up.z = cos(vertex.p.theta)
417 endif
418 if (debug(4)) write(6,*)'comp_ev: at 2'
419
420 ! First finish off the e side
421 jones 1.1 ! Calculate scattered electron energy for hydrogen/deuterium (e,e'p)
422
423 if (doing_hyd_elast) then
424 vertex.e.E = vertex.Ein*Mh/(Mh+vertex.Ein*(1.-vertex.ue.z))
425
426 if (vertex.e.E.gt.vertex.Ein) return
427 vertex.e.P = vertex.e.E
428 vertex.e.delta = (vertex.e.P - spec.e.P)*100./spec.e.P
429 if (debug(4)) write(6,*)'comp_ev: at 3'
430 endif
431
432 ! The q vector
433
434 if (debug(5)) write(6,*)'comp_ev: Ein,E,uez=',vertex.Ein,vertex.e.E,vertex.ue.z
435
436 vertex.omega = vertex.Ein - vertex.e.E
437 vertex.Q2 = 2*vertex.Ein*vertex.e.E*(1.-vertex.ue.z)
438 vertex.q = sqrt(vertex.Q2 + vertex.omega**2)
439
440 vertex.uq.x = - vertex.e.P*vertex.ue.x / vertex.q
441 vertex.uq.y = - vertex.e.P*vertex.ue.y / vertex.q
442 jones 1.1 vertex.uq.z =(vertex.Ein - vertex.e.P*vertex.ue.z)/ vertex.q
443 if (abs(vertex.uq.x**2+vertex.uq.y**2+vertex.uq.z**2-1).gt.0.01)
444 > stop 'Error in q vector normalization'
445 if (debug(4)) write(6,*)'comp_ev: at 5'
446
447 ! Determine vertex.Pm, vertex.Em, vertex.Mrec (needed to complete hadron).
448 ! For D(e,e'p) and A(e,e'p), vertex.Pm has to wait until after getting hadron.
449 ! For A(e,e'p) case, Em and Mrec have to wait until after calculating Pm.
450 ! For electroproduction, Em was generated at same time as pfermi (=0 for A=1).
451
452 if (doing_hyd_elast) then
453 vertex.Em = 0.0
454 vertex.Pm = 0.0
455 vertex.Mrec = 0.0
456 else if (doing_deuterium) then
457 vertex.Em = Mp + Mn - targ.M !2.2 MeV
458 vertex.Mrec = targ.M - targ.Mtar_struck + vertex.Em !neutron mass.
459 else if (doing_pion .or. doing_kaon .or. doing_dvcs) then
460 vertex.Pm = pfer
461 vertex.Mrec = targ.M - targ.Mtar_struck + vertex.Em
462 endif
463 jones 1.1
464 ! Now complete the p side.
465
466 if (doing_hyd_elast) then !p = q
467
468 vertex.up.x = vertex.uq.x
469 vertex.up.y = vertex.uq.y
470 vertex.up.z = vertex.uq.z
471 vertex.p.P = vertex.q
472 vertex.p.theta = acos(vertex.up.z)
473 if (sin(vertex.p.theta).eq.0) write(6,*) 'sin(vertex.p.theta)=',sin(vertex.p.theta)
474 if (abs(vertex.up.x/sin(vertex.p.theta)).gt.1)
475 > write(6,*) 'cos(phi)=',vertex.up.x/sin(vertex.p.theta)
476 vertex.p.phi = acos(min(1.d0,max(-1.d0,vertex.up.x/sin(vertex.p.theta))))
477
478 !...Get the cosine of the angle between the central spectrometer setting
479 !...and the real event (dot product of p and p0 unit vectors). For in-plane
480 !...spectrometers, xptar (=dx/dz) is just the x component of p (=dx)
481 !...over the component of p along the spectrometer direction (=dz=cos(dtheta)).
482 !...Knowing xptar, get yptar from cos(dtheta)=1/sqrt(1+xptar**2+yptar**2)
483 !...Choose sign of yptar by comparing y for the spectrometer to y of the
484 jones 1.1 !...event (projected to the xptar-yptar plane).
485
486 cos_dtheta = vertex.up.z * cos(spec.p.theta)
487 > + vertex.up.y * sin(spec.p.theta)*sin(spec.p.phi)
488 > + vertex.up.x * sin(spec.p.theta)*cos(spec.p.phi)
489 vertex.p.xptar = vertex.up.x / cos_dtheta
490 if ((1/cos_dtheta**2 - (1 + vertex.p.xptar**2)).lt.0) write (6,*)
491 > 'complete_ev: yptar**2=',(1/cos_dtheta**2 - (1 + vertex.p.xptar**2))
492 vertex.p.yptar = sqrt(max(0.d0,1/cos_dtheta**2-(1+vertex.p.xptar**2)))
493 y_spec = sin(spec.p.theta)*sin(spec.p.phi)
494 y_event = vertex.up.y/cos_dtheta !projected to plane perp. to spectrometer.
495 if (y_event .lt. y_spec) vertex.p.yptar = -vertex.p.yptar
496
497 vertex.p.E = sqrt(vertex.p.P**2+Mh2)
498 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
499 if (debug(4)) write(6,*)'comp_ev: at 6'
500
501 elseif (doing_deuterium) then !need Ep, and a jacobian.
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 jones 1.1 c = vertex.omega + 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 if (vertex.p.E.le.Mh) return
514 vertex.p.P = sqrt(vertex.p.E**2 - Mh2)
515 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
516
517 ! ........ the Jacobian here is |dEp'/dEm|
518 main.jacobian = (t*(c-vertex.p.E) + 2*c*vertex.p.E*(vertex.p.E-c)) /
519 1 (2*(a**2-c**2)*vertex.p.E + c*t)
520 main.jacobian = abs(main.jacobian)
521
522 elseif (doing_pion .or. doing_kaon .or. doing_dvcs) then
523 ! Hydrogen case. Coherent production (bound final state) is treated as
524 ! hydrogen, but with targ.Mtar_struck=targ.M, targ.Mtar_rec=bound final state.
525
526 jones 1.1 a = -1.*vertex.q*(vertex.uq.x*vertex.up.x+vertex.uq.y*vertex.up.y+vertex.uq.z*vertex.up.z)
527 b = vertex.q**2
528 c = vertex.omega + targ.M
529
530 ! For nuclei, correct for fermi motion and missing energy.
531
532 if (doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon
533 1 .or.doing_deutdvcs.or.doing_hedvcs) then
534 a = a - abs(pfer)*(pferx*vertex.up.x+pfery*vertex.up.y+pferz*vertex.up.z)
535 b = b + pfer**2 + 2*vertex.q*abs(pfer)*
536 > (pferx*vertex.uq.x+pfery*vertex.uq.y+pferz*vertex.uq.z)
537 c = c - sqrt(vertex.Mrec**2+pfer**2)
538 endif
539
540 t = c**2 - b + Mh2 - targ.Mrec_struck**2
541 QA = 4.*(a**2 - c**2)
542 QB = 4.*c*t
543 QC = -4.*a**2*Mh2 - t**2
544 radical = QB**2 - 4.*QA*QC
545 if (radical.lt.0) return
546 vertex.p.E = (-QB - sqrt(radical))/2./QA
547 jones 1.1 if (vertex.p.E.le.Mh) return
548 vertex.p.P = sqrt(vertex.p.E**2 - Mh2)
549 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
550
551 elseif (doing_phsp) then
552
553 vertex.p.P = spec.p.P !????? single arm phsp??
554 vertex.p.E= sqrt(Mh2+vertex.p.P**2)
555 vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P
556 if (debug(4)) write(6,*)'comp_ev: at 7.5',Mh2,vertex.p.E
557
558 endif
559
560 if (debug(4)) write(6,*)'comp_ev: at 7'
561
562 ! Compute some pion and kaon stuff
563
564 if (doing_pion .or. doing_kaon.or.doing_dvcs) then
565 main.epsilon=1./(1. + 2.*(1+vertex.omega**2/vertex.Q2)*tan(vertex.e.theta/2.)**2)
566 main.theta_pq=acos(vertex.up.x*vertex.uq.x+vertex.up.y*vertex.uq.y+vertex.up.z*vertex.uq.z)
567
568 jones 1.1 ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE.
569 ! Therefore, define a new system with the z axis parallel to q, and
570 ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q
571 ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane
572 ! particles closer to the downstream beamline than q.
573 ! phi=90 is above the horizontal plane when q points to the right, and
574 ! below the horizontal plane when q points to the left.
575 ! Also take into account the different definitions of x, y and z in
576 ! replay and SIMC:
577 ! As seen looking downstream: replay SIMC (old_simc)
578 ! x right down (left)
579 ! y down left (up)
580 ! z all have z pointing downstream
581 !
582 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
583
584 qx = -vertex.uq.y !convert to 'replay' coord. system
585 qy = vertex.uq.x
586 qz = vertex.uq.z
587 px = -vertex.up.y
588 py = vertex.up.x
589 jones 1.1 pz = vertex.up.z
590
591 dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
592 new_x_x = -qx*qz/dummy
593 new_x_y = -qy*qz/dummy
594 new_x_z = (qx**2 + qy**2)/dummy
595
596 dummy = sqrt(qx**2 + qy**2)
597 new_y_x = qy/dummy
598 new_y_y = -qx/dummy
599 new_y_z = 0.0
600
601 p_new_x = px*new_x_x + py*new_x_y + pz*new_x_z
602 p_new_y = px*new_y_x + py*new_y_y + pz*new_y_z
603
604 main.phi_pq = atan2(p_new_y,p_new_x) !atan2(y,x)=atan(y/x)
605 if(main.phi_pq .lt. 0.) main.phi_pq = main.phi_pq + 2.*pi
606 ! if ((p_new_x**2+p_new_y**2).eq.0.) then
607 ! main.phi_pq = 0.0
608 ! else
609 ! main.phi_pq = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2))
610 jones 1.1 ! endif
611 ! if (p_new_y.lt.0.) then
612 ! main.phi_pq = 2*pi - main.phi_pq
613 ! endif
614 !
615 if (debug(2)) then
616 write(6,*)'comp_ev: omega =',vertex.omega
617 write(6,*)'comp_ev: Q2 =',vertex.Q2
618 write(6,*)'comp_ev: theta_e =',vertex.e.theta
619 write(6,*)'comp_ev: epsilon =',main.epsilon
620 write(6,*)'comp_ev: theta_pq =',main.theta_pq
621 write(6,*)'comp_ev: phi_pq =',main.phi_pq
622 endif
623 endif !end of pion/kaon specific stuff.
624
625 ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND TARGET POLARIZATION.
626 ! Therefore, define a new system with the z axis parallel to q, and
627 ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q
628 ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane
629 ! particles closer to the downstream beamline than q.
630 ! phi=90 is above the horizontal plane when q points to the right, and
631 jones 1.1 ! below the horizontal plane when q points to the left.
632 ! Also take into account the different definitions of x, y and z in
633 ! replay and SIMC:
634 ! As seen looking downstream: replay SIMC (old_simc)
635 ! x right down (left)
636 ! y down left (up)
637 ! z all have z pointing downstream
638 !
639 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
640
641 qx = -vertex.uq.y !convert to 'replay' coord. system
642 qy = vertex.uq.x
643 qz = vertex.uq.z
644 targx = -targ_pol*sin(targ_Bangle) ! replay coordinates
645 targy = 0.0
646 targz = targ_pol*cos(targ_Bangle)
647
648 dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
649 new_x_x = -qx*qz/dummy
650 new_x_y = -qy*qz/dummy
651 new_x_z = (qx**2 + qy**2)/dummy
652 jones 1.1
653 dummy = sqrt(qx**2 + qy**2)
654 new_y_x = qy/dummy
655 new_y_y = -qx/dummy
656 new_y_z = 0.0
657
658 p_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z
659 p_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z
660
661 main.phi_targ = atan2(p_new_y,p_new_x) !atan2(y,x)=atan(y/x)
662 if(main.phi_targ.lt.0.) main.phi_targ = 2.*pi+main.phi_targ
663
664
665 ! CALCULATE ANGLE BETA BETWEEN REACTION PLANE AND TRANSVERSE TARGET
666 ! POLARIZATION.
667 ! Therefore, define a new system with the z axis parallel to q, and
668 ! with the y-direction defined by q x p_meson and x-axis given
669 ! by y x z
670 ! Also take into account the different definitions of x, y and z in
671 ! replay and SIMC:
672 ! As seen looking downstream: replay SIMC (old_simc)
673 jones 1.1 ! x right down (left)
674 ! y down left (up)
675 ! z all have z pointing downstream
676 !
677 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
678
679 qx = -vertex.uq.y !convert to 'replay' coord. system
680 qy = vertex.uq.x
681 qz = vertex.uq.z
682
683 pol = 1.0
684 targx = -targ_pol*sin(targ_Bangle) ! 'replay' coordinates
685 targy = 0.0
686 targz = targ_pol*cos(targ_Bangle)
687
688 px = -vertex.up.y
689 py = vertex.up.x
690 pz = vertex.up.z
691
692 dummy = sqrt((qy*pz-qz*py)**2 + (qz*px-qx*pz)**2 + (qx*py-qy*px)**2)
693 new_y_x = (qy*pz-qz*py)/dummy
694 jones 1.1 new_y_y = (qz*px-qx*pz)/dummy
695 new_y_z = (qx*py-qy*px)/dummy
696
697 dummy = sqrt((new_y_y*qz-new_y_z*qy)**2 + (new_y_z*qx-new_y_x*qz)**2
698 1 + (new_y_x*qy-new_y_y*qx)**2)
699
700 new_x_x = (new_y_y*qz - new_y_z*qy)/dummy
701 new_x_y = (new_y_z*qx - new_y_x*qz)/dummy
702 new_x_z = (new_y_x*qy - new_y_y*qx)/dummy
703
704
705 targ_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z
706 targ_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z
707
708
709 main.beta = atan2(targ_new_y,targ_new_x)
710 if(main.beta .lt. 0.) main.beta = 2*pi+main.beta
711
712 c if ((targ_new_x**2+targ_new_y**2).eq.0.) then
713 c main.beta = 0.0
714 c else
715 jones 1.1 c main.beta = acos(targ_new_x/sqrt(targ_new_x**2+targ_new_y**2))
716 c endif
717 c if (targ_new_y.lt.0.) then
718 c main.beta = 2*3.1415926536 - main.beta
719 c endif
720
721 dummy = sqrt((qx**2+qy**2+qz**2))*sqrt((targx**2+targy**2+targz**2))
722 main.theta_tarq = acos((qx*targx+qy*targy+qz*targz)/dummy)
723
724
725 if (debug(4)) write(6,*)'comp_ev: at 8'
726
727 ! Compute the Pm vector in in SIMC LAB system, with x down, and y to the left.
728 ! Computer Parallel, Perpendicular, and Out of Plane componenants.
729 ! Parallel is component along q_hat. Out/plane is component along
730 ! (e_hat) x (q_hat) (oop_x,oop_y are components of out/plane unit vector)
731 ! Perp. component is what's left: along (q_hat) x (oop_hat).
732 ! So looking along q, out of plane is down, perp. is left.
733
734 vertex.Pmx = vertex.p.P*vertex.up.x - vertex.q*vertex.uq.x
735 vertex.Pmy = vertex.p.P*vertex.up.y - vertex.q*vertex.uq.y
736 jones 1.1 vertex.Pmz = vertex.p.P*vertex.up.z - vertex.q*vertex.uq.z
737 vertex.Pmiss = sqrt(vertex.Pmx**2+vertex.Pmy**2+vertex.Pmz**2)
738 vertex.Emiss = vertex.omega + targ.M - vertex.p.E
739
740 !STILL NEED SIGN FOR PmPer!!!!!!
741
742 oop_x = -vertex.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat)
743 oop_y = vertex.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat)
744 vertex.PmPar = (vertex.Pmx*vertex.uq.x + vertex.Pmy*vertex.uq.y + vertex.Pmz*vertex.uq.z)
745 vertex.PmOop = (vertex.Pmx*oop_x + vertex.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2)
746 vertex.PmPer = sqrt( max(0.d0, vertex.Pm**2 - vertex.PmPar**2 - vertex.PmOop**2 ) )
747 c
748
749 c
750 if (debug(4)) write(6,*)'comp_ev: at 9',vertex.Pmx,vertex.Pmy,vertex.Pmz
751
752 ! Calculate Em, Pm, Mrec, Trec for all cases not already done.
753 ! For doing_heavy, get Mrec from nu+M=Ep+Erec, and Erec**2=Mrec**2+Pm**2
754 ! For (e,e'pi/K), could go back and determine momentum of recoil struck
755 ! particle, and get Mrec and Trec seperately for struck nucleon(hyperon)
756 ! and A-1 system. For now, just take Trec for the A-1 system, and ignore
757 jones 1.1 ! the recoiling struck nucleon (hyperon), so Trec=0 for hydrogen target.
758
759 if (doing_hyd_elast) then
760 vertex.Trec = 0.0
761 else if (doing_deuterium) then
762 vertex.Pm = vertex.Pmiss
763 vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec
764 else if (doing_heavy) then
765 vertex.Pm = vertex.Pmiss
766 vertex.Mrec = sqrt(vertex.Emiss**2-vertex.Pmiss**2)
767 vertex.Em = targ.Mtar_struck + vertex.Mrec - targ.M
768 vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec
769 else if (doing_hydpi .or. doing_hydkaon .or. doing_dvcs) then
770 vertex.Trec = 0.0
771 else if (doing_deutpi.or.doing_deutkaon.or.doing_hepi.or.doing_hekaon) then
772 vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec
773 endif
774
775 if (debug(5)) write(6,*) 'vertex.Pm,vertex.Trec,vertex.Em',vertex.Pm,vertex.Trec,vertex.Em
776 if (debug(4)) write(6,*)'comp_ev: at 10'
777
778 jones 1.1
779 ! calculate krel for deuteron/heavy pion(kaon). Deuteron is straightforward.
780 ! A>2 case is some approximation for 3He (DJG).
781
782 if (doing_deutpi .or. doing_deutkaon .or. doing_deutdvcs) then
783 if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) write(6,*) 'BAD MM!!!!! Emiss,Pmiss=',vertex.Emiss, vertex.Pmiss
784 MM = sqrt(max(0.d0,vertex.Emiss**2-vertex.Pmiss**2))
785 krel = sqrt( max(0.d0,MM**2-4.*targ.Mrec_struck**2) )
786 else if (doing_hepi .or. doing_hekaon .or. doing_hedvcs) then
787 if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) write(6,*) 'BAD MM!!!!! Emiss,Pmiss=',vertex.Emiss, vertex.Pmiss
788 MM = sqrt(max(0.d0,vertex.Emiss**2-vertex.Pmiss**2))
789 krelx = vertex.Pmx + 1.5*pferx*pfer
790 krely = vertex.Pmy + 1.5*pfery*pfer
791 krelz = vertex.Pmz + 1.5*pferz*pfer
792 krel = sqrt(krelx**2+krely**2+krelz**2)
793 if(vertex.Em.lt.6.0) krel = -krel !bound state test???
794 endif
795 ntup.krel = krel
796
797 ! Determine PHYSICS scattering angles theta/phi for the two spectrometer
798 ! vectors, and the Jacobian which we must include in our xsec computation
799 jones 1.1 ! to take into account the fact that we're not generating in the physics
800 ! solid angle d(omega) but in 'spectrometer angles' th, ph.
801
802 call physics_angles(spec.e.theta,spec.e.phi,
803 & vertex.e.xptar,vertex.e.yptar,vertex.e.theta,vertex.e.phi)
804 call physics_angles(spec.p.theta,spec.p.phi,
805 & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi)
806
807 r = sqrt(1.+vertex.e.yptar**2+vertex.e.xptar**2) !r=cos(theta-theta0)
808 main.jacobian = main.jacobian / r**3 !1/cos**3(theta-theta0)
809
810 if (doing_heavy .or. doing_pion .or. doing_kaon .or. doing_dvcs) then
811 r = sqrt(1.+vertex.p.yptar**2+vertex.p.xptar**2)
812 main.jacobian = main.jacobian / r**3 !1/cos**3(theta-theta0)
813 endif
814
815 if (debug(4)) write(6,*)'comp_ev: at 11'
816
817 ! The effective target thickness that the scattered particles see, and the
818 ! resulting energy losses
819
820 jones 1.1 call trip_thru_target (2, main.target.z-targ.zoffset, vertex.e.E,
821 > vertex.e.theta, main.target.Eloss(2), main.target.teff(2),Me,1)
822 call trip_thru_target (3, main.target.z-targ.zoffset, vertex.p.E,
823 > vertex.p.theta, main.target.Eloss(3), main.target.teff(3),Mh,1)
824 if (.not.using_Eloss) then
825 main.target.Eloss(2) = 0.0
826 main.target.Eloss(3) = 0.0
827 endif
828 if (debug(4)) write(6,*)'comp_ev: at 12'
829
830 ! Initialize parameters necessary for radiative corrections
831
832 call radc_init_ev(main,vertex)
833
834 ! Success code
835
836 success = .true.
837 if(debug(2)) write(6,*)'comp_ev: at end...'
838 return
839 end
840
841 jones 1.1 !---------------------------------------------------------------------
842
843 subroutine complete_recon_ev(recon,success)
844
845 implicit none
846 include 'simulate.inc'
847
848 real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z
849 real*8 targ_new_x,targ_new_y
850 real*8 new_y_x,new_y_y,new_y_z,dummy
851 real*8 px,py,pz,qx,qy,qz
852 real*8 targx,targy,targz,pol
853 real*8 W2
854 real*8 oop_x,oop_y
855 real*8 mm,mmA,mm2,mmA2,t
856
857 logical success
858 record /event/ recon
859
860 !-----------------------------------------------------------------------
861 ! Calculate everything left in the /event/ structure, given
862 jones 1.1 ! recon.Ein, and the following for both recon.e AND recon.p:
863 ! delta,xptar,yptar, z,P,E,theta,phi (with E,P corrected for eloss, if desired).
864 !
865 ! The SINGLE element of /event/ NOT computed here is sigcc (for (e,e'p)).
866 !
867 !-----------------------------------------------------------------------
868
869 ! Initialize
870
871 success = .false.
872
873 recon.Ein = Ebeam_vertex_ave !lowered by most probable eloss (init.f)
874
875 ! ... unit vector components of outgoing e,p
876 ! ... z is DOWNSTREAM, x is DOWN and y is LEFT looking downstream.
877
878 if (debug(4)) write(6,*)'comp_rec_ev: at 1'
879 recon.ue.x = sin(recon.e.theta)*cos(recon.e.phi)
880 recon.ue.y = sin(recon.e.theta)*sin(recon.e.phi)
881 recon.ue.z = cos(recon.e.theta)
882 recon.up.x = sin(recon.p.theta)*cos(recon.p.phi)
883 jones 1.1 recon.up.y = sin(recon.p.theta)*sin(recon.p.phi)
884 recon.up.z = cos(recon.p.theta)
885 if (debug(4)) write(6,*)'comp_rec_ev: at 2'
886
887 ! The q vector
888
889 if (debug(5)) write(6,*)'comp_rec_ev: Ein,E,uez=',recon.Ein,recon.e.E,recon.ue.z
890 recon.omega = recon.Ein - recon.e.E
891 recon.Q2 = 2*recon.Ein*recon.e.E*(1-recon.ue.z)
892 recon.q = sqrt(recon.Q2 + recon.omega**2)
893 recon.uq.x = - recon.e.P*recon.ue.x / recon.q
894 recon.uq.y = - recon.e.P*recon.ue.y / recon.q
895 recon.uq.z =(recon.Ein - recon.e.P*recon.ue.z)/ recon.q
896 W2 = targ.M**2 + 2.*targ.M*recon.omega - recon.Q2
897 recon.W = sqrt(abs(W2)) * W2/abs(W2)
898 if (debug(4)) write(6,*)'comp_rec_ev: at 5'
899
900 ! Now complete the p side
901
902 if(doing_phsp)then
903 recon.p.P=spec.p.P
904 jones 1.1 recon.p.E=sqrt(Mh2+recon.p.P**2)
905 if (debug(4)) write(6,*)'comp_rec_ev: at 7.5',Mh2,recon.p.E
906 endif
907
908 ! Compute some pion/kaon stuff.
909
910 recon.epsilon=1./(1. + 2.*(1+recon.omega**2/recon.Q2)*tan(recon.e.theta/2.)**2)
911 recon.theta_pq=acos(recon.up.x*recon.uq.x+recon.up.y*recon.uq.y+recon.up.z*recon.uq.z)
912
913 ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE.
914 ! Therefore, define a new system with the z axis parallel to q, and
915 ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q
916 ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane
917 ! particles closer to the downstream beamline than q.
918 ! phi=90 is above the horizontal plane when q points to the right, and
919 ! below the horizontal plane when q points to the left.
920 ! Also take into account the different definitions of x, y and z in
921 ! replay and SIMC:
922 ! As seen looking downstream: replay SIMC (old_simc)
923 ! x right down (left)
924 ! y down left (up)
925 jones 1.1 ! z all have z pointing downstream
926 !
927 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
928
929 qx = -recon.uq.y !convert to 'replay' coord. system
930 qy = recon.uq.x
931 qz = recon.uq.z
932 px = -recon.up.y
933 py = recon.up.x
934 pz = recon.up.z
935
936 dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
937 new_x_x = -qx*qz/dummy
938 new_x_y = -qy*qz/dummy
939 new_x_z = (qx**2 + qy**2)/dummy
940
941 dummy = sqrt(qx**2 + qy**2)
942 new_y_x = qy/dummy
943 new_y_y = -qx/dummy
944 new_y_z = 0.0
945
946 jones 1.1 p_new_x = px*new_x_x + py*new_x_y + pz*new_x_z
947 p_new_y = px*new_y_x + py*new_y_y + pz*new_y_z
948
949 recon.phi_pq = atan2(p_new_y,p_new_x)
950 if(recon.phi_pq.lt.0.) recon.phi_pq = 2.*pi + recon.phi_pq
951 c if ((p_new_x**2+p_new_y**2).eq.0.) then
952 c recon.phi_pq = 0.0
953 c else
954 c recon.phi_pq = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2))
955 c endif
956 c if (p_new_y.lt.0.) then
957 c recon.phi_pq = 2*3.1415926536 - recon.phi_pq
958 c endif
959
960 ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND TARGET POL.
961 ! Therefore, define a new system with the z axis parallel to q, and
962 ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q
963 ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane
964 ! particles closer to the downstream beamline than q.
965 ! phi=90 is above the horizontal plane when q points to the right, and
966 ! below the horizontal plane when q points to the left.
967 jones 1.1 ! Also take into account the different definitions of x, y and z in
968 ! replay and SIMC:
969 ! As seen looking downstream: replay SIMC (old_simc)
970 ! x right down (left)
971 ! y down left (up)
972 ! z all have z pointing downstream
973 !
974 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
975
976 qx = -recon.uq.y !convert to 'replay' coord. system
977 qy = recon.uq.x
978 qz = recon.uq.z
979 pol = 1.0
980 targx = -targ_pol*sin(targ_Bangle) ! 'replay' coordinates
981 targy = 0.0
982 targz = targ_pol*cos(targ_Bangle)
983
984 dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
985 new_x_x = -qx*qz/dummy
986 new_x_y = -qy*qz/dummy
987 new_x_z = (qx**2 + qy**2)/dummy
988 jones 1.1
989 dummy = sqrt(qx**2 + qy**2)
990 new_y_x = qy/dummy
991 new_y_y = -qx/dummy
992 new_y_z = 0.0
993
994 p_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z
995 p_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z
996
997
998 recon.phi_targ = atan2(p_new_y,p_new_x)
999 if(recon.phi_targ.lt.0.) recon.phi_targ = 2.*pi + recon.phi_targ
1000 c if ((p_new_x**2+p_new_y**2).eq.0.) then
1001 c recon.phi_targ = 0.0
1002 c else
1003 c recon.phi_targ = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2))
1004 c endif
1005 c if (p_new_y.lt.0.) then
1006 c recon.phi_targ = 2*3.1415926536 - recon.phi_targ
1007 c endif
1008
1009 jones 1.1 c recon.phi_targ = atan2(p_new_y,p_new_x)
1010 ! CALCULATE ANGLE BETA BETWEEN REACTION PLANE AND TRANSVERSE TARGET
1011 ! POLARIZATION.
1012 ! Therefore, define a new system with the z axis parallel to q, and
1013 ! with the y-direction defined by q x p_meson and x-axis given
1014 ! by y x z
1015 ! Also take into account the different definitions of x, y and z in
1016 ! replay and SIMC:
1017 ! As seen looking downstream: replay SIMC (old_simc)
1018 ! x right down (left)
1019 ! y down left (up)
1020 ! z all have z pointing downstream
1021 !
1022 ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc
1023
1024 qx = -recon.uq.y !convert to 'replay' coord. system
1025 qy = recon.uq.x
1026 qz = recon.uq.z
1027
1028 pol = 1.0
1029 targx = -targ_pol*sin(targ_Bangle) ! 'replay' coordinates
1030 jones 1.1 targy = 0.0
1031 targz = targ_pol*cos(targ_Bangle)
1032 px = -recon.up.y
1033 py = recon.up.x
1034 pz = recon.up.z
1035
1036 dummy = sqrt((qy*pz-qz*py)**2 + (qz*px-qx*pz)**2 + (qx*py-qy*px)**2)
1037 new_y_x = (qy*pz-qz*py)/dummy
1038 new_y_y = (qz*px-qx*pz)/dummy
1039 new_y_z = (qx*py-qy*px)/dummy
1040
1041 dummy = sqrt((new_y_y*qz-new_y_z*qy)**2 + (new_y_z*qx-new_y_x*qz)**2
1042 1 + (new_y_x*qy-new_y_y*qx)**2)
1043
1044 new_x_x = (new_y_y*qz - new_y_z*qy)/dummy
1045 new_x_y = (new_y_z*qx - new_y_x*qz)/dummy
1046 new_x_z = (new_y_x*qy - new_y_y*qx)/dummy
1047
1048
1049 targ_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z
1050 targ_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z
1051 jones 1.1
1052 recon.beta = atan2(targ_new_y,targ_new_x)
1053 if(recon.beta.lt.0.) recon.beta = 2.*pi + recon.beta
1054
1055 c if ((targ_new_x**2+targ_new_y**2).eq.0.) then
1056 c recon.beta = 0.0
1057 c else
1058 c recon.beta = acos(targ_new_x/sqrt(targ_new_x**2+targ_new_y**2))
1059 c endif
1060 c if (targ_new_y.lt.0.) then
1061 c recon.beta = 2*3.1415926536 - recon.beta
1062 c endif
1063
1064
1065 dummy = sqrt((qx**2+qy**2+qz**2))*sqrt((targx**2+targy**2+targz**2))
1066 recon.theta_tarq = acos((qx*targx+qy*targy+qz*targz)/dummy)
1067 c write(6,*) 'cheesy poofs', recon.beta, p_new_x, p_new_y,qx,targy
1068 c if (debug(2)) then
1069 c write(6,*)'comp_rec_ev: omega =',recon.omega
1070 c write(6,*)'comp_rec_ev: Q2 =',recon.Q2
1071 c write(6,*)'comp_rec_ev: theta_e =',recon.e.theta
1072 jones 1.1 c write(6,*)'comp_rec_ev: epsilon =',recon.epsilon
1073 c write(6,*)'comp_rec_ev: theta_pq =',recon.theta_pq
1074 c write(6,*)'comp_rec_ev: phi_pq =',recon.phi_pq
1075 c endif
1076
1077 if (debug(4)) write(6,*)'comp_rec_ev: at 8'
1078
1079 ! Compute the Pm vector in in SIMC LAB system, with x down, and y to the left.
1080 ! Computer Parallel, Perpendicular, and Out of Plane componenants.
1081 ! Parallel is component along q_hat. Out/plane is component along
1082 ! (e_hat) x (q_hat) (oop_x,oop_y are components of out/plane unit vector)
1083 ! Perp. component is what's left: along (q_hat) x (oop_hat).
1084 ! So looking along q, out of plane is down, perp. is left.
1085
1086 recon.Pmx = recon.p.P*recon.up.x - recon.q*recon.uq.x
1087 recon.Pmy = recon.p.P*recon.up.y - recon.q*recon.uq.y
1088 recon.Pmz = recon.p.P*recon.up.z - recon.q*recon.uq.z
1089 recon.Pm = sqrt(recon.Pmx**2+recon.Pmy**2+recon.Pmz**2)
1090
1091 !STILL NEED SIGN FOR PmPer!!!!!!
1092
1093 jones 1.1 oop_x = -recon.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat)
1094 oop_y = recon.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat)
1095 recon.PmPar = (recon.Pmx*recon.uq.x + recon.Pmy*recon.uq.y + recon.Pmz*recon.uq.z)
1096 recon.PmOop = (recon.Pmx*oop_x + recon.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2)
1097 recon.PmPer = sqrt( max(0.d0, recon.Pm**2 - recon.PmPar**2 - recon.PmOop**2 ) )
1098
1099 if (doing_pion .or. doing_kaon .or. doing_dvcs) then
1100 recon.Em = recon.omega + targ.Mtar_struck - recon.p.E
1101 mm2 = recon.Em**2 - recon.Pm**2
1102 mm = sqrt(abs(mm2)) * abs(mm2)/mm2
1103 mmA2= (recon.omega + targ.M - recon.p.E)**2 - recon.Pm**2
1104 mmA = sqrt(abs(mmA2)) * abs(mmA2)/mmA2
1105 t = recon.Q2 - Mh2
1106 > + 2*(recon.omega*recon.p.E - recon.p.P*recon.q*cos(recon.theta_pq))
1107 ntup.mm = mm
1108 ntup.mmA = mmA
1109 ntup.t = t
1110 endif
1111
1112 ! Calculate Trec, Em. Trec for (A-1) system (eep), or for struck nucleon (pi/K)
1113 ! Note that there are other ways to calculate 'Em' for the pion/kaon case.
1114 jones 1.1
1115 if (doing_hyd_elast) then
1116 recon.Trec = 0
1117 recon.Em = recon.Ein + Mp - recon.e.E - recon.p.E - recon.Trec
1118 else if (doing_deuterium .or. doing_heavy) then
1119 recon.Trec = sqrt(recon.Pm**2+targ.Mrec**2) - targ.Mrec
1120 recon.Em = recon.Ein + Mp - recon.e.E - recon.p.E - recon.Trec
1121 else if (doing_pion) then
1122 recon.Em = recon.omega + targ.Mtar_struck - recon.p.E
1123 else if (doing_kaon) then
1124 recon.Em = recon.omega + targ.Mtar_struck - recon.p.E
1125 else if (doing_dvcs) then
1126 recon.Em = recon.omega + targ.Mtar_struck - recon.p.E
1127 endif
1128
1129 if (debug(5)) write(6,*) 'recon.Pm,recon.Trec,recon.Em',recon.Pm,recon.Trec,recon.Em
1130 if (debug(4)) write(6,*)'comp_rec_ev: at 10'
1131
1132 success=.true.
1133 return
1134 end
1135 jones 1.1
1136 !------------------------------------------------------------------------
1137
1138 subroutine complete_main(force_sigcc,main,vertex,recon,success)
1139
1140 implicit none
1141 include 'simulate.inc'
1142
1143 integer i, iPm1
1144 real*8 a, b, r, frac, peepi, peeK
1145 real*8 survivalprob
1146 real*8 weight, width, sigep, deForest, tgtweight
1147 logical force_sigcc, success
1148 record /event_main/ main
1149 record /event/ vertex, recon
1150
1151 !-----------------------------------------------------------------------
1152 ! Calculate everything left in the /main/ structure that hasn't been
1153 ! required up til now. This routine is called ONLY after we
1154 ! know an event IS going to contribute, don't want to be doing these
1155 ! calculations if they're not needed!
1156 jones 1.1 !
1157 ! A small anomaly is that sigcc from the /event/ structure is computed
1158 ! here since it's not needed til now, and was not calculated by
1159 ! COMPLETE_ev.
1160 !-----------------------------------------------------------------------
1161
1162
1163 ! Initialize success code
1164
1165 if (debug(2)) write(6,*)'comp_main: entering...'
1166 success = .false.
1167
1168 ! The spectral function weighting
1169
1170 if (doing_hyd_elast.or.doing_pion.or.doing_kaon.or.doing_dvcs.or.doing_phsp) then !no SF.
1171 main.SF_weight=1.0
1172 else if (doing_deuterium .or. doing_heavy) then
1173 main.SF_weight = 0.0
1174 do i=1,nrhoPm
1175 weight = 0.0
1176
1177 jones 1.1 ! ... use linear interpolation to determine rho(pm)
1178
1179 r = (vertex.Pm-Pm_theory(i).min)/Pm_theory(i).bin
1180 if (r.ge.0 .and. r.le.Pm_theory(i).n) then
1181 iPm1 = nint(r)
1182 if (iPm1.eq.0) iPm1 = 1
1183 if (iPm1.eq.Pm_theory(i).n) iPm1 = Pm_theory(i).n-1
1184 frac = r+0.5-float(iPm1)
1185 b = theory(i,iPm1)
1186 a = theory(i,iPm1+1)-b
1187 weight = a*frac + b
1188 endif
1189
1190 if (doing_heavy) then
1191
1192 width=Emsig_theory(i)/2.0
1193 if (vertex.Em.lt.E_Fermi) weight = 0.0
1194 weight=weight/pi/Em_int_theory(i) * width/
1195 > ((vertex.Em-Em_theory(i))**2+width**2)
1196 endif
1197 main.SF_weight = main.SF_weight+weight*nprot_theory(i)
1198 jones 1.1 enddo ! <i=1,nrhoPm>
1199 endif
1200
1201 ! ... if we have come up with weight<=, quit now (avoid weight<0 in ntuple)
1202 ! ... unless FORCE_SIGCC is on (to get central sigcc).
1203
1204 if (debug(5))write(6,*)'comp_main: calculating cross section'
1205 if (main.SF_weight.le.0 .and. .not.force_sigcc) return
1206
1207 ! ... Cross section, for BOTH vertex and recon (where possible)
1208 ! ... Note that all of these are cross section per nucleon. For A(e,e'p)
1209 ! ... this becomes cross section per nucleus because of the weighting of
1210 ! ... the spectral function. For pion/kaon production, we explicitly
1211 ! ... apply a weight for the number of nucleons involved (protons for pi+ or Y0,
1212 ! ... neutrons for pi- or Y-) (Y=hyperon)
1213
1214 tgtweight = 1.0
1215
1216 if (doing_phsp) then
1217 main.sigcc = 1.0
1218 main.sigcc_recon = 1.0
1219 jones 1.1
1220 elseif (doing_hyd_elast) then
1221 main.sigcc = sigep(vertex)
1222 main.sigcc_recon = sigep(recon)
1223
1224 elseif (doing_deuterium.or.doing_heavy) then
1225 main.sigcc = deForest(vertex)
1226 main.sigcc_recon = deForest(recon)
1227
1228 elseif (doing_pion) then
1229 main.sigcc = peepi(vertex,main)
1230 main.sigcc_recon = 1.0
1231 if (which_pion.eq.1 .or. which_pion.eq.11) then !OK for coherent???
1232 tgtweight = targ.N
1233 else
1234 tgtweight = targ.Z
1235 endif
1236
1237 elseif (doing_kaon) then
1238 main.sigcc = peeK(vertex,main,survivalprob)
1239 main.sigcc_recon = 1.0
1240 jones 1.1 if (which_kaon.eq.2 .or. which_kaon.eq.12) then !OK for coherent???
1241 tgtweight = targ.N
1242 else
1243 tgtweight = targ.Z
1244 endif
1245 elseif (doing_dvcs) then
1246 main.sigcc = 1.0
1247 main.sigcc_recon = 1.0
1248 if (which_dvcs.eq.0) then !OK for coherent???
1249 tgtweight = targ.Z
1250 else if(which_dvcs.eq.10) then
1251 tgtweight = targ.Z
1252 else
1253 tgtweight = targ.N
1254 endif
1255 else
1256 main.sigcc = 1.0
1257 main.sigcc_recon = 1.0
1258 endif
1259
1260 if (debug(3)) then
1261 jones 1.1 write(6,*)'======================================'
1262 write(6,*)'complete_main:'
1263 write(6,*)'theta_e =',vertex.e.theta
1264 write(6,*)'q mag =',vertex.q
1265 write(6,*)'omega =',vertex.omega
1266 write(6,*)'Q2 =',vertex.Q2/1000000.
1267 write(6,*)'Ein =',vertex.Ein
1268 write(6,*)'hadron mom =',vertex.p.P
1269 write(6,*)'e mom =',vertex.e.P
1270 write(6,*)'mass =',Mh
1271 write(6,*)'epsilon =',main.epsilon
1272 write(6,*)'phi_pq =',main.phi_pq
1273 write(6,*)'theta_pq =',main.theta_pq
1274 write(6,*)'======================================'
1275 endif
1276
1277 ! The total contributing weight from this event -- this weight is
1278 ! proportional to # experimental counts represented by the event.
1279 ! Apply survival probability to kaons if we're not modeling decay.
1280
1281 main.weight = main.SF_weight*main.jacobian*main.gen_weight*main.sigcc
1282 jones 1.1 main.weight = main.weight * tgtweight !correct for #/nucleons involved
1283 if (doing_kaon .and. .not.doing_decay)
1284 > main.weight = main.weight*survivalprob
1285 if (debug(5))write(6,*) 'gen_weight = ',main.gen_weight,
1286 > main.jacobian,main.sigcc
1287
1288 success = .true.
1289
1290 if (debug(2)) write(6,*)'comp_main: ending, success =',success
1291 return
1292 end
1293
1294
1295 subroutine physics_angles(theta0,phi0,dx,dy,theta,phi)
1296
1297 !Generate physics angles in lab frame. Theta is angle from beamline.
1298 !phi is projection on x-y plane (so phi=0 is down, sos=pi/2, hms=3*pi/2.
1299 !
1300 !theta=acos( (cos(theta0)-dy*sin(theta0)*sin(phi0))/ sqrt(1+dx**2+dy**2) )
1301 !phi=atan( (dy*cos(theta0)+sin(theta0)*sin(phi0)) / (sin(theta0)*cos(phi0)+dx) )
1302 !
1303 jones 1.1 ! Note that these formulae assume phi0=pi/2 or 3*pi/2 (thus the sin(phi0)
1304 ! gives a -/+ sign for the HMS/SOS). Thus, we set the cos(phi0) term to zero.
1305
1306 real*8 dx,dy !dx/dy (xptar/yptar) for event.
1307 real*8 theta0,phi0 !central physics angles of spectrometer.
1308 real*8 theta,phi !physics angles for event.
1309 real*8 r,sinth,costh,sinph,cosph !intermediate variables.
1310 real*8 tmp
1311
1312 include 'simulate.inc'
1313
1314 costh = cos(theta0)
1315 sinth = sin(theta0)
1316 sinph = sin(phi0)
1317 cosph = cos(phi0)
1318 r = sqrt( 1. + dx**2 + dy**2 )
1319
1320 if (abs(cosph).gt.0.0001) then !phi not at +/- pi/2
1321 write(6,*) 'theta,phi will be incorrect if phi0 <> pi/2 or 3*pi/2'
1322 write(6,*) 'phi0=',phi0,'=',phi0*180/pi,'degrees'
1323 endif
1324 jones 1.1
1325 tmp=(costh - dy*sinth*sinph) / r
1326 if (abs(tmp).gt.1) write(6,*) 'tmp=',tmp
1327 theta = acos( (costh - dy*sinth*sinph) / r )
1328 if (dx.ne.0.0) then
1329 phi = atan( (dy*costh + sinth*sinph) / dx ) !gives -90 to 90 deg.
1330 if (phi.le.0) phi=phi+pi !make 0 to 180 deg.
1331 if (sinph.lt.0.) phi=phi+pi !add pi to phi for HMS
1332 else
1333 phi = phi0
1334 endif
1335
1336 return
1337 end
|
Return to
event_orig.f
CVS log
Up to [HallC] / Poltar