Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dpm25qelpo.f
Go to the documentation of this file.
1  SUBROUTINE qel_pol(ENU, LTYP,P21,P22,P23,P24,P25)
2  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3  common/req/k2
4 C COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
5  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
6  common/pol/polarx(4),pmodul
7  parameter(nmxhkk= 89998)
8 c PARAMETER (NMXHKK=25000)
9  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
10  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
11  +(4,nmxhkk)
12 C
13 C REAL*4 P,V
14  COMMON /cbad/ lbad, nbad
15  COMMON /cnuc/ xmn,xmn2,pfermi,efermi,ebind,eb2,c0
16  COMMON /cconlun/ luna
17  DATA inipri/0/
18 C WRITE(*,*)'Input Enu, Nu-typ, N-ev '
19 C READ(*,*) ENU,LTYP,NUM
20  CALL mass_ini
21 C num=1000
22 C ltyp=5 !select neutrino tau (anutau=6, numu=3, anumu = 4)
23 C DO j=5,30 !loop for different neutrino energies
24 C enu=j
25  enup=enu
26 C DO K2=1,NUM
27 C CALL GEN_QEL(ENU,LTYP)
28  CALL gen_qel(enu, ltyp,p21,p22,p23,p24,p25)
29  IF(inipri.LT.10)THEN
30  inipri=inipri+1
31  WRITE(6,23) enup,(polarx(j1),j1=1,3)
32  WRITE(6,*) enup,p(4,4),p(5,4)
33  WRITE(6,*)p21,p22,p23,p24,p25
34  ENDIF
35 C WHKK(1,4)=POLARX(1)
36 C WHKK(2,4)=POLARX(2)
37 C WHKK(3,4)=POLARX(3)
38 C WHKK(4,4)=POLARX(4)
39  23 FORMAT(4(1x,f10.6))
40 C END DO
41 C END DO
42  RETURN
43  END
44 
45 
46 C==================================================================
47 C. Generation of a Quasi-Elastic neutrino scattering
48 C==================================================================
49 
50 C SUBROUTINE GEN_QEL (ENU, LTYP)
51  SUBROUTINE gen_qel (ENU, LTYP,P21,P22,P23,P24,P25)
52 C...Generate a quasi-elastic neutrino/antineutrino
53 C. Interaction on a nuclear target
54 C. INPUT : LTYP = neutrino type (1,...,6)
55 C. ENU (GeV) = neutrino energy
56 C----------------------------------------------------
57  IMPLICIT DOUBLE PRECISION (a-h,o-z)
58  common/req/k2
59 C COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
60  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
61  dimension rot(3,3),pi(3),po(3)
62 C REAL*4 P,V
63  COMMON /cbad/ lbad, nbad
64  COMMON /cnuc/ xmn,xmn2,pfermi,efermi,ebind,eb2,c0
65  COMMON /cconlun/ luna
66 CJR+
67  COMMON /nucimp/ prmom(5,248),tamom(5,248), prmfep,prmfen,tamfep,
68  +tamfen, prefep,prefen,taefep,taefen, prepot(210),taepot(210),
69  +prebin,taebin,fermod,etacou
70  COMMON /tautau/q(4,5),
71  + etl,pxl,pyl,pzl,
72  + etb,pxb,pyb,pzb,etn,pxn,pyn,pzn
73  COMMON /neutyy/neutyp,neudec
74  COMMON /nucros/dsigsu,dsigmc,ndsig
75  COMMON /qhelp/eeenu,kktyp
76  DATA ininu/0/
77 CJR-
78 C REAL*8 DBETA(3)
79 C REAL*8 MN(2), ML0(6), ML, ML2, MI, MI2, MF, MF2
80  dimension dbeta(3),dbetb(3),amn(2),aml0(6)
81  DATA amn /0.93827231d0, 0.93956563d0/
82  DATA aml0 /2*0.51100d-03,2*0.105659d0, 2*1.777d0/
83 C DATA PFERMI/0.22D0/
84  DATA inipri/0/
85 CGB+...Binding Energy
86  DATA ebind/0.008d0/
87 CGB-...
88  ininu=ininu+1
89  IF(ininu.EQ.1)ndsig=0
90  lbad = 0
91  enu0=enu
92 c write(*,*) enu0
93 C...Lepton mass
94  aml = aml0(ltyp) ! massa leptoni
95  aml2 = aml**2 ! massa leptoni **2
96 C...Particle labels (LUND)
97  n = 5
98  k(1,1) = 21
99  k(2,1) = 21
100  k(3,1) = 21
101  k(3,3) = 1
102  k(4,1) = 1
103  k(4,3) = 1
104  k(5,1) = 1
105  k(5,3) = 2
106  k0 = (ltyp-1)/2 ! 2
107  k1 = ltyp/2 ! 2
108  ka = 12 + 2*k0 ! 16
109  is = -1 + 2*ltyp - 4*k1 ! -1 +10 -8 = 1
110  k(1,2) = is*ka
111  k(4,2) = is*(ka-1)
112  k(3,2) = is*24
113  lnu = 2 - ltyp + 2*k1 ! 2 - 5 + 2 = - 1
114  IF (lnu .EQ. 2) THEN
115  k(2,2) = 2212
116  k(5,2) = 2112
117  ami = amn(1)
118  amf = amn(2)
119 CJR+
120  pfermi=tamfen
121 CJR-
122  ELSE
123  k(2,2) = 2112
124  k(5,2) = 2212
125  ami = amn(2)
126  amf = amn(1)
127 CJR+
128  pfermi=tamfep
129 CJR-
130  ENDIF
131  ami2 = ami**2
132  amf2 = amf**2
133 
134  DO igb=1,5
135  p(3,igb) = 0.
136  p(4,igb) = 0.
137  p(5,igb) = 0.
138  END DO
139 
140  ntry = 0
141 CGB+...
142  efmax = sqrt(pfermi**2 + ami2) -ami ! max. Fermi Energy
143  enwell = efmax + ebind ! depth of nuclear potential well
144 CGB-...
145 
146  100 CONTINUE
147 
148 C...4-momentum initial lepton
149  p(1,5) = 0. ! massa
150  p(1,4) = enu0 ! energia
151  p(1,1) = 0. ! px
152  p(1,2) = 0. ! py
153  p(1,3) = enu0 ! pz
154 
155 C PF = PFERMI*PYR(0)**(1./3.)
156 c write(23,*) PYR(0)
157 c write(*,*) 'Pfermi=',PF
158 c PF = 0.
159  ntry=ntry+1
160  IF(ntry.GT.2) WRITE(*,*) ntry,enu0,k2
161  IF (ntry .GT. 500) THEN
162  lbad = 1
163  lbad = 1
164  WRITE ( 6,1001) nbad, enu
165  RETURN
166  ENDIF
167 C CT = -1. + 2.*PYR(0)
168 c CT = -1.
169 C ST = SQRT(1.-CT*CT)
170 C F = 2.*3.1415926*PYR(0)
171 c F = 0.
172 
173 C P(2,4) = SQRT(PF*PF + MI2) - EBIND ! energia
174 C P(2,1) = PF*ST*COS(F) ! px
175 C P(2,2) = PF*ST*SIN(F) ! py
176 C P(2,3) = PF*CT ! pz
177 C P(2,5) = SQRT(P(2,4)**2-PF*PF) ! massa
178  p(2,1) = p21
179  p(2,2) = p22
180  p(2,3) = p23
181  p(2,4) = p24
182  p(2,5) = p25
183 c call lulist(1)
184  beta1=-p(2,1)/p(2,4)
185  beta2=-p(2,2)/p(2,4)
186  beta3=-p(2,3)/p(2,4)
187  n=2
188 C WRITE(6,*)' before transforming into target rest frame'
189 C call PYlist(1)
190  CALL pyrobo(0,0,0.,0.,beta1,beta2,beta3)
191 C print*,' nucl. rest fram ( fermi incl.) prima della rotazione'
192 C call PYlist(1)
193  n=5
194 
195  phi11=atan(p(1,2)/p(1,3))
196  pi(1)=p(1,1)
197  pi(2)=p(1,2)
198  pi(3)=p(1,3)
199 
200  CALL testrot1(pi,po,phi11)
201  DO ll=1,3
202  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
203  END DO
204 c WRITE(*,*) po
205  p(1,1)=po(1)
206  p(1,2)=po(2)
207  p(1,3)=po(3)
208  phi12=atan(p(1,1)/p(1,3))
209 
210  pi(1)=p(1,1)
211  pi(2)=p(1,2)
212  pi(3)=p(1,3)
213  CALL testrot2(pi,po,phi12)
214  DO ll=1,3
215  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
216  END DO
217 c WRITE(*,*) po
218  p(1,1)=po(1)
219  p(1,2)=po(2)
220  p(1,3)=po(3)
221 
222 
223 C call PYlist(1)
224 
225 
226  enu=p(1,4)
227 
228 C...Kinematical limits in Q**2
229 c S = P(2,5)**2 + 2.*ENU*(P(2,4)-P(2,3)) ! ????
230  s = p(2,5)**2 + 2.*enu*p(2,5)
231  sqs = sqrt(s) ! E centro massa
232  IF (sqs .LT. (aml + amf + 3.e-03)) goto 100
233  elf = (s-amf2+aml2)/(2.*sqs) ! energia leptone finale p
234  pstar = (s-p(2,5)**2)/(2.*sqs) ! p* neutrino nel c.m.
235  plf = sqrt(elf**2-aml2) ! 3-momento leptone finale
236  q2min = -aml2 + 2.*pstar*(elf-plf) ! + o -
237  q2max = -aml2 + 2.*pstar*(elf+plf) ! according con cos(theta)
238  IF (q2min .LT. 0.) q2min = 0. ! ??? non fisico
239 
240 C...Generate Q**2
241  dsigmax = dsqel_q2(ltyp,enu, q2min)
242  200 q2 = q2min + (q2max-q2min)*pyr(0)
243  dsig = dsqel_q2(ltyp,enu, q2)
244  IF (dsig .LT. dsigmax*pyr(0)) goto 200
245  CALL qgaus(q2min,q2max,dsigev,enu,ltyp)
246  ndsig=ndsig+1
247 C WRITE(6,*)' Q2,Q2min,Q2MAX,DSIGEV',
248 C &Q2,Q2min,Q2MAX,DSIGEV
249 
250 
251 C...c.m. frame. Neutrino along z axis
252  detot = (p(1,4)) + (p(2,4)) ! e totale
253  dbeta(1) = ((p(1,1)) + (p(2,1)))/detot ! px1+px2/etot = beta_x
254  dbeta(2) = ((p(1,2)) + (p(2,2)))/detot !
255  dbeta(3) = ((p(1,3)) + (p(2,3)))/detot !
256 c WRITE(*,*)
257 c WRITE(*,*)
258 C WRITE(*,*) 'Input values laboratory frame'
259 C CALL PYLIST(1)
260  n=2
261  CALL pyrobo(0,0,0.,0.,-dbeta(1),-dbeta(2),-dbeta(3))
262  n=5
263 c STHETA = ULANGL(P(1,3),P(1,1))
264 c write(*,*) 'stheta' ,stheta
265 c stheta=0.
266 c CALL PYROBO (0,0,-STHETA,0.,0.D0,0.D0,0.D0)
267 c WRITE(*,*)
268 c WRITE(*,*)
269 C WRITE(*,*) 'Output values cm frame'
270 C CALL PYLIST(1)
271 C...Kinematic in c.m. frame
272  ctstar = elf/plf - (q2 + aml2)/(2.*pstar*plf) ! cos(theta) cm
273  ststar = sqrt(1.-ctstar**2)
274  phi = 6.28319*pyr(0) ! random phi tra 0 e 2*pi
275  p(4,5) = aml ! massa leptone
276  p(4,4) = elf ! e leptone
277  p(4,3) = plf*ctstar ! px
278  p(4,1) = plf*ststar*cos(phi) ! py
279  p(4,2) = plf*ststar*sin(phi) ! pz
280 
281 
282  p(5,5) = amf ! barione
283  p(5,4) = (s+amf2-aml2)/(2.*sqs)! e barione
284  p(5,3) = -p(4,3) ! px
285  p(5,1) = -p(4,1) ! py
286  p(5,2) = -p(4,2) ! pz
287 
288 
289  p(3,5) = -q2
290  p(3,1) = p(1,1)-p(4,1)
291  p(3,2) = p(1,2)-p(4,2)
292  p(3,3) = p(1,3)-p(4,3)
293  p(3,4) = p(1,4)-p(4,4)
294 
295 C...Transform back to laboratory frame
296 c WRITE(*,*)
297 c WRITE(*,*)
298 C WRITE(*,*) 'before going back to nucl rest frame'
299 C CALL PYLIST(1)
300 c CALL PYROBO (0,0,STHETA,0.,0.D0,0.D0,0.D0)
301  n=5
302  CALL pyrobo(0,0,0.,0.,dbeta(1),dbeta(2),dbeta(3))
303 c WRITE(*,*)
304 c WRITE(*,*)
305 C WRITE(*,*) 'Now back in nucl rest frame'
306 C CALL PYLIST(1)
307  IF(ltyp.GE.3) CALL prepola(q2,ltyp,enu)
308 
309 c********************************************
310 
311  DO kw=1,5
312  pi(1)=p(kw,1)
313  pi(2)=p(kw,2)
314  pi(3)=p(kw,3)
315  CALL testrot3(pi,po,phi12)
316  DO ll=1,3
317  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
318  END DO
319  p(kw,1)=po(1)
320  p(kw,2)=po(2)
321  p(kw,3)=po(3)
322  END DO
323 c********************************************
324 
325 C call PYlist(1)
326 
327 
328 
329 
330 c********************************************
331 
332  DO kw=1,5
333  pi(1)=p(kw,1)
334  pi(2)=p(kw,2)
335  pi(3)=p(kw,3)
336  CALL testrot4(pi,po,phi11)
337  DO ll=1,3
338  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
339  END DO
340  p(kw,1)=po(1)
341  p(kw,2)=po(2)
342  p(kw,3)=po(3)
343  END DO
344 
345 c********************************************
346 
347 C WRITE(*,*) 'Now back in lab frame'
348 C call PYlist(1)
349  CALL pyrobo(1,5,0.,0.,-beta1,-beta2,-beta3)
350 CGB+...
351 C...test (on final momentum of nucleon) if Fermi-blocking
352 C...is operating
353  enucl = sqrt(p(5,1)**2 + p(5,2)**2 + p(5,3)**2 + p(5,5)**2)
354  & - p(5,5)
355  IF (enucl.LT. efmax) THEN
356  IF(inipri.LT.10)THEN
357  inipri=inipri+1
358  WRITE(6,*).LT.' qel: Pauli ENUCLEFMAX ', enucl,efmax
359 C...the interaction is not possible due to Pauli-Blocking and
360 C...it must be resampled
361  ENDIF
362  goto 100
363  ELSE IF (enucl.LT.enwell.and.enucl.GE.efmax) THEN
364  IF(inipri.LT.10)THEN
365  inipri=inipri+1
366  WRITE(6,*).LT.' qel: inside ENUCLENWELL ', enucl,enwell
367  ENDIF
368 C Reject (J:R) here all these events
369 C are otherwise rejected in dpmjet
370  goto 100
371 C...the interaction is possible, but the nucleon remains inside
372 C...the nucleus. The nucleus is therefore left excited.
373 C...We treat this case as a nucleon with 0 kinetic energy.
374 C P(5,5) = AMF
375 C P(5,4) = AMF
376 C P(5,1) = 0.
377 C P(5,2) = 0.
378 C P(5,3) = 0.
379  ELSE IF (enucl.GE.enwell) THEN
380 C WRITE(6,*)' qel ENUCL.GE.ENWELL ',ENUCL,ENWELL
381 C...the interaction is possible, the nucleon can exit the nucleus
382 C...but the nuclear well depth must be subtracted. The nucleus could be
383 C...left in an excited state.
384  pstart = sqrt(p(5,1)**2 + p(5,2)**2 + p(5,3)**2)
385 C P(5,4) = ENUCL-ENWELL + AMF
386  pnucl = sqrt(p(5,4)**2-amf**2)
387 C...The 3-momentum is scaled assuming that the direction remains
388 C...unaffected
389  p(5,1) = p(5,1) * pnucl/pstart
390  p(5,2) = p(5,2) * pnucl/pstart
391  p(5,3) = p(5,3) * pnucl/pstart
392 C WRITE(6,*)' qel new P(5,4) ',P(5,4)
393  ENDIF
394 CGB-...
395  dsigsu=dsigsu+dsigev
396 
397 C gbsumx = pxl + pxb + pxn
398 C gbsumy = pyl + pyb + pyn
399 C gbsumz = pzl + pzb + pzn
400 C gbsume = etl + etb + etn
401 C write(*,*) 'GB output to LEPDCY 1',ETL,PXL,PYL,PZL
402 C write(*,*) 'GB output to LEPDCY 2',ETb,PXb,PYb,PZb
403 C write(*,*) 'GB output to LEPDCY 3',ETn,PXn,PYn,PZn
404 C write(*,*) 'GB output to LEPDCY e',gbsumx,gbsumy,gbsumz,gbsume
405 C WRITE(6,*)' P(4,1-5) ',P(4,1),P(4,2),P(4,3),P(4,4),P(4,5)
406 C WRITE(6,*)' Q(4,1-5) ',Q(4,1),Q(4,2),Q(4,3),Q(4,4),Q(4,5)
407  ga=p(4,4)/p(4,5)
408  bgx=p(4,1)/p(4,5)
409  bgy=p(4,2)/p(4,5)
410  bgz=p(4,3)/p(4,5)
411 C WRITE(6,*)' g,bg1-3 ',GA,BGX,BGY,BGZ
412 C CALL DALTRA(GA,BGX,BGY,BGZ,PXL,PYL,PZL,ETL,
413 C & PPL,PPXL,PPYL,PPZL,EETL)
414 C CALL DALTRA(GA,BGX,BGY,BGZ,PXB,PYB,PZB,ETB,
415 C & PPB,PPXB,PPYB,PPZB,EETB)
416 C CALL DALTRA(GA,BGX,BGY,BGZ,PXN,PYN,PZN,ETN,
417 C & PPN,PPXN,PPYN,PPZN,EETN)
418 C rrsumx = ppxl + ppxb + ppxn
419 C rrsumy = ppyl + ppyb + ppyn
420 C rrsumz = ppzl + ppzb + ppzn
421 C rrsume = eetl + eetb + eetn
422 C write(*,*) 'rr output to LEPDCY 1',EETL,PPXL,PPYL,PPZL
423 C write(*,*) 'rr output to LEPDCY 2',EETb,PPXb,PPYb,PPZb
424 C write(*,*) 'rr output to LEPDCY 3',EETn,PPXn,PPYn,PPZn
425 C write(*,*) 'rr output to LEPDCY e',rrsumx,rrsumy,rrsumz,rrsume
426  dbetb(1)=bgx/ga
427  dbetb(2)=bgy/ga
428  dbetb(3)=bgz/ga
429  IF(neudec.EQ.1.OR.neudec.EQ.2) THEN
430  CALL pyrobo(6,8,0.,0.,dbetb(1),dbetb(2),dbetb(3))
431  ENDIF
432 c
433 C PRINT*,' FINE EVENTO '
434 C call PYlist(1)
435  enu=enu0
436  RETURN
437 
438  1001 FORMAT(2x, 'GEN_QEL : event rejected ', i5, g10.3)
439  END
440 
441 
442 C====================================================================
443 C. Masses
444 C====================================================================
445 
446  SUBROUTINE mass_ini
447 C...Initialize the kinematics for the quasi-elastic cross section
448  IMPLICIT DOUBLE PRECISION (a-h,o-z)
449  COMMON /cmaqel/ eml(6), emlsq(6), emn1(6), emn2(6), etqe(6)
450  & ,emn1sq(6),emn2sq(6),emprot,emneut,emn,emprotsq,emneutsq,emnsq
451  eml(1) = 0.51100d-03 ! e-
452  eml(2) = eml(1) ! e+
453  eml(3) = 0.105659d0 ! mu-
454  eml(4) = eml(3) ! mu+
455  eml(5) = 1.7777d0 ! tau-
456  eml(6) = eml(5) ! tau+
457  emprot = 0.93827231d0 ! p
458  emneut = 0.93956563d0 ! n
459  emprotsq = emprot**2
460  emneutsq = emneut**2
461  emn = (emprot + emneut)/2.
462  emnsq = emn**2
463  DO j=1,3
464  j0 = 2*(j-1)
465  emn1(j0+1) = emneut
466  emn1(j0+2) = emprot
467  emn2(j0+1) = emprot
468  emn2(j0+2) = emneut
469  ENDDO
470  DO j=1,6
471  emlsq(j) = eml(j)**2
472  etqe(j) = ((emn2(j)+ eml(j))**2-emn1(j)**2)/(2.*emn1(j))
473  ENDDO
474  RETURN
475  END
476 
477  FUNCTION dsqel_q2 (JTYP,ENU, Q2)
478 C...differential cross section for Quasi-Elastic scattering
479 C. nu + N -> l + N'
480 C. From Llewellin Smith Phys.Rep. 3C, 261, (1971).
481 C.
482 C. INPUT : JTYP = 1,...,6 nu_e, ...., nubar_tau
483 C. ENU (GeV) = Neutrino energy
484 C. Q2 (GeV**2) = (Transfer momentum)**2
485 C.
486 C. OUTPUT : DSQEL_Q2 = differential cross section :
487 C. dsigma/dq**2 (10**-38 cm+2/GeV**2)
488 C------------------------------------------------------------------
489  IMPLICIT DOUBLE PRECISION (a-h,o-z)
490  COMMON /caxial/ fa0, axial2
491  COMMON /cmaqel/ eml(6), emlsq(6), emn1(6), emn2(6), etqe(6)
492  & ,emn1sq(6),emn2sq(6),emprot,emneut,emn,emprotsq,emneutsq,emnsq
493  dimension ss(6)
494  DATA c0 /0.17590d0 / ! G_F**2 cos(theta_c)**2 M**2 /(8 pi) 10**-38 cm+2
495  DATA ss /1.d0, -1.d0, 1.d0, -1.d0, 1.d0, -1.d0/
496  DATA axial2 /1.03d0/ ! to be checked
497  fa0=-1.253d0
498  csi = 3.71d0 ! ???
499  gve = 1.d0/ (1.d0 + q2/0.84d0**2)**2 ! G_e(q**2)
500  gvm = (1.d0+csi)*gve ! G_m (q**2)
501  x = q2/(emn*emn) ! emn=massa barione
502  xa = x/4.d0
503  fv1 = 1.d0/(1.d0+xa)*(gve+xa*gvm)
504  fv2 = 1.d0/(1.d0+xa)*(gvm-gve)
505  fa = fa0/(1.d0 + q2/axial2)**2
506  ffa = fa*fa
507  ffv1 = fv1*fv1
508  ffv2 = fv2*fv2
509  rm = emlsq(jtyp)/(emn*emn) ! emlsq(jtyp)
510  a1 = (4.d0+x)*ffa - (4.d0-x)*ffv1 + x*ffv2*(1.d0-xa)+4*x*fv1*fv2
511  a2 = -rm * ((fv1 + fv2)**2 + ffa)
512  aa = (xa+0.25d0*rm)*(a1 + a2)
513  bb = -x*fa*(fv1 + fv2)
514  cc = 0.25d0*(ffa + ffv1 + xa*ffv2)
515  su = (4.d0*enu*emn - q2 - emlsq(jtyp))/(emn*emn)
516  dsqel_q2 = c0*(aa + ss(jtyp)*bb*su + cc*su*su) / (enu*enu) !
517  IF(dsqel_q2 .LT. 0.d0) dsqel_q2 = 0.d0
518  RETURN
519  END
520 
521  SUBROUTINE prepola(Q2,JTYP,ENU)
522  IMPLICIT DOUBLE PRECISION (a-h,o-z)
523 c
524 c By G. Battistoni and E. Scapparone (sept. 1997)
525 c According to:
526 c Albright & Jarlskog, Nucl Phys B84 (1975) 467
527 c
528 c
529 C COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
530  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
531 C REAL*4 P,V
532  COMMON /caxial/ fa0, axial2
533  COMMON /cmaqel/ eml(6), emlsq(6), emn1(6), emn2(6), etqe(6)
534  & ,emn1sq(6),emn2sq(6),emprot,emneut,emn,emprotsq,emneutsq,emnsq
535  REAL*8 mumass,pol(4,4),bb2(3),nu
536  common/pol/polarx(4),pmodul
537  COMMON /tautau/q(4,5),
538  + etl,pxl,pyl,pzl,
539  + etb,pxb,pyb,pzb,etn,pxn,pyn,pzn
540  COMMON /neutyy/neutyp,neudec
541  dimension ss(6)
542  DATA c0 /0.17590d0 / ! G_F**2 cos(theta_c)**2 M**2 /(8 pi) 10**-38 cm+2
543  DATA ss /1.d0, -1.d0, 1.d0, -1.d0, 1.d0, -1.d0/
544 C DATA AXIAL2 /1.03D0/ ! to be checked
545  rml=p(4,5)
546  rmm=0.93960d+00
547  fm2 = rmm**2
548  mpi = 0.135d+00
549  oldq2=q2
550  fa0=-1.253d+00
551  csi = 3.71d+00 !
552  gve = 1.d0/ (1.d0 + q2/(0.84d+00)**2)**2 ! G_e(q**2)
553  gvm = (1.d0+csi)*gve ! G_m (q**2)
554  x = q2/(emn*emn) ! emn=massa barione
555  xa = x/4.d0
556  fv1 = 1.d0/(1.d0+xa)*(gve+xa*gvm)
557  fv2 = 1.d0/(1.d0+xa)*(gvm-gve)
558  fa = fa0/(1.d0 + q2/axial2**2)**2
559  ffa = fa*fa
560  ffv1 = fv1*fv1
561  ffv2 = fv2*fv2
562  fp=2.d0*fa*rmm/(mpi**2 + q2)
563  rm = emlsq(jtyp)/(emn*emn) ! emlsq(jtyp)
564  a1 = (4.d0+x)*ffa-(4.d0-x)*ffv1+x*ffv2*(1.d0-xa)+4.d0*x*fv1*fv2
565  a2 = -rm * ((fv1 + fv2)**2 + ffa)
566  aa = (xa+0.25d+00*rm)*(a1 + a2)
567  bb = -x*fa*(fv1 + fv2)
568  cc = 0.25d+00*(ffa + ffv1 + xa*ffv2)
569  su = (4.d+00*enu*emn - q2 - emlsq(jtyp))/(emn*emn)
570 
571  omega1=ffa+xa*(ffa+(fv1+fv2)**2 ) ! articolo di ll...-smith
572  omega2=4.d+00*cc
573  omega3=2.d+00*fa*(fv1+fv2)
574  omega4p=(-(fv1+fv2)**2-(fa+2*fp)**2+(4.0d+00+
575  1 (q2/fm2))*fp**2)
576  omega5=omega2
577  omega4=(omega4p-omega2+2*omega5)/4.d+00
578  ww1=2.d+00*omega1*emn**2
579  ww2=2.d+00*omega2*emn**2
580  ww3=2.d+00*omega3*emn**2
581  ww4=2.d+00*omega4*emn**2
582  ww5=2.d+00*omega5*emn**2
583 
584  DO i=1,3
585  bb2(i)=-p(4,i)/p(4,4)
586  END DO
587 c WRITE(*,*)
588 c WRITE(*,*)
589 c WRITE(*,*) 'Prepola: ready to transform to lepton rest frame'
590 c CALL LULIST(1)
591  n=5
592  CALL pyrobo(0,0,0.,0.,bb2(1),bb2(2),bb2(3) )
593 * NOW PARTICLES ARE IN THE SCATTERED LEPTON REST FRAME
594 c WRITE(*,*)
595 c WRITE(*,*)
596 c WRITE(*,*) 'Prepola: now in lepton rest frame'
597 c CALL LULIST(1)
598  ee=enu
599  qm2=q2+rml**2
600  u=q2/(2.*rmm)
601  frac=qm2*ww1 + (2.d+00*ee*(ee-u) - 0.5d+00*qm2)*ww2 - ss(jtyp)*
602  + (0.5d+00/(rmm**2))*(2.d+00*rmm*ee*q2 - u*qm2)*ww3 +
603  + ((rml**2)/(2.d+00*fm2))*(qm2*ww4-2.d+00*rmm*ee*ww5) !<=FM2 inv di RMM!!
604 
605  factk=2.d+00*ww1 -ww2 -ss(jtyp)*(ee/rmm)*ww3 +((ee-u)/rmm)*ww5
606  + - ((rml**2)/fm2)*ww4 !<=FM2 inv di RMM!!
607 
608  factp=2.d+00*ee/rmm*ww2 - (qm2/(2.d+00*rmm**2))*(ss(jtyp)*ww3+ww5)
609 
610  DO i=1,3
611  pol(4,i)=rml*ss(jtyp)*(factk*p(1,i)+factp*p(2,i))/frac
612  polarx(i)=pol(4,i)
613  END DO
614 
615 
616  pmodul=0.d0
617  DO i=1,3
618  pmodul=pmodul+pol(4,i)**2
619  END DO
620 
621  IF(jtyp.GT.4.AND.neudec.GT.0) THEN
622  IF(neudec.EQ.1) THEN
623  CALL lepdcyp(eml(jtyp),eml(jtyp-2),polarx(3),
624  + etl,pxl,pyl,pzl,
625  + etb,pxb,pyb,pzb,etn,pxn,pyn,pzn)
626 c
627 c Tau has decayed in muon
628 c
629  ENDIF
630  IF(neudec.EQ.2) THEN
631  CALL lepdcyp(eml(jtyp),eml(jtyp-4),polarx(3),
632  + etl,pxl,pyl,pzl,
633  + etb,pxb,pyb,pzb,etn,pxn,pyn,pzn)
634 c
635 c Tau has decayed in electron
636 c
637  ENDIF
638  k(4,1)=15
639  k(4,4) = 6
640  k(4,5) = 8
641  n=n+3
642 c
643 c fill common for muon(electron)
644 c
645  p(6,1)=pxl
646  p(6,2)=pyl
647  p(6,3)=pzl
648  p(6,4)=etl
649  k(6,1)=1
650  IF(jtyp.EQ.5) THEN
651  IF(neudec.EQ.1) THEN
652  p(6,5)=eml(jtyp-2)
653  k(6,2)=13
654  ELSEIF(neudec.EQ.2) THEN
655  p(6,5)=eml(jtyp-4)
656  k(6,2)=11
657  ENDIF
658  ELSEIF(jtyp.EQ.6) THEN
659  IF(neudec.EQ.1) THEN
660  k(6,2)=-13
661  ELSEIF(neudec.EQ.2) THEN
662  k(6,2)=-11
663  ENDIF
664  END IF
665  k(6,3)=4
666  k(6,4)=0
667  k(6,5)=0
668 c
669 c fill common for tau_(anti)neutrino
670 c
671  p(7,1)=pxb
672  p(7,2)=pyb
673  p(7,3)=pzb
674  p(7,4)=etb
675  p(7,5)=0.
676  k(7,1)=1
677  IF(jtyp.EQ.5) THEN
678  k(7,2)=16
679  ELSEIF(jtyp.EQ.6) THEN
680  k(7,2)=-16
681  END IF
682  k(7,3)=4
683  k(7,4)=0
684  k(7,5)=0
685 c
686 c Fill common for muon(electron)_(anti)neutrino
687 c
688  p(8,1)=pxn
689  p(8,2)=pyn
690  p(8,3)=pzn
691  p(8,4)=etn
692  p(8,5)=0.
693  k(8,1)=1
694  IF(jtyp.EQ.5) THEN
695  IF(neudec.EQ.1) THEN
696  k(8,2)=-14
697  ELSEIF(neudec.EQ.2) THEN
698  k(8,2)=-12
699  ENDIF
700  ELSEIF(jtyp.EQ.6) THEN
701  IF(neudec.EQ.1) THEN
702  k(8,2)=14
703  ELSEIF(neudec.EQ.2) THEN
704  k(8,2)=12
705  ENDIF
706  END IF
707  k(8,3)=4
708  k(8,4)=0
709  k(8,5)=0
710  ENDIF
711 c WRITE(*,*)
712 c WRITE(*,*)
713 
714 c IF(PMODUL.GE.1.D+00) THEN
715 c WRITE(*,*) 'Pol',(POLARX(I),I=1,3)
716 c write(*,*) pmodul
717 c DO I=1,3
718 c POL(4,I)=POL(4,I)/PMODUL
719 c POLARX(I)=POL(4,I)
720 c END DO
721 c PMODUL=0.
722 c DO I=1,3
723 c PMODUL=PMODUL+POL(4,I)**2
724 c END DO
725 c WRITE(*,*) 'Pol',(POLARX(I),I=1,3)
726 c
727 c ENDIF
728 
729 c WRITE(*,*) 'PMODUL = ',PMODUL
730 
731 c WRITE(*,*)
732 c WRITE(*,*)
733 c WRITE(*,*) 'prepola: Now back to nucl rest frame'
734  CALL pyrobo(1,5,0.,0.,-bb2(1),-bb2(2),-bb2(3))
735 c call lulist(1)
736  xdc = v(4,1)+v(4,5)*p(4,1)/p(4,5)
737  ydc = v(4,2)+v(4,5)*p(4,2)/p(4,5)
738  zdc = v(4,3)+v(4,5)*p(4,3)/p(4,5)
739  DO ndc =6,8
740  v(ndc,1) = xdc
741  v(ndc,2) = ydc
742  v(ndc,3) = zdc
743  END DO
744 
745  RETURN
746  END
747 
748  SUBROUTINE testrot1(PI,PO,PHI)
749  IMPLICIT DOUBLE PRECISION (a-h,o-z)
750  dimension rot(3,3),pi(3),po(3)
751  rot(1,1)=1.d0
752  rot(1,2)=0.d0
753  rot(1,3)=0.d0
754  rot(2,1)=0.d0
755  rot(2,2)=cos(phi)
756  rot(2,3)=-sin(phi)
757  rot(3,1)=0.d0
758  rot(3,2)=sin(phi)
759  rot(3,3)=cos(phi)
760  DO 140 j=1,3
761  po(j)=rot(j,1)*pi(1)+rot(j,2)*pi(2)+rot(j,3)*pi(3)
762  140 CONTINUE
763  RETURN
764  END
765 
766 
767  SUBROUTINE testrot2(PI,PO,PHI)
768  IMPLICIT DOUBLE PRECISION (a-h,o-z)
769  dimension rot(3,3),pi(3),po(3)
770  rot(1,1)=0.d0
771  rot(1,2)=1.d0
772  rot(1,3)=0.d0
773  rot(2,1)=cos(phi)
774  rot(2,2)=0.d0
775  rot(2,3)=-sin(phi)
776  rot(3,1)=sin(phi)
777  rot(3,2)=0.d0
778  rot(3,3)=cos(phi)
779  DO 140 j=1,3
780  po(j)=rot(j,1)*pi(1)+rot(j,2)*pi(2)+rot(j,3)*pi(3)
781  140 CONTINUE
782  RETURN
783  END
784 
785  SUBROUTINE testrot3(PI,PO,PHI)
786  IMPLICIT DOUBLE PRECISION (a-h,o-z)
787  dimension rot(3,3),pi(3),po(3)
788  rot(1,1)=0.d0
789  rot(2,1)=1.d0
790  rot(3,1)=0.d0
791  rot(1,2)=cos(phi)
792  rot(2,2)=0.d0
793  rot(3,2)=-sin(phi)
794  rot(1,3)=sin(phi)
795  rot(2,3)=0.d0
796  rot(3,3)=cos(phi)
797  DO 140 j=1,3
798  po(j)=rot(j,1)*pi(1)+rot(j,2)*pi(2)+rot(j,3)*pi(3)
799  140 CONTINUE
800  RETURN
801  END
802 
803  SUBROUTINE testrot4(PI,PO,PHI)
804  IMPLICIT DOUBLE PRECISION (a-h,o-z)
805  dimension rot(3,3),pi(3),po(3)
806  rot(1,1)=1.d0
807  rot(2,1)=0.d0
808  rot(3,1)=0.d0
809  rot(1,2)=0.d0
810  rot(2,2)=cos(phi)
811  rot(3,2)=-sin(phi)
812  rot(1,3)=0.d0
813  rot(2,3)=sin(phi)
814  rot(3,3)=cos(phi)
815  DO 140 j=1,3
816  po(j)=rot(j,1)*pi(1)+rot(j,2)*pi(2)+rot(j,3)*pi(3)
817  140 CONTINUE
818  RETURN
819  END
820 
821  SUBROUTINE lepdcyp (AMA,AML,POL,ETL,PXL,PYL,PZL,
822  & etb,pxb,pyb,pzb,etn,pxn,pyn,pzn)
823 C
824 C-----------------------------------------------------------------
825 C
826 C Author :- G. Battistoni 10-NOV-1995
827 C
828 C=================================================================
829 C
830 C Purpose : performs decay of polarized lepton in
831 C its rest frame: a => b + l + anti-nu
832 C (Example: mu- => nu-mu + e- + anti-nu-e)
833 C Polarization is assumed along Z-axis
834 C WARNING:
835 C 1) b AND anti-nu ARE ASSUMED TO BE NEUTRINOS
836 C OF NEGLIGIBLE MASS
837 C 2) RADIATIVE CORRECTIONS ARE NOT CONSIDERED
838 C IN THIS VERSION
839 C
840 C Method : modifies phase space distribution obtained
841 C by routine EXPLOD using a rejection against the
842 C matrix element for unpolarized lepton decay
843 C
844 C Inputs : Mass of a : AMA
845 C Mass of l : AML
846 C Polar. of a: POL
847 C (Example: fully polar. mu- decay: AMA=AMMUON, AML=AMELCT,
848 C POL = -1)
849 C
850 C Outputs : kinematic variables in the rest frame of decaying lepton
851 C ETL,PXL,PYL,PZL 4-moment of l
852 C ETB,PXB,PYB,PZB 4-moment of b
853 C ETN,PXN,PYN,PZN 4-moment of anti-nu
854 C
855 C============================================================
856 C +
857 C Declarations.
858 C -
859 *
860 *=== dblprc ==========================================================*
861 * *
862  IMPLICIT DOUBLE PRECISION (a-h,o-z)
863  parameter( kalgnm = 2 )
864  parameter( anglgb = 5.0d-16 )
865  parameter( anglsq = 2.5d-31 )
866  parameter( axcssv = 0.2d+16 )
867  parameter( andrfl = 1.0d-38 )
868  parameter( avrflw = 1.0d+38 )
869  parameter( ainfnt = 1.0d+30 )
870  parameter( azrzrz = 1.0d-30 )
871  parameter( einfnt = +69.07755278982137 d+00 )
872  parameter( ezrzrz = -69.07755278982137 d+00 )
873  parameter( onemns = 0.999999999999999 d+00 )
874  parameter( onepls = 1.000000000000001 d+00 )
875  parameter( csnnrm = 2.0d-15 )
876  parameter( dmxtrn = 1.0d+08 )
877  parameter( zerzer = 0.d+00 )
878  parameter( oneone = 1.d+00 )
879  parameter( twotwo = 2.d+00 )
880  parameter( thrthr = 3.d+00 )
881  parameter( foufou = 4.d+00 )
882  parameter( fivfiv = 5.d+00 )
883  parameter( sixsix = 6.d+00 )
884  parameter( sevsev = 7.d+00 )
885  parameter( eigeig = 8.d+00 )
886  parameter( aninen = 9.d+00 )
887  parameter( tenten = 10.d+00 )
888  parameter( hlfhlf = 0.5d+00 )
889  parameter( onethi = oneone / thrthr )
890  parameter( twothi = twotwo / thrthr )
891  parameter( pipipi = 3.1415926535897932270 d+00 )
892  parameter( eneper = 2.7182818284590452354 d+00 )
893  parameter( sqrent = 1.6487212707001281468 d+00 )
894  parameter( clight = 2.99792458 d+10 )
895  parameter( avogad = 6.0221367 d+23 )
896  parameter( amelgr = 9.1093897 d-28 )
897  parameter( plckbr = 1.05457266 d-27 )
898  parameter( elccgs = 4.8032068 d-10 )
899  parameter( elcmks = 1.60217733 d-19 )
900  parameter( amugrm = 1.6605402 d-24 )
901  parameter( ammumu = 0.113428913 d+00 )
902  parameter( alpfsc = 7.2973530791728595 d-03 )
903  parameter( fscto2 = 5.3251361962113614 d-05 )
904  parameter( fscto3 = 3.8859399018437826 d-07 )
905  parameter( fscto4 = 2.8357075508200407 d-09 )
906  parameter( plabrc = 0.197327053 d+00 )
907  parameter( amelct = 0.51099906 d-03 )
908  parameter( amugev = 0.93149432 d+00 )
909  parameter( ammuon = 0.105658389 d+00 )
910  parameter( rclsel = 2.8179409183694872 d-13 )
911  parameter( gevmev = 1.0 d+03 )
912  parameter( emvgev = 1.0 d-03 )
913  parameter( algvmv = 6.90775527898214 d+00 )
914  parameter( raddeg = 180.d+00 / pipipi )
915  parameter( degrad = pipipi / 180.d+00 )
916 C +
917 C variables for EXPLOD
918 C -
919  parameter( kpmx = 10 )
920  dimension amexpl(kpmx), pxexpl(kpmx), pyexpl(kpmx),
921  & pzexpl(kpmx), ptexpl(kpmx), etexpl(kpmx),
922  & amhelp(kpmx)
923 C +
924 C test variables
925 C -
926  COMMON /gbatnu/ elerat,ntry
927 C +
928 C Initializes test variables
929 C -
930  ntry = 0
931  elerat = 0.d+00
932 C +
933 C Maximum value for matrix element
934 C -
935  elemax = ( ama**2 + aml**2 )**2 / ama**2 * ( ama**2 - aml**2 +
936  & sqrt( ama**4 + aml**4 - 3.d+00 * ama**2 * aml**2 ) )
937 C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
938 C Inputs for EXPLOD
939 C part. no. 1 is l (e- in mu- decay)
940 C part. no. 2 is b (nu-mu in mu- decay)
941 C part. no. 3 is anti-nu (anti-nu-e in mu- decay)
942 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
943  npexpl = 3
944  etotex = ama
945  amexpl(1) = aml
946  amexpl(2) = 0.d+00
947  amexpl(3) = 0.d+00
948 C +
949 C phase space distribution
950 C -
951  100 CONTINUE
952  ntry = ntry + 1
953  CALL explod( npexpl, amexpl, etotex, etexpl, pxexpl,
954  & pyexpl, pzexpl )
955 C +
956 C Calculates matrix element:
957 C 64*GF**2{[P(a)-ama*S(a)]*P(anti-nu)}{P(l)*P(b)}
958 C Here CTH is the cosine of the angle between anti-nu and Z axis
959 C -
960  cth = pzexpl(3) / sqrt( pxexpl(3)**2 + pyexpl(3)**2 +
961  & pzexpl(3)**2 )
962  prod1 = etexpl(3) * ama * (1.d+00 - pol * cth)
963  prod2 = etexpl(1) * etexpl(2) - pxexpl(1)*pxexpl(2) -
964  & pyexpl(1)*pyexpl(2) - pzexpl(1)*pzexpl(2)
965  elemat = 16.d+00 * prod1 * prod2
966  IF(elemat.GT.elemax) THEN
967  WRITE(6,*) 'Problems in LEPDCY',elemax,elemat
968  stop
969  ENDIF
970 C +
971 C Here performs the rejection
972 C -
973  test = rndm(etotex) * elemax
974  IF ( test .GT. elemat ) go to 100
975 C +
976 C final assignment of variables
977 C -
978  elerat = elemat/elemax
979  etl = etexpl(1)
980  pxl = pxexpl(1)
981  pyl = pyexpl(1)
982  pzl = pzexpl(1)
983  etb = etexpl(2)
984  pxb = pxexpl(2)
985  pyb = pyexpl(2)
986  pzb = pzexpl(2)
987  etn = etexpl(3)
988  pxn = pxexpl(3)
989  pyn = pyexpl(3)
990  pzn = pzexpl(3)
991  999 RETURN
992  END
993 
994 C==================================================================
995 C. Generation of Delta resonance events
996 C==================================================================
997 
998  SUBROUTINE gen_delta(ENU,LLEP,LTARG,JINT,P21,P22,P23,P24,P25)
999  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1000 C...Generate a Delta-production neutrino/antineutrino
1001 C. CC-interaction on a nucleon
1002 C
1003 C. INPUT ENU (GeV) = Neutrino Energy
1004 C. LLEP = neutrino type
1005 C. LTARG = nucleon target type 1=p, 2=n.
1006 C. JINT = 1:CC, 2::NC
1007 C.
1008 C. OUTPUT PPL(4) 4-monentum of final lepton
1009 C----------------------------------------------------
1010 C COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
1011  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
1012  COMMON /kinqel/ ppl(5), ppn(5)
1013  COMMON /cbad/ lbad, nbad
1014  COMMON /cconlun/ luna
1015  dimension rot(3,3),pi(3),po(3)
1016 C REAL*4 AMD0, AMD, AMN(2), AML0(6), AML, AML2, AMDMIN
1017  dimension aml0(6),amn(2)
1018  DATA amd0 /1.231/, gamd /0.12/, deld/0.169/, amdmin/1.084/
1019  DATA amn /0.93827231, 0.93956563/
1020  DATA aml0 /2*0.51100e-03,2*0.105659, 2*1.777/
1021 
1022 c WRITE(6,*)' GEN_DEL',ENU,LLEP,LTARG,JINT,P21,P22,P23,P24,P25
1023  lbad = 0
1024 C...Final lepton mass
1025  IF (jint.EQ.1) THEN
1026  aml = aml0(llep)
1027  ELSE
1028  aml = 0.
1029  ENDIF
1030  aml2 = aml**2
1031 
1032 C...Particle labels (LUND)
1033  n = 5
1034  k(1,1) = 21
1035  k(2,1) = 21
1036  k(3,1) = 21
1037  k(4,1) = 1
1038  k(3,3) = 1
1039  k(4,3) = 1
1040  IF (ltarg .EQ. 1) THEN
1041  k(2,2) = 2212
1042  ELSE
1043  k(2,2) = 2112
1044  ENDIF
1045  k0 = (llep-1)/2
1046  k1 = llep/2
1047  ka = 12 + 2*k0
1048  is = -1 + 2*llep - 4*k1
1049  lnu = 2 - llep + 2*k1
1050  k(1,2) = is*ka
1051  k(5,1) = 1
1052  k(5,3) = 2
1053  IF (jint .EQ. 1) THEN ! CC interactions
1054  k(3,2) = is*24
1055  k(4,2) = is*(ka-1)
1056  IF(lnu.EQ.1) THEN
1057  IF (ltarg .EQ. 1) THEN
1058  k(5,2) = 2224
1059  ELSE
1060  k(5,2) = 2214
1061  ENDIF
1062  ELSE
1063  IF (ltarg .EQ. 1) THEN
1064  k(5,2) = 2114
1065  ELSE
1066  k(5,2) = 1114
1067  ENDIF
1068  ENDIF
1069  ELSE
1070  k(3,2) = 23 ! NC (Z0) interactions
1071  k(4,2) = k(1,2)
1072  IF (ltarg .EQ. 1) THEN
1073  k(5,2) = 2114
1074  ELSE
1075  k(5,2) = 2214
1076  ENDIF
1077  ENDIF
1078 
1079 C...4-momentum initial lepton
1080  p(1,5) = 0.
1081  p(1,4) = enu
1082  p(1,1) = 0.
1083  p(1,2) = 0.
1084  p(1,3) = enu
1085 C...4-momentum initial nucleon
1086  p(2,5) = amn(ltarg)
1087 C P(2,4) = P(2,5)
1088 C P(2,1) = 0.
1089 C P(2,2) = 0.
1090 C P(2,3) = 0.
1091  p(2,1) = p21
1092  p(2,2) = p22
1093  p(2,3) = p23
1094  p(2,4) = p24
1095  p(2,5) = p25
1096  n=2
1097 C call PYlist(1)
1098  beta1=-p(2,1)/p(2,4)
1099  beta2=-p(2,2)/p(2,4)
1100  beta3=-p(2,3)/p(2,4)
1101  n=2
1102  CALL pyrobo(0,0,0.,0.,beta1,beta2,beta3)
1103 
1104 C print*,' nucl. rest fram ( fermi incl.) prima della rotazione'
1105 C call PYlist(1)
1106 
1107  phi11=atan(p(1,2)/p(1,3))
1108  pi(1)=p(1,1)
1109  pi(2)=p(1,2)
1110  pi(3)=p(1,3)
1111 
1112  CALL testrot1(pi,po,phi11)
1113  DO ll=1,3
1114  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
1115  END DO
1116  p(1,1)=po(1)
1117  p(1,2)=po(2)
1118  p(1,3)=po(3)
1119  phi12=atan(p(1,1)/p(1,3))
1120 
1121  pi(1)=p(1,1)
1122  pi(2)=p(1,2)
1123  pi(3)=p(1,3)
1124  CALL testrot2(pi,po,phi12)
1125  DO ll=1,3
1126  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
1127  END DO
1128  p(1,1)=po(1)
1129  p(1,2)=po(2)
1130  p(1,3)=po(3)
1131 C call PYlist(1)
1132 
1133  enuu=p(1,4)
1134 
1135 C...Generate the Mass of the Delta
1136  ntry = 0
1137 100 r = pyr(0)
1138  amd=amd0+0.5*gamd*tan((2.*r-1.)*atan(2.*deld/gamd))
1139  ntry = ntry + 1
1140  IF (ntry .GT. 1000) THEN
1141  lbad = 1
1142  WRITE ( 6,1001) nbad, enuu,amd,amdmin,amd0,gamd,et
1143  WRITE (luna,1001) nbad, enuu
1144  RETURN
1145  ENDIF
1146  IF (amd .LT. amdmin) goto 100
1147  et = ((amd+aml)**2 - amn(ltarg)**2)/(2.*amn(ltarg))
1148  IF (enuu .LT. et) goto 100
1149 
1150 C...Kinematical limits in Q**2
1151  s = amn(ltarg)**2 + 2.*amn(ltarg)*enuu
1152  sqs = sqrt(s)
1153  pstar = (s - amn(ltarg)**2)/(2.*sqs)
1154  elf = (s - amd**2 + aml2)/(2.*sqs)
1155  plf = sqrt(elf**2 - aml2)
1156  q2min = -aml2 + 2.*pstar*(elf-plf)
1157  q2max = -aml2 + 2.*pstar*(elf+plf)
1158  IF (q2min .LT. 0.) q2min = 0.
1159 
1160  dsigmax = dsigma_delta(lnu,-q2min, s, aml, amd)
1161 200 q2 = q2min + (q2max-q2min)*pyr(0)
1162  dsig = dsigma_delta(lnu,-q2, s, aml, amd)
1163  IF (dsig .LT. dsigmax*pyr(0)) goto 200
1164 
1165 
1166 C...Generate the kinematics of the final particles
1167  eistar = (s + amn(ltarg)**2)/(2.*sqs)
1168  gam = eistar/amn(ltarg)
1169  bet = pstar/eistar
1170  ctstar = elf/plf - (q2 + aml2)/(2.*pstar*plf)
1171  el = gam*(elf + bet*plf*ctstar)
1172  plz = gam*(plf*ctstar + bet*elf)
1173  pl = sqrt(el**2 - aml2)
1174  plt = sqrt(max(1.e-06,(pl*pl - plz*plz)))
1175  phi = 6.28319*pyr(0)
1176  p(4,1) = plt*cos(phi)
1177  p(4,2) = plt*sin(phi)
1178  p(4,3) = plz
1179  p(4,4) = el
1180  p(4,5) = aml
1181 
1182 C...4-momentum of Delta
1183  p(5,1) = -p(4,1)
1184  p(5,2) = -p(4,2)
1185  p(5,3) = enuu-p(4,3)
1186  p(5,4) = enuu+amn(ltarg)-p(4,4)
1187  p(5,5) = amd
1188 
1189 C...4-momentum of intermediate boson
1190  p(3,5) = -q2
1191  p(3,4) = p(1,4)-p(4,4)
1192  p(3,1) = p(1,1)-p(4,1)
1193  p(3,2) = p(1,2)-p(4,2)
1194  p(3,3) = p(1,3)-p(4,3)
1195  n=5
1196 C call PYlist(1)
1197 
1198  DO kw=1,5
1199  pi(1)=p(kw,1)
1200  pi(2)=p(kw,2)
1201  pi(3)=p(kw,3)
1202  CALL testrot3(pi,po,phi12)
1203  DO ll=1,3
1204  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
1205  END DO
1206  p(kw,1)=po(1)
1207  p(kw,2)=po(2)
1208  p(kw,3)=po(3)
1209  END DO
1210 
1211 c********************************************
1212 C call PYlist(1)
1213 
1214 c********************************************
1215 
1216  DO kw=1,5
1217  pi(1)=p(kw,1)
1218  pi(2)=p(kw,2)
1219  pi(3)=p(kw,3)
1220  CALL testrot4(pi,po,phi11)
1221  DO ll=1,3
1222  IF(abs(po(ll)).LT.1.d-07) po(ll)=0.
1223  END DO
1224  p(kw,1)=po(1)
1225  p(kw,2)=po(2)
1226  p(kw,3)=po(3)
1227  END DO
1228 c********************************************
1229 C call PYlist(1)
1230 
1231 C transform back into Lab.
1232  CALL pyrobo(0,0,0.,0.,-beta1,-beta2,-beta3)
1233 
1234 C WRITE(6,*)' Lab fram ( fermi incl.) '
1235 C call PYlist(1)
1236  n=5
1237  CALL pyexec
1238 C call PYlist(1)
1239 
1240 
1241  RETURN
1242 1001 FORMAT(2x, 'GEN_DELTA : event rejected ', i5, 6g10.3)
1243  END
1244 
1245  FUNCTION dsigma_delta (LNU, QQ, S, AML, MD)
1246  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1247 C...Reaction nu + N -> lepton + Delta
1248 C. returns the cross section
1249 C. dsigma/dt
1250 C. INPUT LNU = 1, 2 (neutrino-antineutrino)
1251 C. QQ = t (always negative) GeV**2
1252 C. S = (c.m energy)**2 GeV**2
1253 C. OUTPUT = 10**-38 cm+2/GeV**2
1254 C-----------------------------------------------------
1255  REAL*8 mn, mn2, mn4, md,md2, md4
1256  DATA mn /0.938/
1257  DATA pi /3.1415926/
1258  gf = (1.1664 * 1.97)
1259  gf2 = gf*gf
1260  mn2 = mn*mn
1261  mn4 = mn2*mn2
1262  md2 = md*md
1263  md4 = md2*md2
1264  aml2 = aml*aml
1265  aml4 = aml2*aml2
1266  vq = (mn2 - md2 - qq)/2.
1267  vpi = (mn2 + md2 - qq)/2.
1268  vk = (s + qq - mn2 - aml2)/2.
1269  pik = (s - mn2)/2.
1270  qk = (aml2 - qq)/2.
1271  piq = (qq + mn2 - md2)/2.
1272  q = sqrt(-qq)
1273  c3v = 2.07*sqrt(exp(-6.3*q)*(1.+9*q))
1274  c3 = sqrt(3.)*c3v/mn
1275  c4 = -c3/md ! attenzione al segno
1276  c5a = 1.18/(1.-qq/0.4225)**2
1277  c32 = c3**2
1278  c42 = c4**2
1279  c5a2 = c5a**2
1280 
1281  IF (lnu .EQ. 1) THEN
1282  ans3=-md2*vpi*qk*qq*c32+md2*vpi*qk*c5a2+2.*md2*vq*
1283  . pik*qk*c32+2.*md2*vq*qk*piq*c32+md4*vpi*qk*qq*c42-
1284  . 2.*vk**2*vpi*qq*c32+2.*vk**2*vpi*c5a2+4.*vk*vpi*vq*
1285  . qk*c32+2.*vk*vpi*vq*c5a2+2.*vpi*vq**2*qk*c32
1286  ans2=2.*mn*md*md2*vk**2*qq*c42-4.*mn*md*md2*vk*vq*qk
1287  . *c42-2.*mn*md*md2*vq**2*qk*c42-2.*mn*md*md2*qk**2*
1288  . c32-3.*mn*md*md2*qk*qq*c32+mn*md*md2*qk*c5a2-mn*md*
1289  . md4*qk*qq*c42+2.*mn*md*vk**2*c5a2+2.*mn*md*vk*vq*
1290  . c5a2+4.*mn*c3*c4*md2*vk**2*qq-8.*mn*c3*c4*md2*vk*vq
1291  . *qk-4.*mn*c3*c4*md2*vq**2*qk-2.*mn*c3*c4*md4*qk*qq-
1292  . 4.*mn*c3*c5a*md2*vk*qq+4.*mn*c3*c5a*md2*vq*qk-2.*md*
1293  . c3*c4*md2*vk*pik*qq+2.*md*c3*c4*md2*vk*qk*piq+2.*md
1294  . *c3*c4*md2*vpi*qk*qq+2.*md*c3*c4*md2*vq*pik*qk+2.*
1295  . md*c3*c4*md2*vq*qk*piq-2.*md*c3*c4*vk**2*vpi*qq+4.*
1296  . md*c3*c4*vk*vpi*vq*qk+2.*md*c3*c4*vpi*vq**2*qk-md*
1297  . c3*c5a*md2*pik*qq+md*c3*c5a*md2*qk*piq-3.*md*c3*c5a
1298  . *vk*vpi*qq+md*c3*c5a*vk*vq*piq+3.*md*c3*c5a*vpi*vq*
1299  . qk-md*c3*c5a*vq**2*pik+c4*c5a*md2*vk*vpi*qq+c4*c5a*
1300  . md2*vk*vq*piq-c4*c5a*md2*vpi*vq*qk-c4*c5a*md2*vq**2
1301  . *pik-c4*c5a*md4*pik*qq+c4*c5a*md4*qk*piq-2.*md2*vk
1302  . **2*vpi*qq*c42+4.*md2*vk*vpi*vq*qk*c42-2.*md2*vk*
1303  . pik*qq*c32+2.*md2*vk*qk*piq*c32+2.*md2*vpi*vq**2*qk
1304  . *c42-2.*md2*vpi*qk**2*c32+ans3
1305  ELSE
1306  ans3=-md2*vpi*qk*qq*c32+md2*vpi*qk*c5a2+2.*md2*vq*
1307  . pik*qk*c32+2.*md2*vq*qk*piq*c32+md4*vpi*qk*qq*c42-
1308  . 2.*vk**2*vpi*qq*c32+2.*vk**2*vpi*c5a2+4.*vk*vpi*vq*
1309  . qk*c32+2.*vk*vpi*vq*c5a2+2.*vpi*vq**2*qk*c32
1310  ans2=2.*mn*md*md2*vk**2*qq*c42-4.*mn*md*md2*vk*vq*qk
1311  . *c42-2.*mn*md*md2*vq**2*qk*c42-2.*mn*md*md2*qk**2*
1312  . c32-3.*mn*md*md2*qk*qq*c32+mn*md*md2*qk*c5a2-mn*md*
1313  . md4*qk*qq*c42+2.*mn*md*vk**2*c5a2+2.*mn*md*vk*vq*
1314  . c5a2+4.*mn*c3*c4*md2*vk**2*qq-8.*mn*c3*c4*md2*vk*vq
1315  . *qk-4.*mn*c3*c4*md2*vq**2*qk-2.*mn*c3*c4*md4*qk*qq+
1316  . 4.*mn*c3*c5a*md2*vk*qq-4.*mn*c3*c5a*md2*vq*qk-2.*md*
1317  . c3*c4*md2*vk*pik*qq+2.*md*c3*c4*md2*vk*qk*piq+2.*md
1318  . *c3*c4*md2*vpi*qk*qq+2.*md*c3*c4*md2*vq*pik*qk+2.*
1319  . md*c3*c4*md2*vq*qk*piq-2.*md*c3*c4*vk**2*vpi*qq+4.*
1320  . md*c3*c4*vk*vpi*vq*qk+2.*md*c3*c4*vpi*vq**2*qk+md*
1321  . c3*c5a*md2*pik*qq-md*c3*c5a*md2*qk*piq+3.*md*c3*c5a
1322  . *vk*vpi*qq-md*c3*c5a*vk*vq*piq-3.*md*c3*c5a*vpi*vq*
1323  . qk+md*c3*c5a*vq**2*pik-c4*c5a*md2*vk*vpi*qq-c4*c5a*
1324  . md2*vk*vq*piq+c4*c5a*md2*vpi*vq*qk+c4*c5a*md2*vq**2
1325  . *pik+c4*c5a*md4*pik*qq-c4*c5a*md4*qk*piq-2.*md2*vk
1326  . **2*vpi*qq*c42+4.*md2*vk*vpi*vq*qk*c42-2.*md2*vk*
1327  . pik*qq*c32+2.*md2*vk*qk*piq*c32+2.*md2*vpi*vq**2*qk
1328  . *c42-2.*md2*vpi*qk**2*c32+ans3
1329  ENDIF
1330  ans1=32.*ans2
1331  ans=ans1/(3.*md2)
1332  p1cm = (s-mn2)/(2.*sqrt(s))
1333  dsigma_delta = gf2/2. * ans/(64.*pi*s*p1cm**2)
1334  RETURN
1335  END
1336 
1337  SUBROUTINE qgaus(A,B,SS,ENU,LTYP)
1338  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1339  dimension x(5),w(5)
1340  DATA x/.1488743389d0,.4333953941d0,
1341  & .6794095682d0,.8650633666d0,.9739065285d0
1342  */
1343  DATA w/.2955242247d0,.2692667193d0,
1344  & .2190863625d0,.1494513491d0,.0666713443d0
1345  */
1346  xm=0.5d0*(b+a)
1347  xr=0.5d0*(b-a)
1348  ss=0
1349  DO 11 j=1,5
1350  dx=xr*x(j)
1351  ss=ss+w(j)*(dsqel_q2(ltyp,enu,xm+dx)+
1352  * dsqel_q2(ltyp,enu,xm-dx))
1353 11 CONTINUE
1354  ss=xr*ss
1355  RETURN
1356  END
pyexec
subroutine pyexec
Definition: pythia61.f:31539
prepola
subroutine prepola(Q2, JTYP, ENU)
Definition: dpm25qelpo.f:521
gen_qel
subroutine gen_qel(ENU, LTYP, P21, P22, P23, P24, P25)
Definition: dpm25nonu.f:7
s
const XML_Char * s
Definition: include/expat.h:262
z
G4double z
Definition: TRTMaterials.hh:39
p
const char * p
Definition: xmltok.h:285
dsigma_delta
function dsigma_delta(LNU, QQ, S, AML, MD)
Definition: dpm25qelpo.f:1245
HepGeom::Transform3D::dx
double dx() const
Definition: Transform3D.h:279
c2_factory::tan
static c2_tan_p< float_type > & tan()
make a *new object
Definition: c2_factory.hh:136
a
G4double a
Definition: TRTMaterials.hh:39
HepGeom::Plane3D::d
T d() const
Definition: Plane3D.h:86
nmxhkk
const int nmxhkk
Definition: G4DPMJET2_5Interface.hh:78
qgaus
subroutine qgaus(A, B, SS, ENU, LTYP)
Definition: dpm25qelpo.f:1337
ebind
double precision function ebind(IA, IZ)
Definition: dpm25nuc1.f:5789
dbeta
double precision function dbeta(X1, X2, BET)
Definition: dpm25nuc7.f:2672
qel_pol
subroutine qel_pol(ENU, LTYP, P21, P22, P23, P24, P25)
Definition: dpm25nonu.f:1
gen_delta
subroutine gen_delta(ENU, LLEP, LTARG, JINT, P21, P22, P23, P24, P25)
Definition: dpm25nonu.f:19
pyjets
#define pyjets
Definition: pythia6_common_address.c:39
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
rndm
double precision function rndm(RDUMMY)
Definition: dpm25nulib.f:1460
testrot3
subroutine testrot3(PI, PO, PHI)
Definition: dpm25qelpo.f:785
testrot4
subroutine testrot4(PI, PO, PHI)
Definition: dpm25qelpo.f:803
lepdcyp
subroutine lepdcyp(AMA, AML, POL, ETL, PXL, PYL, PZL, ETB, PXB, PYB, PZB, ETN, PXN, PYN, PZN)
Definition: dpm25qelpo.f:821
CLHEP::HepLorentzVector::et
double et() const
c2_factory::sqrt
static c2_sqrt_p< float_type > & sqrt()
make a *new object
Definition: c2_factory.hh:142
dsqel_q2
function dsqel_q2(JTYP, ENU, Q2)
Definition: dpm25qelpo.f:477
testrot2
subroutine testrot2(PI, PO, PHI)
Definition: dpm25qelpo.f:767
c2_factory::cos
static c2_cos_p< float_type > & cos()
make a *new object
Definition: c2_factory.hh:134
test
program test
Definition: Main_HIJING.f:1
testrot1
subroutine testrot1(PI, PO, PHI)
Definition: dpm25qelpo.f:748
pyrobo
subroutine pyrobo(IMI, IMA, THE, PHI, BEX, BEY, BEZ)
Definition: pythia61.f:36748
mass_ini
subroutine mass_ini
Definition: dpm25nonu.f:14
pyr
function pyr(IDUMMY)
Definition: pythia61.f:36587
c2_factory::sin
static c2_sin_p< float_type > & sin()
make a *new object
Definition: c2_factory.hh:132
c2_factory::exp
static c2_exp_p< float_type > & exp()
make a *new object
Definition: c2_factory.hh:140
Generated on Wed Apr 30 2014 15:54:32 for Geant4.10 by   doxygen 1.8.7