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